我现在遇到这么个问题：用Python处理CFD导出的二维数据，数据是在网格质心记录的，但我想把它插值到规则间距的点网格上展示。本来想着用scipy.interpolate.griddata来做，但卡在了怎么给非矩形平面生成规则点网格这一步——总不能手动拆成小矩形吧？有没有现成的函数能搞定这事？

Update1：我跟着评论里Christian Karcher的建议试了试：先生成均匀间距的矩形网格，再把CFD数据插值上去，但结果有点奇怪——线性插值没能排除那些没有数据的区域点。

Update2：后来看了Christian在回答下的补充评论，发现是「远距离插值」搞的鬼！解决办法是用scipy.spatial.KDTree过滤掉矩形网格里离原始数据点太远的点，只保留邻近区域的点。我把修正后的代码贴出来了，还附了实际数据、网格和插值结果的对比图。

Update3：又遇到新问题了——用matplotlib的contourf画等高线的时候，发现原来用meshgrid生成的xx、yy是结构化的（形状是[[x1,x2...xn],[x1,x2...xn]...]），但过滤后的xx_filtered变成了一维数组，网格结构没了。有没有办法过滤数据的同时保留原始网格的结构啊？

完整代码（已修正） #

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_magnitude'])


#---------------------- 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
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

# scipy.spatial.KDTree function is used to lookup nearest-neighbour for defined points

#define input coordinates (imporeted data point coordinates) as numpy arrays

print(df_results['X'].to_numpy())
print(type(df_results['X'].to_numpy()))

#input_coords = np.concatenate([df_results['X'].to_numpy()[:, None], df_results['Y'].to_numpy()[:, None]], axis =1)
input_coords = np.array(list(zip(df_results['X'].to_numpy(), df_results['Y'].to_numpy())))
tree = scipy.spatial.KDTree(input_coords)
dist, idx = tree.query(xy_mesh)
radius = (rect_mesh_spacing + rect_mesh_spacing) / 2
distance_filter = dist <= radius

xx_filtered, yy_filtered = xy_mesh[distance_filter].T



#-------------------- INTERPOLATE DATA TO RECTANGULAR GRID ----------

# 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
interpolated_data_nearest = scipy.interpolate.griddata((df_results['X'], df_results['Y']), df_results['p_stat'],
                                               (xx_filtered, yy_filtered), method = 'nearest')
interpolated_data_linear = scipy.interpolate.griddata((df_results['X'], df_results['Y']), df_results['p_stat'],
                                               (xx_filtered, yy_filtered), method = 'linear')
interpolated_data_cubic = scipy.interpolate.griddata((df_results['X'], df_results['Y']), df_results['p_stat'],
                                               (xx_filtered, yy_filtered), method = 'cubic')

#------------------------- PLOT ---------------------------------
fig, axs = plt.subplots(1,5)
fig.set_size_inches(13, 4)

# Plot non-interpolated data
z = df_results['p_stat']
axs[0].scatter(df_results['X'], df_results['Y'], c = z,
              norm = matplotlib.colors.Normalize(vmin = np.min(z), vmax = np.max(z)), cmap = 'bwr', s=2)


# Plot rectangular uniform empty grid
axs[1].plot(xx, yy, marker='.', color = 'k', linestyle='none')

# Plot interpolated data
z1 = interpolated_data_nearest
z2 = interpolated_data_linear
z3 = interpolated_data_cubic

axs[2].scatter(xx_filtered, yy_filtered, c =  z1, norm = matplotlib.colors.Normalize(vmin = np.min(z1), vmax = np.max(z1)),
            cmap = 'bwr', s=2)
axs[3].scatter(xx_filtered, yy_filtered, c =  z2, norm = matplotlib.colors.Normalize(vmin = np.min(z2), vmax = np.max(z2)),
            cmap = 'bwr', s=2)
axs[4].scatter(xx_filtered, yy_filtered, c =  z3, norm = matplotlib.colors.Normalize(vmin = np.min(z3), vmax = np.max(z3)),
            cmap = 'bwr', s=2)



plt.show()

插值结果对比图 #

备注：内容来源于stack exchange，提问作者user2882635