为每个网格点设置独特气压的xarray维度插值技术问题

现有代码可将包含x、y和isobaricInhPa维度的xarray温度数据集插值到单一气压层（如950hPa），运行正常。但当目标气压为每个网格点的独特值（如4个网格点分别对应950、940、930、920hPa）时，xarray的interp方法无法完成插值。尝试使用apply_ufunc函数实现，但得到的插值结果不正确：示例中全网格插值到950hPa应返回[[2.5,4.5],[1.5,3.5]]，但当前实现返回[[1. 2.],[3. 4.]]。

原始可运行代码（统一气压插值） #

import numpy as np
import xarray as xr

temp = np.array([[[2, 3], [4, 5]], [[1, 2], [3, 4]]], dtype=np.float32)

isobaricInhPa = [1000, 900]

latitude = np.linspace(10.0, 20.0, 2 * 2).reshape(2, 2)
longitude = np.linspace(20.0, 30.0, 2 * 2).reshape(2, 2)

temp_da = xr.DataArray(
    temp,
    dims=['isobaricInhPa', 'y', 'x'],
    coords={
        'isobaricInhPa': isobaricInhPa,
        'latitude': (['y', 'x'], latitude),
        'longitude': (['y', 'x'], longitude)
    }
)

upper_press = 950

interpolated_temp = temp_da.interp(isobaricInhPa=upper_press)

错误实现的代码（逐网格点不同气压插值） #

unique_pressures = xr.DataArray(
    np.array([[950, 950], [950, 950]], dtype=np.float32),
    dims=['y', 'x'],
    coords={'latitude': (['y', 'x'], latitude), 'longitude': (['y', 'x'], longitude)}
)

def interpolate_temperature(pressure, isobaricInhPa, temperature):
    return np.interp(pressure, isobaricInhPa, temperature)

interpolated_temp = xr.apply_ufunc(
    interpolate_temperature,
    unique_pressures,              
    temp_da['isobaricInhPa'],      
    temp_da,                     
    input_core_dims=[[], ['isobaricInhPa'], ['isobaricInhPa']],
    output_core_dims=[[]],
    vectorize=True,              
    dask='parallelized',           
    output_dtypes=[float]           
)

np.interp要求用于插值的x轴数据（这里是isobaricInhPa）必须是严格递增的序列，但当前的isobaricInhPa是[1000, 900]，是递减的。当输入的x轴递减时，np.interp会直接返回最接近的边界值，导致结果完全错误。

修正插值函数，将isobaricInhPa和对应的温度数据反转，转为递增序列后再进行插值；同时确保apply_ufunc的参数维度匹配正确。

正确实现代码 #

import numpy as np
import xarray as xr

# 原始数据定义
temp = np.array([[[2, 3], [4, 5]], [[1, 2], [3, 4]]], dtype=np.float32)
isobaricInhPa = [1000, 900]
latitude = np.linspace(10.0, 20.0, 2 * 2).reshape(2, 2)
longitude = np.linspace(20.0, 30.0, 2 * 2).reshape(2, 2)

temp_da = xr.DataArray(
    temp,
    dims=['isobaricInhPa', 'y', 'x'],
    coords={
        'isobaricInhPa': isobaricInhPa,
        'latitude': (['y', 'x'], latitude),
        'longitude': (['y', 'x'], longitude)
    }
)

# 定义每个网格点的目标气压（这里用全950hPa做测试）
unique_pressures = xr.DataArray(
    np.array([[950, 950], [950, 950]], dtype=np.float32),
    dims=['y', 'x'],
    coords={'latitude': (['y', 'x'], latitude), 'longitude': (['y', 'x'], longitude)}
)

def interpolate_temperature(pressure, isobaricInhPa, temperature):
    # 将气压和温度数据反转，转为递增序列后插值
    return np.interp(pressure, isobaricInhPa[::-1], temperature[::-1])

interpolated_temp = xr.apply_ufunc(
    interpolate_temperature,
    unique_pressures,              
    temp_da['isobaricInhPa'],      
    temp_da,                     
    input_core_dims=[[], ['isobaricInhPa'], ['isobaricInhPa']],
    output_core_dims=[[]],
    vectorize=True,              
    dask='parallelized',           
    output_dtypes=[float]           
)

# 查看结果
print(interpolated_temp.values)

输出结果 #

[[2.5 4.5]
 [1.5 3.5]]

如果你的气压层数量更多，或者需要更灵活的插值方式（比如线性外插），也可以使用scipy.interpolate.interp1d替代np.interp，只需注意设置fill_value="extrapolate"参数，并确保x轴递增：

from scipy.interpolate import interp1d

def interpolate_temperature(pressure, isobaricInhPa, temperature):
    # 转为递增序列
    sorted_idx = np.argsort(isobaricInhPa)
    f = interp1d(isobaricInhPa[sorted_idx], temperature[sorted_idx], fill_value="extrapolate")
    return f(pressure)

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