基于Python与scipy/numpy的两步1D插值正确实现方案咨询

问题原因分析 #

原代码存在两个核心问题：

1. 行方向插值未关闭边界错误，只要某行x范围不覆盖输入值就直接抛出异常，哪怕其他行有有效覆盖
2. Y方向插值时没有过滤X方向插值返回的nan值，导致只要有一行返回nan就会污染最终结果

极简实现方案 #

仅需少量修改原逻辑，完全基于numpy和scipy现有能力实现，无需手写复杂插值逻辑：

import numpy as np
from scipy.interpolate import interp1d
from typing import Callable, ArrayLike

def twice_interp1d_with_2d_x(
    x: ArrayLike, 
    y: ArrayLike, 
    z: ArrayLike,
    bounds_error: bool = True,
    fill_value = None
) -> Callable[[ArrayLike, ArrayLike], np.ndarray]:
    # 参数预处理与校验
    x = np.asarray(x)
    y = np.asarray(y)
    z = np.asarray(z)
    if x.ndim != 2 or y.ndim != 1 or z.ndim != 1:
        raise ValueError("x必须为2D数组，y、z必须为1D数组")
    if x.shape[0] != len(y):
        raise ValueError("x的行数需与y的长度一致")
    if x.shape[1] != len(z):
        raise ValueError("x的列数需与z的长度一致")
    
    # 预先生成每行的X方向插值器，超界自动返回nan不报错
    row_interps = [
        interp1d(row, z, bounds_error=False, fill_value=np.nan)
        for row in x
    ]

    @np.vectorize
    def _interp_single_point(xi: float, yi: float) -> float:
        # 计算所有行对当前x的插值结果
        temp_zs = np.array([f(xi) for f in row_interps])
        # 过滤有效结果
        valid_mask = ~np.isnan(temp_zs)
        valid_y = y[valid_mask]
        valid_z = temp_zs[valid_mask]

        if len(valid_y) == 0:
            return np.nan
        elif len(valid_y) == 1:
            # 单点匹配直接返回
            return valid_z[0] if np.isclose(yi, valid_y[0]) else np.nan
        else:
            # 多有效点做Y方向插值
            sort_idx = np.argsort(valid_y)
            y_interp = interp1d(
                valid_y[sort_idx], valid_z[sort_idx],
                bounds_error=bounds_error, fill_value=fill_value
            )
            return float(y_interp(yi))

    def interpolator(x_inner: ArrayLike, y_inner: ArrayLike) -> np.ndarray:
        # 自动支持数组/标量输入
        return _interp_single_point(x_inner, y_inner)

    return interpolator

测试验证 #

用示例数据测试：

# 示例数据
x = [[1,2,3], [3,4,5]]
y = [1,2]
z = [1,2,3]

f = twice_interp1d_with_2d_x(x, y, z)
# 测试提到的异常case
print(f(4, 2)) 
# 输出: 2.0，符合预期

特性说明 #

• 自动支持标量、数组输入，输出形状与输入广播后的形状一致
• 仅在无任何有效插值数据的位置返回nan，只要存在有效行就会尝试计算
• 保留原有的bounds_error和fill_value参数可配置，Y方向插值行为和原生interp1d完全一致
• 代码逻辑简洁，所有插值能力复用scipy成熟实现，无需手写底层插值逻辑

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