如何在xarray中利用非索引坐标插值提取单点时间序列？

你加载的ARCO-ERA5数据集采用了扁平化存储结构：latitude和longitude并非独立维度，而是绑定在values维度上的坐标变量。这就是直接用.sel()或.interp()报错的核心原因——这两个方法要求目标坐标对应独立的数据维度。

以下是两种高效可行的解决方案：

方法一：提取最邻近格点 #

通过计算目标坐标与所有格点的距离，找到最接近的点的values索引，直接提取该索引对应的时间序列。这种方法计算量小，适合无需插值的场景：

import xarray as xr
import numpy as np

# 设定目标地理坐标
target_lat = 60
target_lon = 20

# 加载数据集（保留原分块配置）
data = xr.open_zarr(
    'gs://gcp-public-data-arco-era5/co/single-level-reanalysis.zarr/',
    chunks={'time': 48},
    consolidated=True,
)

# 计算每个格点与目标坐标的平面距离（高精度场景可替换为球面距离公式）
lat_diff = data.latitude - target_lat
lon_diff = data.longitude - target_lon
distance = np.sqrt(lat_diff**2 + lon_diff**2)

# 获取距离最小的格点索引
closest_idx = distance.argmin().compute()

# 提取目标点的CAPE时间序列
cape_time_series = data['cape'].sel(values=closest_idx)

方法二：对values维度执行插值 #

如果需要线性插值结果（而非最邻近点），可以用xarray.apply_ufunc对每个时间步的values维度执行插值操作，将latitude和longitude作为插值坐标传入：

from scipy.interpolate import interp1d

# 提前加载坐标数据（避免重复计算）
lats = data.latitude.compute()
lons = data.longitude.compute()

# 定义单时间步插值函数
def interpolate_point(arr, target_lat, target_lon):
    # 过滤无效数据点（可选）
    valid_mask = ~np.isnan(arr)
    valid_lats = lats[valid_mask]
    valid_lons = lons[valid_mask]
    valid_data = arr[valid_mask]
    
    # 基于经纬度执行线性插值，支持外推
    interp_func = interp1d(
        np.column_stack((valid_lats, valid_lons)),
        valid_data,
        kind='linear',
        fill_value='extrapolate'
    )
    return interp_func([[target_lat, target_lon]])[0]

# 对整个时间序列应用插值
cape_interpolated = xr.apply_ufunc(
    interpolate_point,
    data['cape'],
    input_core_dims=[['values']],
    output_core_dims=[[]],
    vectorize=True,
    kwargs={'target_lat': target_lat, 'target_lon': target_lon},
    dask='parallelized',
    output_dtypes=[np.float32]
)

注意事项 #

• 由于数据集采用Dask延迟加载，所有操作默认仅生成计算图，调用.compute()才会实际执行计算。
• 方法二中的插值方式可按需替换（例如用scipy.interpolate.griddata处理不规则格点）。
• 批量提取多点时，建议一次性整理所有目标坐标批量处理，提升计算效率。

内容的提问来源于stack exchange，提问作者Pörripeikko