Python griddata与MATLAB scatteredInterpolant结果差异及替代方案咨询

问题背景 #

我正在将一段涉及三维插值的MATLAB代码迁移至Python。原本认为MATLAB的scatteredInterpolant与Python中scipy.interpolate.griddata底层机制一致（均基于网格数据的Delaunay三角剖分及线性插值），但用相同数据测试时，部分插值结果存在差异（部分结果一致）。推测这是由于Delaunay三角剖分不唯一，存在多种有效解导致。但需求是获得完全一致的结果，请问如何解决该问题，或是否有更合适的Python库？已尝试其他Python方法，但均不符合需求。

MATLAB测试代码 #

%MATLAB code on some made up data (real data behaves similarly):

depth = [10;2;0;2;1;0;0;10;1;1.1;7.5;7;9;9.9;8;5];
lat = [5;5;0;0;4;2.7;0;0;10;6;3;4.4;3.3;2;9;9];
lon = [0;0;8;4;3;2.5;0;10;5;4.4;7;9.9;1;6.5;8;7];

C_grid = [99;87;88;89;78;79;65;67;75;75;75;75;75;75;75;75];

lat_grid = [0;0;0;0;10;10;10;10;-50;-50;-50;-50;50;50;50;50];
lon_grid = [0;10;0;10;0;10;0;10;-50;50;-50;50;-50;50;-50;50];
depth_grid = [0;0,;10;10;0;0;10;10;-50;-50;50;50;-50;-50;50;50];

interp = scatteredInterpolant(lon_grid, lat_grid, depth_grid, C_grid);
interpolation = interp(lon, lat, depth);

%RESULT: 76.5, 86.3, 89.4, 92.0, 85.9, 90.3, 99.0, 89.0, 77.2, 80.2, 81.3, 82.9, 80.3, 84.0, 71.2, 74.6

Python测试代码（相同数据） #

import numpy as np

depth = np.array((10, 2, 0, 2, 1, 0, 0, 10, 1, 1.1, 7.5, 7, 9, 9.9, 8, 5))
lat = np.array((5,5, 0, 0, 4, 2.7, 0, 0, 10, 6, 3, 4.4, 3.3, 2, 9, 9))
lon = np.array((0, 0, 8, 4, 3, 2.5, 0, 10, 5, 4.4, 7, 9.9, 1, 6.5, 8, 7))

lat_grid = np.array((0,0,0,0,10,10,10,10,-50,-50,-50,-50,50,50,50,50))
lon_grid = np.array((0,10,0,10,0,10,0,10,50,50,-50,50,-50,50,-50,50))
depth_grid = np.array((0,0,10,10,0,0,10,10,-50,-50,50,50,-50,-50,50,50))
C_grid = np.array((99,87,88,89,78,79,65,67,75,75,75,75,75,75,75,75))

points = (lon_grid, lat_grid, depth_grid)
values = C_grid
xi = (lon, lat, depth)

import scipy
from scipy.interpolate import griddata

interpolation = griddata(points, values, xi, method="linear")
print(interpolation)

#RESULT: 76.5, 86.3, 89.4, 92.0, 87.1, 90.3, 99.0, 89.0, 77.3, 81.23, 81.6, 78.7, 81.9, 84.2, 71.3, 73.5

具体问题 #

1. 两种方法为何部分结果不同？
2. 是否存在可完全复现MATLAB scatteredInterpolant的Python方法？

问题解答 #

1. 结果差异的核心原因 #

你推测的Delaunay三角剖分不唯一性是关键，具体细节包括：

• 剖分实现差异：MATLAB和SciPy依赖的Delaunay剖分算法（均基于Qhull）在实现细节上有区别，比如共面点的处理逻辑、点排序优先级、四面体划分规则，这些差异会导致三维空间中生成的三角网格结构不同。
• 插值权重计算差异：线性插值的权重依赖于四面体的顶点位置，不同的剖分结构会计算出不同的权重，最终导致插值结果偏差。
• 数值精度与边界处理：两者在近邻点判定阈值、浮点精度控制上的细微差别，也会对部分边界点或近边界点的插值结果产生影响。

2. 完全复现MATLAB结果的方案 #

方案一：调用MATLAB引擎执行原生逻辑

这是最可靠的方式，直接复用MATLAB的scatteredInterpolant实现，确保结果完全一致：

import matlab.engine
import numpy as np

# 启动MATLAB引擎
eng = matlab.engine.start_matlab()

# 将NumPy数组转换为MATLAB支持的数组格式
depth_mat = matlab.double(depth.tolist())
lat_mat = matlab.double(lat.tolist())
lon_mat = matlab.double(lon.tolist())
C_grid_mat = matlab.double(C_grid.tolist())
lat_grid_mat = matlab.double(lat_grid.tolist())
lon_grid_mat = matlab.double(lon_grid.tolist())
depth_grid_mat = matlab.double(depth_grid.tolist())

# 将变量传入MATLAB工作区
eng.workspace['lon_grid'] = lon_grid_mat
eng.workspace['lat_grid'] = lat_grid_mat
eng.workspace['depth_grid'] = depth_grid_mat
eng.workspace['C_grid'] = C_grid_mat

# 执行MATLAB插值逻辑
eng.eval("interp = scatteredInterpolant(lon_grid, lat_grid, depth_grid, C_grid);", nargout=0)

# 传入待插值点并计算结果
eng.workspace['lon'] = lon_mat
eng.workspace['lat'] = lat_mat
eng.workspace['depth'] = depth_mat
interpolation_mat = eng.eval("interp(lon, lat, depth);")

# 转换回NumPy数组
interpolation = np.array(interpolation_mat)
print(interpolation)

# 关闭MATLAB引擎
eng.quit()

该方案需要本地安装MATLAB及对应版本的MATLAB Python引擎，适合对结果一致性要求极高的场景。

方案二：手动对齐剖分与插值逻辑（难度较高）

若无法依赖MATLAB，需要深入分析MATLAB的scatteredInterpolant底层逻辑：

1. 研究MATLAB的delaunayTriangulation三维剖分规则，特别是共面点、边界点的处理方式。
2. 在SciPy中使用scipy.spatial.Delaunay时，通过调整qhull_options参数（如'QJ'合并共面点、'QbB'启用边界约束）尝试匹配MATLAB的剖分结果。
3. 手动实现线性插值的权重计算，确保与MATLAB的插值公式完全一致。
   但这种方案需要对两种工具的底层实现有深入理解，且难以覆盖所有场景的完全对齐。

内容的提问来源于stack exchange，提问作者Larissa Dias