如何对含缺失时间步的三维(lat,lon,time)数组做时间线性插值

问题背景 #

有一组维度为lat、lon、time的xarray格式网格数据，其中某一时间步全为NaN，但该时间步的前后时间步均有有效数据。原尝试使用scipy.interpolate.griddata做全域插值，但该方法会同时利用空间和时间邻点，导致插值结果不符合仅用时间维度线性插值的需求。

测试数据 #

import numpy as np

# 三维数据：[time, lat, lon]
test_data = np.array([
    [
        [0.47942554, 0.48971277, 0.47942554, 0.98971277, 1.23971277],
        [0.2171174 , 0.4671174 , 0.7171174 , 0.9671174 , 1.2171174 ],
        [0.18102272, 0.43102272, 0.68102272, 0.93102272, 1.18102272],
        [0.47942554, 0.38367293, 0.63367293, 0.88367293, 1.13367293]
    ],
    [
        [np.nan, np.nan, np.nan, np.nan, np.nan],
        [np.nan, np.nan, np.nan, np.nan, np.nan],
        [np.nan, np.nan, np.nan, np.nan, np.nan],
        [np.nan, np.nan, np.nan, np.nan, np.nan],
        [np.nan, np.nan, np.nan, np.nan, np.nan]
    ],
    [
        [0.47942554, 0.47942554, 0.97942554, 1.22942554, 1.47942554],
        [0.46452136, 0.71452136, 0.96452136, 1.21452136, 1.46452136],
        [0.42073549, 0.67073549, 0.92073549, 1.17073549, 0.47942554],
        [0.35079033, 0.60079033, 0.85079033, 1.10079033, 1.35079033],
        [0.47942554, 0.50903472, 0.75903472, 0.47942554, 1.25903472]
    ]
])

原方法问题 #

原代码将三维数据拉平为点云后使用griddata插值，会同时考虑空间（lat/lon）和时间维度的邻近点，导致插值引入空间信息，不符合仅用时间维度插值的要求：

import scipy.interpolate

# 原尝试代码（存在问题）
points = da.data.ravel()
valid = ~np.isnan(points)
points_valid = points[valid]

xx, yy, zz = np.meshgrid(x, y, z)
xx, yy, zz = xx.ravel(), yy.ravel(), zz.ravel()
xxv = xx[valid]
yyv = yy[valid]
zzv = zz[valid]

interpolated = scipy.interpolate.griddata(np.stack([xxv, yyv, zzv]).T, points_valid, (xx, yy, zz), method="linear")
da.values = interpolated.reshape(da.shape)

正确解决方案：xarray内置时间维度插值 #

使用xarray的interpolate_na方法，指定dim='time'后，会对每个(lat, lon)网格点的时间序列单独做线性插值，完全不涉及空间邻点：

import xarray as xr
import pandas as pd

# 构造xarray DataArray
da = xr.DataArray(
    test_data,
    dims=['time', 'lat', 'lon'],
    coords={
        'time': pd.date_range('2020-01-01', periods=3),  # 时间坐标需单调
        'lat': np.linspace(90, -90, test_data.shape[1]),
        'lon': np.linspace(-180, 180, test_data.shape[2])
    }
)

# 核心：仅针对time维度做线性插值
da_interpolated = da.interpolate_na(dim='time', method='linear')

# 查看中间缺失时间步的插值结果
print(da_interpolated.isel(time=1))

关键说明 #

• interpolate_na方法会自动遍历每个空间网格点，基于其前后时间步的有效数据做线性插值，严格限制在时间维度内计算。
• 要求时间坐标是单调递增/递减的（xarray默认生成的时间序列均满足此条件）。
• 相比scipy的全域插值，该方法更高效且完全符合需求，无需手动处理掩码或点云转换。

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