如何正确配置scipy.ndimage.map_coordinates实现3D插值？

核心问题分析 #

你当前的代码存在两个致命错误：

• 输入数组格式错误：map_coordinates要求第一个参数是N维数值网格数组（这里对应3D的(t,x,y)网格），但你传入的func_values是(1780020000,4)的「坐标+值」列表，完全不符合函数的输入要求。
• 坐标参数含义误解：coordinates参数需要传入的是数组的索引坐标（比如t轴索引范围是0到网格点数-1），而非原始物理坐标（t_new/x_new/y_new），直接传入物理坐标会导致维度不匹配报错。

正确实现步骤 #

步骤1：将离散函数值转换为3D网格数组 #

首先要把func_values的(t,x,y,value)结构转换成shape为(num_t, num_x, num_y)的3D数组，其中num_t/num_x/num_y分别是t、x、y轴上的网格点数：

import numpy as np

# 从func_values中提取并排序所有唯一网格坐标
t_coords = np.unique(func_values[:, 0])
x_coords = np.unique(func_values[:, 1])
y_coords = np.unique(func_values[:, 2])
num_t, num_x, num_y = len(t_coords), len(x_coords), len(y_coords)

# 初始化3D网格数组，默认填充0
grid_values = np.zeros((num_t, num_x, num_y))

# 向量化填充网格（避免循环，处理大数据更快）
t_idx = np.searchsorted(t_coords, func_values[:, 0])
x_idx = np.searchsorted(x_coords, func_values[:, 1])
y_idx = np.searchsorted(y_coords, func_values[:, 2])
grid_values[t_idx, x_idx, y_idx] = func_values[:, 3]

步骤2：将物理坐标转换为数组索引坐标 #

把新数据集中的t_new/x_new/y_new转换成网格的浮点索引（map_coordinates支持非整数索引，用于插值计算）：

def physical_to_index(physical_vals, grid_coords):
    # 用线性插值将物理坐标映射为网格索引
    return np.interp(physical_vals, grid_coords, np.arange(len(grid_coords)))

# 转换新数据集的物理坐标为索引
t_indices = physical_to_index(df['t_new'].values, t_coords)
x_indices = physical_to_index(df['x_new'].values, x_coords)
y_indices = physical_to_index(df['y_new'].values, y_coords)

# 组合成map_coordinates要求的格式：(3, N)，N为新数据点数量
coordinates = np.vstack([t_indices, x_indices, y_indices])

步骤3：调用map_coordinates完成插值 #

现在可以正确执行插值，还可通过order参数指定插值阶数（0=最近邻，1=线性，3=三次样条，默认是3）：

from scipy.ndimage import map_coordinates

def interpolator_3d(df, func_values):
    # 构建3D网格数组
    t_coords = np.unique(func_values[:, 0])
    x_coords = np.unique(func_values[:, 1])
    y_coords = np.unique(func_values[:, 2])
    num_t, num_x, num_y = len(t_coords), len(x_coords), len(y_coords)
    
    grid_values = np.zeros((num_t, num_x, num_y))
    t_idx = np.searchsorted(t_coords, func_values[:, 0])
    x_idx = np.searchsorted(x_coords, func_values[:, 1])
    y_idx = np.searchsorted(y_coords, func_values[:, 2])
    grid_values[t_idx, x_idx, y_idx] = func_values[:, 3]
    
    # 转换物理坐标为索引坐标
    def physical_to_index(physical_vals, grid_coords):
        return np.interp(physical_vals, grid_coords, np.arange(len(grid_coords)))
    
    t_indices = physical_to_index(df['t_new'].values, t_coords)
    x_indices = physical_to_index(df['x_new'].values, x_coords)
    y_indices = physical_to_index(df['y_new'].values, y_coords)
    coordinates = np.vstack([t_indices, x_indices, y_indices])
    
    # 执行插值
    interpolated_values = map_coordinates(grid_values, coordinates, order=3)
    return interpolated_values

额外优化提示 #

如果你的网格是规则均匀分布的，也可以使用scipy.interpolate.RegularGridInterpolator，它不需要手动转换索引，用法更直观。

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