计算科学 - 插值对于网格是正确的，但对于随机点是错误的 - 吾爱随笔录

我正在使用 scipy 插值构建两个坐标到其他两个坐标的映射 --- 实际上一组坐标是 2D 笛卡尔坐标（称为 X 和 Y），另一组类似于极坐标（径向变量 r 和角度坐标t) 但没有关于笛卡尔坐标的解析表达式。变量 r 标记不同的轮廓，t 是轮廓周围的角度。

第一的：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d, RectBivariateSpline
from skimage import measure

我从一个数据集开始，其中包含 r 作为 X、Y 在特定点的函数。X_data 有 nX 点，Y_data 有 nY 点，而 r_data 是 nX x nY 数组。

X_data = np.linspace(1.0, 3.0, 1000)
Y_data = np.linspace(-3.0, 3.0, 1000)
XX, YY = np.meshgrid(X_data, Y_data)
X_centre = 2.0
Y_centre = 0.0
r_data = ((XX-X_centre)**2 + (YY-Y_centre)**2)**0.5 + 0.25*np.sin(YY)

请注意，这只是原始数据的模型，没有分析形式。我在 (r, t) 的规则网格上构建了一个 X 和 Y 的二维数组：

nR, nT = 250, 250
t_array = np.linspace(0.0, 2.0*np.pi, nT)
r_array = np.linspace(0.01, 1.0, nR)
X_grid = np.zeros([ nR, nT])
Y_grid = np.zeros([ nR, nT])

现在我想在 r_data 中找到 r 的轮廓，使用以下函数，该函数使用 skimage 中的 measure.find_contours 来查找轮廓，然后处理该输出以获得沿 r=constant 的完整轮廓的 X、Y 坐标。然后我计算沿轮廓的角度（Contour_t）。

def GetContour(r, r_data, X_data, Y_data, X_centre, Y_centre):    
    Contours = measure.find_contours(r_data, r) # gives X,Y indices that form contour

    Contour_X, Contour_Y = np.array([]), np.array([])
    for contour in Contours: #combine lists of partial contours for full contour
        Contour_X = np.concatenate((Contour_X, X_data[0] + contour[:,0]/len(X_data)*(X_data[-1] - X_data[0]) )) 
        Contour_Y = np.concatenate((Contour_Y, Y_data[0] + contour[:,1]/len(Y_data)*(Y_data[-1] - Y_data[0]) )) 

    # re-arrange contour so angle is taken from Y_data = 0 on right hand side of X_centre (r=0)
    PossibleInds = np.argsort(np.abs(Contour_Y)) #indices for the best candidates for the start
    j = 0
    while True:
        if Contour_X[PossibleInds[j]] >= X_centre:
            if Contour_Y[PossibleInds[j]] >= Y_centre:
                break
        j += 1
    Ind = PossibleInds[j]

    if Ind!=0: # if it's not already in the right order (by chance)
        Contour_X = np.concatenate(( Contour_X[Ind-1:], Contour_X[:Ind]))
        Contour_Y = np.concatenate(( Contour_Y[Ind-1:], Contour_Y[:Ind]))

    Contour_l = np.zeros(len(Contour_X)) # length along the contour
    Contour_l[1:] = np.sqrt(np.diff(Contour_X)**2 + np.diff(Contour_Y)**2 )
    Contour_l = np.cumsum(Contour_l)
    Contour_t = Contour_l/Contour_l[-1] * 2.0*np.pi  # angle along contour

    return Contour_X, Contour_Y, Contour_t

现在找到 r_array 中每个元素的每个轮廓，然后将它们绘制在数据上以进行完整性检查：

for i, r in enumerate(r_array):
   X, Y, Contour_t = GetContour(r, r_data, X_data, Y_data, X_centre, Y_centre)
   fX = interp1d( Contour_t, X) # interpolate X as a function of angle along contour
   fY = interp1d( Contour_t, Y)
   X_grid[i,:] = fX(t_array)    # now save X at fixed values of t given by t_array
   Y_grid[i,:] = fY(t_array)
   plt.plot(X_grid[i,:], Y_grid[i,:])

plt.contourf(X_data, Y_data, r_data.T)
plt.show()

伟大的。我已成功找到轮廓。现在我使用 RectBivariateSpline 来构造 X(r,t) 和 Y(r,t) 的映射：

X_map = RectBivariateSpline( r_array, t_array, X_grid )
Y_map = RectBivariateSpline( r_array, t_array, Y_grid )

我可以通过使用以下方法绘制地图来检查这项工作：

r_test = np.linspace(0.0, 1.0, 50)
t_test = np.linspace(0.0, 2.0*np.pi, 50)
rr, tt = np.meshgrid(r_test, t_test, indexing='ij')
plt.contourf(rr, tt, X_map(r_test, t_test))
plt.show()
plt.contourf(rr, tt, Y_map(r_test, t_test))
plt.show()

它给出了我所期望的。但是，我实际上想在 r、t 中随机生成点，然后使用地图获取随机生成点的 X、Y 坐标。

N = 100000
r_pts = np.random.uniform(0.01, 1.0, size = N)
t_pts = np.random.uniform(0.0, 2.0*np.pi, size = N)
X_pts = X_map(r_pts, t_pts, grid=False)
Y_pts = Y_map(r_pts, t_pts, grid=False)

plt.hist2d(X_pts, Y_pts, 100)
plt.show()

现在所有的点都聚集在中心，没有达到 r=1.0，而且它们的角度肯定是不均匀的！

问题（S）：我怀疑我搞砸了索引，但我玩过并且看不到明显的错误。你可以吗？为什么绘制地图时看起来不错，但在随机生成的点上使用地图会得到错误的结果？

我尝试了各种插值方法（interp2d、RegularGridInterpolator、map_coordinates），这些都有同样的问题......

谢谢！