基于Python+xarray的数据集插值后对比展示及精度保留问询

背景 #

原始数据集为GFS 2m温度场，网格步长0.25°×0.25°（纬度47.0°_{51.0°共17个点，经度1.0°}4.0°共13个点），通过xarray+cfgrib读取。需插值到目标坐标：纬度[48.5705, 50.4752]、经度[1.875, 3.75]，但遇到两个问题：无法直观对比原始与插值结果、插值后纬度坐标显示精度被截断。

1. 并排展示原始与插值数据集，验证插值正确性 #

可以通过Matplotlib绘制子图标记插值点，同时手动计算双线性插值结果与xarray插值结果对比，双重验证精度。

代码实现 #

import xarray as xr
import matplotlib.pyplot as plt
import numpy as np

# 读取原始数据
ds = xr.open_dataset('gfs.t00z.pgrb2.0p25.f000', engine='cfgrib')
t2m_original = ds['t2m']

# 定义目标插值坐标
target_lats = [48.5705, 50.4752]
target_lons = [1.875, 3.75]

# 执行插值
ds_interp = ds.interp(latitude=target_lats, longitude=target_lons, method="linear")
t2m_interp = ds_interp['t2m']

# 绘制并排对比图
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# 原始场热力图+插值点标记
t2m_original.plot(ax=ax1, cmap='viridis', extend='both')
ax1.scatter(target_lons, target_lats, color='#ff3333', s=120, marker='X', label='插值目标点')
ax1.set_title('原始2m温度场（标记插值目标点）')
ax1.legend(fontsize=10)

# 插值结果可视化+数值标注
t2m_interp.plot(ax=ax2, cmap='viridis', extend='both')
for lon, lat, val in zip(target_lons, target_lats, t2m_interp.values.flatten()):
    ax2.text(lon, lat, f'{val:.3f}', ha='center', va='bottom', color='white', fontsize=10)
ax2.set_title('插值后2m温度场')

plt.tight_layout()
plt.show()

# 手动计算双线性插值验证（以第一个点为例）
# 提取原始数据中相邻的网格点
lat_left, lat_right = 48.5, 48.75
lon_bottom, lon_top = 1.75, 2.0

# 获取四个角的温度值
t_bl = t2m_original.sel(latitude=lat_left, longitude=lon_bottom).values
t_br = t2m_original.sel(latitude=lat_left, longitude=lon_top).values
t_tl = t2m_original.sel(latitude=lat_right, longitude=lon_bottom).values
t_tr = t2m_original.sel(latitude=lat_right, longitude=lon_top).values

# 计算插值比例
lat_ratio = (target_lats[0] - lat_left) / (lat_right - lat_left)
lon_ratio = (target_lons[0] - lon_bottom) / (lon_top - lon_bottom)

# 双线性插值计算
manual_interp = (1 - lat_ratio)*(1 - lon_ratio)*t_bl + \
                (1 - lat_ratio)*lon_ratio*t_br + \
                lat_ratio*(1 - lon_ratio)*t_tl + \
                lat_ratio*lon_ratio*t_tr

# 对比结果
print(f"手动双线性插值结果：{manual_interp:.4f}")
print(f"xarray线性插值结果：{t2m_interp.sel(latitude=target_lats[0], longitude=target_lons[0]).values:.4f}")

验证逻辑 #

• 热力图对比：原始场标记插值点，直观确认插值点在网格中的位置；插值场标注数值，快速核对趋势合理性。
• 数值对比：手动双线性插值结果与xarray结果误差在浮点精度范围内，即可确认插值逻辑正确。

2. 保留插值坐标的精度 #

插值后坐标显示精度被截断是xarray默认格式化显示的问题，实际坐标值并未丢失精度，可通过以下方式处理：

查看原始精度坐标 #

直接读取坐标的numpy数组值，即可看到完整精度：

print("插值后纬度坐标实际值：", ds_interp.latitude.values)
print("插值后经度坐标实际值：", ds_interp.longitude.values)

输出会显示[48.5705 50.4752]，说明精度完全保留。

修改xarray显示精度 #

如果希望print(ds_interp['t2m'])时显示更高精度，设置xarray全局显示选项：

xr.set_options(display_precision=4)
print(ds_interp['t2m'])

此时输出会显示完整小数位，不再截断。

保存时确保精度 #

若需将插值结果保存到文件，指定用float64类型存储坐标，避免精度损失：

ds_interp.to_netcdf('interpolated_t2m.nc', encoding={
    'latitude': {'dtype': 'float64'},
    'longitude': {'dtype': 'float64'}
})

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