Plot constant surface on tricontourf in Python

Question

I currently have a tricontourf contour map and I wish to interpolate the points and plot a line on the contour plot where the z value is equal to 0. Currently the plot looks like this:

And the accompanying code looks like this:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import scipy.interpolate

# Load the 3D data file
data = np.genfromtxt("Ta_parameterspace_2mm.txt", skip_header=14, delimiter="\t", dtype = float)
#print(data)
reflect = data[:,0]
emiss = data[:,1]
tempdiff = data[:,4]

fig, ax = plt.subplots()
cb = ax.tricontourf(reflect, emiss, tempdiff,100, cmap = "seismic")
cbar = plt.colorbar(cb)
cbar.set_label('Temperature (K)', rotation = 270, labelpad = 13)
ax.set_xlabel('Reflectivity')
ax.set_ylabel('Emissivity')
plt.savefig('Ta_parameterplot_diff.pdf', bbox_inches='tight', format='pdf')
plt.savefig('Ta_parameterplot_diff.png', dpi=300, bbox_inches='tight', format='png')
plt.show()

As mentioned at the begining of the question, I would like to interpolate between the points in the dataset to allow us to plot a line on the contour where 'Temperature' is equal to zero. How can I do this? The full data file can be found below:

https://fz-juelich.sciebo.de/s/SjwZyAPB4oEerZE

score 0 · Answer 1 · answered Mar 01 '23 at 08:57

Actually, what you want is the coordinates of intersection point of z-value and zero-values.

I borrowed code from this answer : Find the intersection of two curves given by (x, y) data with high precision in Python

The three arguments in the function interpolated_intercept in your case will be:

x --> y_node (which is the grid coordinated y)
y1 --> z_axis_y (z values along the y-axis)
y2 --> z_zero (an array of zeros has the same size with z_axis_y)

Here is the full code :

import matplotlib.pyplot as plt
import numpy as np


def interpolated_intercepts(x, y1, y2):
    """Find the intercepts of two curves, given by the same x data"""

    def intercept(point1, point2, point3, point4):
        """find the intersection between two lines
        the first line is defined by the line between point1 and point2
        the first line is defined by the line between point3 and point4
        each point is an (x,y) tuple.

        So, for example, you can find the intersection between
        intercept((0,0), (1,1), (0,1), (1,0)) = (0.5, 0.5)

        Returns: the intercept, in (x,y) format
        """    

        def line(p1, p2):
            A = (p1[1] - p2[1])
            B = (p2[0] - p1[0])
            C = (p1[0]*p2[1] - p2[0]*p1[1])
            return A, B, -C

        def intersection(L1, L2):
            D  = L1[0] * L2[1] - L1[1] * L2[0]
            Dx = L1[2] * L2[1] - L1[1] * L2[2]
            Dy = L1[0] * L2[2] - L1[2] * L2[0]

            x = Dx / D
            y = Dy / D
            return x,y

        L1 = line([point1[0],point1[1]], [point2[0],point2[1]])
        L2 = line([point3[0],point3[1]], [point4[0],point4[1]])

        R = intersection(L1, L2)

        return R

    idxs = np.argwhere(np.diff(np.sign(y1 - y2)) != 0)

    xcs = []
    ycs = []

    for idx in idxs:
        xc, yc = intercept((x[idx], y1[idx]),((x[idx+1], y1[idx+1])), ((x[idx], y2[idx])), ((x[idx+1], y2[idx+1])))
        xcs.append(xc)
        ycs.append(yc)
    return np.array(xcs), np.array(ycs)



# Load the 3D data file
data = np.genfromtxt("Ta_parameterspace_2mm.txt", skip_header=14, delimiter="\t", dtype = float)
#print(data)
reflect = data[:,0]
emiss = data[:,1]
tempdiff = data[:,4]

# Find the node list for x and y 
x_node = np.unique(reflect)
y_node = np.unique(emiss)

# Initialization of a array of 0 
z_zero = np.zeros(len(y_node))

# Reshape the tempdiff 
z_values = tempdiff.reshape(40, 40)

# Find the intersection between the z values along the axis y and 0 
x_intersection = []
y_intersection = []

for idx in range(len(z_zero)):
    z_axis_y = z_values[idx, :]
    y_inter_temp, _ = interpolated_intercepts(y_node, z_zero, z_axis_y)
    # Only if there is an intersection 
    if len(y_inter_temp) >0:
        x_intersection.append(x_node[idx])
        y_intersection.append(y_inter_temp[0][0])
    

fig, ax = plt.subplots()
cb = ax.tricontourf(reflect, emiss, tempdiff, 200, cmap = "seismic")
cbar = plt.colorbar(cb)
cbar.set_label('Temperature (K)', rotation = 270, labelpad = 13)
ax.set_xlabel('Reflectivity')
ax.set_ylabel('Emissivity')
plt.plot(x_intersection, y_intersection, 'k-*')
plt.show()

Then you get the figure :

Plot constant surface on tricontourf in Python

1 Answers1