17370845950

`regulargridinterpolator` 不支持含 nan 的输入数据，需改用 `griddata` 等无结构插值器：先剔除 nan 网格点及其对应坐标，再以 `(x_flat, y_flat, z_flat)` 形式调用 `griddata(..., method='cubic')`，即可安全完成三次插值并天然处理边界与空洞区域。

当面对带大量 NaN 的二维规则地理网格（如 1°×1° 的纬度-经度场）并对任意散点坐标（如浮标位置或模型输出点）执行高阶平滑插值时，scipy.interpolate.RegularGridInterpolator 的 'cubic' 模式会直接报错 ValueError: Array must not contain infs or nans.——这是因为其底层依赖 make_interp_spline，而该函数严格要求输入值数组完全有限（finite），不接受任何缺失值。

根本原因在于：RegularGridInterpolator 是为“完整规则网格”设计的；它假设每个 (lat_i, lon_j) 都有定义良好的值，从而构建张量积样条基函数。一旦出现 NaN，不仅破坏了网格完整性，更导致沿任一维度做一维样条拟合时失败（如错误栈中 _do_spline_fit 所示）。

✅ 正确解法是切换到无结构（unstructured）插值范式：将原始网格展平为有效点云，再使用 scipy.interpolate.griddata。该函数原生支持 method='cubic'（基于分片双三次 Hermite 插值），且自动忽略无效点、鲁棒处理空洞邻域，并对查询点外推时返回 NaN（符合预期行为）。

以下是推荐实现流程（适配您的地理数据场景）：

import numpy as np
from scipy.interpolate import griddata

# 假设 lat0 (shape: M), lon0 (shape: N), data_nan (shape: M×N)
# 其中 data_nan 含约 40% NaN，但所有查询点 (lat, lon) 均在网格范围内

# 1. 构建完整坐标网格并展平
latG, lonG = np.meshgrid(lat0, lon0, indexing='ij')  # shape: (M, N)
lat_flat = latG.ravel()
lon_flat = lonG.ravel()
data_flat = data_nan.ravel()

# 2. 筛选非 NaN 有效点（保留坐标与值的一致性）
mask_valid = ~np.isnan(data_flat)
points_valid = np.column_stack((lat_flat[mask_valid], lon_flat[mask_valid]))
values_valid = data_flat[mask_valid]

# 3. 对目标散点 (lat, lon) 执行三次插值
# 注意：griddata 要求查询点为 (N, 2) 形状的数组
xi = np.column_stack((lat, lon))  # shape: (K, 2)

interped_cub = griddata(
    points=points_valid,
    values=values_valid,
    xi=xi,
    method='cubic',
    fill_value=np.nan  # 显式指定越界/空洞处返回 NaN（默认即此行为）
)

# interped_cub.shape == (len(lat),) —— 与 linear 插值结果维度一致！

? 关键优势说明：

• ✅ 自动容错：griddata 在 'cubic' 模式下会基于 Delaunay 三角剖分，仅利用局部有效邻点构造插值多项式，天然跳过 NaN 区域；
• ✅ 语义一致：对靠近 NaN 区域的查询点，若缺乏足够邻近有效支撑点，结果自动为 NaN，与 RegularGridInterpolator(method='linear') 的失败行为逻辑统一；
• ✅ 无需后处理：避免了“先线性 + 再最近邻填补”的两阶段繁琐流程；
• ⚠️ 注意性能：griddata 的 'cubic' 比 'linear' 计算开销大，但对数千个查询点（如您的 ~2000）仍非常高效；若需更高性能，可考虑 scipy.interpolate.CloughTocher2DInterpolator（同样支持 fill_value，且缓存三角剖分）。

? 补充建议：

• 若地理跨度大（如跨极区），建议先将经纬度转换为等距投影坐标（如 pyproj 中的 EPSG:3857 或 EPSG:4326 下的球面距离加权），以提升插值几何合理性；
• 对于超大规模网格（>10⁵ 有效点），可启用 rescale=True 参数缓解坐标尺度差异导致的数值不稳定；
• 验证时，可用小范围人工 NaN 区域对比 griddata 与手动掩膜后 RegularGridInterpolator('linear') 结果，确认空洞边缘过渡平滑性。

至此，您即可在保持代码简洁性的同时，获得比线性插值更光滑、比最近邻更物理合理的三次插值结果，且完全兼容真实地球科学数据中常见的稀疏、不规则缺失模式。