含 NaN 二维网格三次样条插值方法

`RegularGridInterpolator` 不支持含 NaN 的输入数据，需改用 `griddata` 对有效点进行非结构化三次插值，并自动处理边界与空缺区域。

在科学计算与地球系统建模中，常需将规则经纬度网格（如 1°×1°）上的物理量（如温度、降水）插值到任意散点坐标（如观测站点或轨迹点）。当原始数据含大量 NaN（例如海洋区域无观测、云遮挡导致遥感缺失），scipy.interpolate.RegularGridInterpolator 的 'cubic' 方法会直接报错 ValueError: Array must not contain infs or nans.——这是因为其底层调用 make_interp_spline 要求输入值数组完全非空且有限，无法容忍任何缺失值。

而 'linear' 方法能成功运行，是因为线性插值仅依赖局部邻近点的加权平均，RegularGridInterpolator 内部会对 NaN 区域做隐式跳过与边界截断；但三次样条插值需构建全局光滑样条函数，必须基于完整、连续的规则网格，故不兼容稀疏/含缺数据。

✅ 正确解法：转向非结构化插值器 scipy.interpolate.griddata
griddata 专为“散点 → 散点”或“散点 → 规则网格”设计，天然支持从任意子集坐标-值对（即剔除 NaN 后的有效点）出发，执行 'nearest'、'linear' 或 'cubic' 插值。它不要求输入点构成规则网格，也不要求覆盖整个定义域，因此完美适配您的场景。

✅ 推荐实现步骤（含完整示例）

import numpy as np
from scipy.interpolate import griddata

# 1. 构造原始规则网格与含 NaN 数据（同问题复现逻辑）
lat0 = np.linspace(-60, -40, 20)
lon0 = np.linspace(-50, 0.0, 50)
lonG, latG = np.meshgrid(lon0, lat0)
data = np.sin(latG * lonG)
data_nan = data.copy()
data_nan[1:8, 5:12] = np.nan
data_nan[2:6, 12:20] = np.nan

# 2. 提取所有非 NaN 点的坐标与值（展平 + 过滤）
mask = ~np.isnan(data_nan)
lat_valid = latG[mask]   # 1D array of valid latitudes
lon_valid = lonG[mask]   # 1D array of valid longitudes
data_valid = data_nan[mask]  # 1D array of corresponding values

# 3. 准备待插值的“非结构化”目标坐标（例如 2000 个随机点）
lat_q = np.random.uniform(-60, -40, 2000)
lon_q = np.random.uniform(-50, 0.0, 2000)
xi = np.column_stack((lat_q, lon_q))  # shape (2000, 2)

# 4. 执行三次样条插值（关键！）
interped_cubic = griddata(
    points=(lat_valid, lon_valid),  # 注意：顺序需与 xi 一致（此处为 (lat, lon)）
    values=data_valid,
    xi=xi,
    method='cubic',
    fill_value=np.nan  # 显式指定域外/不可插区域返回 NaN
)

# 5. 验证结果（自动包含 NaN 处理）
print(f"插值完成，共 {np.isnan(interped_cubic).sum()} 个 NaN 输出")

⚠️ 关键注意事项

• 坐标顺序一致性：griddata(points=(x, y), ...) 中 x, y 的顺序必须与 xi 的列顺序严格一致。若 xi 是 (lat, lon)，则 points 应为 (lat_valid, lon_valid)；若习惯 (lon, lat)，请同步调整。
• fill_value 行为：默认 fill_value=nan，对超出凸包范围（convex hull）的查询点返回 NaN。您可设为 0 或其他占位值，但建议保留 np.nan 以明确标识插值失败区域。
• 性能提示：griddata(method='cubic') 在高维或大数据量时较慢（时间复杂度约 O(N³)）。若目标点极多（>10⁵），可考虑先插值到中间规则网格，再用 RegularGridInterpolator 二次加速。
• 替代方案：对于地理坐标（经纬度），需注意球面距离畸变。若精度要求极高，应先将 (lat, lon) 投影到平面坐标系（如 UTM），再插值，最后反投影。

✅ 总结

RegularGridInterpolator(..., method='cubic') 本质是规则网格专用的张量积样条，拒绝 NaN；而 griddata(..., method='cubic') 是基于 Delaunay 三角剖分的散点样条，天生适配稀疏数据。只需三步：① 展平并过滤 NaN 点；② 组织 (x, y, z) 三元组；③ 调用 griddata —— 即可获得与线性插值一致的容错行为（有效区三次光滑，无效区自动置 NaN），且代码简洁、语义清晰。