I managed to finally put together a piece of code that does what I want. The problem was that the filtering of the data removes the 2D structure of the coordinates and returns a vector form of the x/y coordinates of the remaining points.
I then transformed the filtered data back in 2D form and replaced the filtered out coordinates with NaN values. This was also done for the intepolated pressure data. This was then used for plotting with .contourf without issues.
Bellow is the complete code, if someone wants to test it, they can use the input data posted bellow. I left all the comments from my code, if moderators think thats too much, I will remove them.
import pandas as pd
import numpy as np
import scipy
import matplotlib.pyplot as plt
import matplotlib
#-------------------------- INPUT ------------------------------
# Read data
input_filename = 'period_rot_shadow_vel-mag_p-stat.txt'
df_results = pd.read_csv(input_filename,
sep =',', skiprows = 1,
names = ['nodenumber', 'X', 'Y', 'Z','p_stat', 'vel_mag'])
#---------------------- CREATE RECTANGULAR GRID ----------------
# Rectangular grid is created with a numpy.meshgrid function, which takes 2 1D arrays to define the
#size of the rectangular grid and the spacing between the points
#get the size of the data to generate the same size array
x_data_min = df_results['X'].min()
x_data_max = df_results['X'].max()
y_data_min = df_results['Y'].min()
y_data_max = df_results['Y'].max()
#from max and min values calculate the range (size) for both axis
x_size = np.abs(x_data_max) + np.abs(x_data_min)
y_size = np.abs(y_data_max) + np.abs(y_data_min)
#create 2 numpy arrays to define the size and the spacing of the rectangular grid, arrays
#are defined as lists, so we can use the list comprehension to use the data size parameters and the
#desired spacing
#all values are in meters
rect_mesh_spacing = 0.001
#number of points on each edge is calculated from the x_size and spacing, +1 is added so that we always
#cover the whole area covered in exported data
N_x_steps = int(x_size/rect_mesh_spacing) + 1
N_y_steps = int(y_size/rect_mesh_spacing) + 1
#create an empty list to fill with the values
rect_mesh_x_axis = np.array([])
rect_mesh_y_axis = np.array([])
#for loop to fill the list of evenly spaced value in the x axis
for i in range(0, N_x_steps):
x = x_data_min + i * rect_mesh_spacing
rect_mesh_x_axis = np.append(rect_mesh_x_axis, x)
for i in range(0, N_y_steps):
y = y_data_min + i * rect_mesh_spacing
rect_mesh_y_axis = np.append(rect_mesh_y_axis, y)
#use meshgrid function to create a rectangular grid, meshgrid returns a tuple of ndarrays
xx, yy = np.meshgrid(rect_mesh_x_axis, rect_mesh_y_axis)
#create an array from meshgrid points
xy_mesh = np.array((xx,yy)).T
#---------------- FILTER RECTANGULAR GRID BASED ON DISTANCE TO ORIGINAL POINTS ----------------
# If we would not filter out only the points that are in some vicinity to the original data points
#exported from Fluent, the interpolation function would interpolate also to the points of rectangular grid
#that are far away from the actual exported data points
# Define input coordinates (imported data point coordinates) as numpy arrays
input_coords = np.array(list(zip(df_results['X'].to_numpy(), df_results['Y'].to_numpy())))
# scipy.spatial.KDTree function is used to lookup nearest-neighbour for defined points
tree = scipy.spatial.KDTree(input_coords)
dist, idx = tree.query(xy_mesh)
#define the radius of lookup of vicinity of points
radius = (rect_mesh_spacing + rect_mesh_spacing) / 2
distance_filter = dist <= radius
# Filter out the points in rectangular grid that are close to the exported data points
#this type of filter produces vectors, 1D arrays for each coordinate, filled with x/y coordinates for each
#point that was not filtered out. The filtered points are not included
xx_filt, yy_filt = xy_mesh[distance_filter].T
#-------------------- INTERPOLATE DATA -------------------------------------------------------
# scipy.interpolate.griddata(points, values, xi)
#points = data point coordinates of the values we want to interpolate
#values = the values at the points
#xi = points at which to interpolate data
#we will interpolate using the NEAREST method
#because the input data is in vector form, p_stat_interpol will also be in vector form
p_stat_interpol= scipy.interpolate.griddata((df_results['X'], df_results['Y']), df_results['p_stat'],
(xx_filt, yy_filt), method = 'nearest')
#----------------------- STRUCTURING DATA TO GRID -----------------------------------------------
# As mentioned, the filtered and interpolated data is in vector form and has no connection to the grid and
#this is the problem for contourf plot function, which needs pressure data in 2D form. One possible way is to
#transform the vector form to the 2D form and replace the missing values (of the points that we filtered out
#before) with NaN
# Numpy unique function is used to get the unique values for coordinates, this way we get the grid axes
x_unique = np.unique(xx_filt)
y_unique = np.unique(yy_filt)
# We construct a new grid of NaN values
p_stat_grid = np.full((len(y_unique), len(x_unique)), np.nan)
# Fill the grid data with the p_stat_interpol values where possible, other values will stay NaN
for xi, yi, zi in zip(xx_filt, yy_filt, p_stat_interpol):
ix = np.where(x_unique == xi)[0][0]
iy = np.where(y_unique == yi)[0][0]
p_stat_grid[iy, ix] = zi
# Now we create a new mesh for plotting, this mesh will be in structured 2D form (we could probably just use the
#first xx, yy mesh from the beginning of the code
X_grid, Y_grid = np.meshgrid(x_unique, y_unique)
#------------------------- PLOT PRESSURE CONTOURS ---------------------------------
fig, axs = plt.subplots(1,4)
fig.set_size_inches(18, 6)
# Define number of contour levels
contour_levels = 50
# Plot non-interpolated data
axs[0].set_title('xx,yy')
axs[0].scatter(xx, yy, color = 'red')
axs[1].set_title('xx_f,yy_f')
axs[1].scatter(xx_filt, yy_filt, color = 'red')
axs[2].set_title('scat xx_f,yy_f,p_int')
axs[2].scatter(xx_filt, yy_filt, c = p_stat_interpol,
norm = matplotlib.colors.Normalize(vmin = np.min(p_stat_interpol),
vmax = np.max(p_stat_interpol)),
cmap = 'bwr')
axs[3].set_title('cont xx,yy, p_int')
axs[3].contourf(X_grid, Y_grid, p_stat_grid, contour_levels, cmap='bwr')
plt.show()
Data Sample:
nodenumber, x-coordinate, y-coordinate, z-coordinate, pressure,velocity-magnitude 1,-2.745435307E-02, 2.104994697E-03,-2.185633883E-02,-8.094133146E+03, 4.686428765E+01 2,-2.738908254E-02, 3.158549336E-03,-2.185629997E-02,-8.104864467E+03, 4.606212337E+01 3,-2.757460451E-02, 5.167262860E-03,-2.185623337E-02,-8.093051020E+03, 4.495193692E+01 4,-2.733931890E-02, 5.656880572E-03,-2.185622490E-02,-8.116433300E+03, 4.405860916E+01 5,-2.759126430E-02, 6.070293390E-03,-2.185618260E-02,-8.084872243E+03, 4.428866265E+01 6,-2.737704207E-02, 6.507508463E-03,-2.185627796E-02,-8.106663765E+03, 4.361527022E+01 7,-2.762655064E-02, 0.000000000E+00,-2.185611682E-02,-8.057108514E+03, 7.720657592E+01 8,-2.762655053E-02, 1.920716437E-05,-2.185609461E-02,-8.057029042E+03, 7.591235710E+01 9,-2.762654491E-02, 4.228439350E-05,-2.185606795E-02,-8.057013286E+03, 7.368296305E+01 10,-2.762653069E-02, 6.997596792E-05,-2.185603598E-02,-8.057068875E+03, 7.113144845E+01 11,-2.762650289E-02, 1.032044066E-04,-2.185599766E-02,-8.057152704E+03, 6.841504839E+01 12,-2.762645410E-02, 1.430766767E-04,-2.185595174E-02,-8.057232638E+03, 6.565721945E+01 13,-2.740227795E-02, 1.197447086E-03,-2.185604704E-02,-8.086002417E+03, 4.769658095E+01 14,-2.729995772E-02, 4.555508762E-03,-2.185613769E-02,-8.113811400E+03, 4.473595941E+01 15,-2.762637329E-02, 1.909744910E-04,-2.185589666E-02,-8.057315245E+03, 6.283751408E+01 16,-2.762624634E-02, 2.484489809E-04,-2.185583068E-02,-8.057450247E+03, 5.993119264E+01 17,-2.762465252E-02, 6.184882281E-04,-2.185591455E-02,-8.059810059E+03, 5.089516358E+01 18,-2.762605091E-02, 3.174151315E-04,-2.185575167E-02,-8.057686244E+03, 5.707278224E+01 19,-2.762531127E-02, 4.994699682E-04,-2.185577109E-02,-8.058658000E+03, 5.248349769E+01 20,-2.762366891E-02, 7.612957389E-04,-2.185583175E-02,-8.062030940E+03, 4.956215663E+01