How to avoid using zeros in interpolation for contourf, for plotting CFD velocity data

Question

From a CFD simulation in Ansys Fluent, I have a dataframe from a csv containing 4 columns: x,y,z (nodal coordinates of irregular grid), and u (streamwise velocity magnitude). I want to plot a contour plot of u over the x-z plane. The solid body around which air flows has u=0, and I want this to be masked for the interpolation, i.e. I want to create a contour plot with an irregular boundary. The first attached image is more or less what I'd like to get with Python, obtained with Ansys CFD-Post -- the geometry of the solid body is outlined in black:

This is what the data looks like. Note column values are not unique or sorted.

### geometric parameters, in m
L = 0.3048 
hwall = 0.0095
Lin = L-hwall*2
H = 0.073 
Lpad = Lin/10
Hpad = H/10

### get df: x,y,z coordinates and instantaneous streamwise velocity
df = pd.read_csv('test2.csv',names=['x','y','z','u'])
dfsub = df.loc[(df.x>-Lin/2-Lpad)&(df.x<Lin/2+Lpad)&(df.z<H+Hpad)] # limit to smaller region

dfsub
x   y   z   u
2   -0.141747   2.106994e-11    0.010645    0.106242
3   -0.139540   -6.458060e-12   0.010615    0.456511
5   -0.132303   -1.308423e-12   0.010602    1.138072
8   0.141747    -2.106730e-11   0.010650    -0.154759
9   0.139543    1.165080e-11    0.010619    -0.150316
... ... ... ... ...
3597    0.159718    -2.567698e-11   0.027614    -0.532316
3598    0.159934    2.542068e-13    0.027258    -0.544210
3599    0.159501    1.936669e-12    0.027613    -0.520241
3600    0.157081    -6.854833e-12   0.035597    -0.392042
3601    0.156974    1.765435e-11    0.027732    -0.382951

Here is a scatterplot to show the data:

umax = max(np.abs(dfsub.u))

fig,ax=plt.subplots(figsize=(16,3))
dfsub.plot.scatter('x','z',c='u',cmap ='seismic',vmin=-umax,vmax=umax,ax=ax)
plt.show()

And here is the basic code for plotting, which I don't expect to work since the solid region is just read as velocity = 0 for interpolation:

### contour plot
crange = 4
fig,ax=plt.subplots(figsize=(16,3))
tcf = ax.tricontourf(dfsub.x,dfsub.z,dfsub.u,levels=np.arange(
crange,crange,0.1),cmap='seismic',extend='both')
plt.colorbar(tcf,ax=ax,label='$u$ (m/s)')
ax.set_title('basic',size=14)
plt.show()

and here is the output:

I've tried to mask these values in several ways, e.g. replace 0s with nans, which gives "ValueError: z array must not contain non-finite values within the triangulation":

### get df: x,y,z coordinates and instantaneous streamwise velocity
df = pd.read_csv('test2.csv',names=['x','y','z','u'])
df.replace(0, np.nan, inplace=True)
dfsub = df.loc[(df.x>-Lin/2-Lpad)&(df.x<Lin/2+Lpad)&(df.z<H+Hpad)] # limit to smaller region

or to remove the 0s from the dataframe:

### get df: x,y,z coordinates and instantaneous streamwise velocity
df = pd.read_csv('test2.csv',names=['x','y','z','u'])
df = df.loc[df.u != 0]
dfsub = df.loc[(df.x>-Lin/2-Lpad)&(df.x<Lin/2+Lpad)&(df.z<H+Hpad)] # limit to smaller region

The output for this is in the attached image. It still interpolates in the regions I mean to exclude.

Answer 1

If I understand your question correctly, you want to use a masked array. When plotting the contours of masked arrays, the masked regions are ignored. Without the actual data to understand how it is formatted, I leave it to you to figure out how to create your masking array.

import numpy as np
import matplotlib.pyplot as plt

plt.close("all")

x = np.linspace(-1, 1, 100)
y = np.linspace(0, 1, 100)
X, Y = np.meshgrid(x, y)

U = X + 2*Y

# create a region where we do not want the contour to exist
solid = (X >= -0.25) & (X <= 0) & (Y >= 0) & (Y <= 0.75)

# the mask argument takes a boolean array that says which elements are masked
U_masked = np.ma.array(U, mask=solid)

fig, ax = plt.subplots()
# contourf will not plot the masked regions
p = ax.contourf(X, Y, U_masked, corner_mask=False)
fig.colorbar(p, ax=ax, label="u")
ax.set_xlabel("x")
ax.set_ylabel("y")

---

Edit: Since the data is irregularly spaced, using tricontourf might work instead. The only issue with this method is that corners will be cut off (or at least I am not sure how to deal with that).

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri

plt.close("all")

x = np.linspace(-1, 1, 100)
y = np.linspace(0, 1, 100)
X, Y = np.meshgrid(x, y)

U = X + 2*Y

# create a region where we do not want the contour to exist
solid = (X >= -0.25) & (X <= 0) & (Y >= 0) & (Y <= 0.75)

# create the triangles; I used the flattened X and Y from before,
# but I think using the x and y values from the data will work
triang = tri.Triangulation(X.flatten(), Y.flatten())

# masking triangles based on the `solid` boolean matrix
triang.set_mask(np.any(solid.flatten()[triang.triangles], axis=1))

fig, ax = plt.subplots()
p = ax.tricontourf(triang, U.flatten())
fig.colorbar(p, ax=ax, label="u")
ax.set_xlabel("x")
ax.set_ylabel("y")

Answer 2

The solution from @jared somewhat works, but I also followed this solution to create a mask for this data.

# 3 geometric conditions for mask 
cond1 = (df.x > -L/2) & (df.x < L/2) & (df.z < hbase)
cond2 = (df.x > -L/2) & (df.x < -Lin/2) & (df.z < H)
cond3 = (df.x > Lin/2) & (df.x < L/2) & (df.z < H)

df.loc[cond1 | cond2 | cond3, 'u'] = 0 # set column u = 0 in those regions

# Extract the values from the DataFrame columns
x = df.x.values
z = df.z.values
u = df.u.values

# Set up the mask based on zero values in the u array
isbad = u == 0

# triangulate, and evaluate indices of triangles
triang = tri.Triangulation(x, z)
mask = np.all(np.where(isbad[triang.triangles], True, False), axis=1)
triang.set_mask(mask)

# Plot the contour using tricontourf with a mask
crange = 4
fig, ax = plt.subplots(figsize=(16, 3))
tcf = ax.tricontourf(triang, u, levels=np.arange(-crange, crange, 0.005), cmap='seismic', extend='both')
plt.colorbar(tcf, ax=ax, label='$u$ (m/s)')
ax.set_ylim(0,0.08)
ax.set_xlim(-0.165,0.165)
plt.xlabel('x')
plt.ylabel('z')
plt.show()

It does seem to prevent interpolation over the zeros, but the triangulation is very coarse, so the mask isn't great.

I guess the solution is to refine the grid + interpolate somehow?