xarray的interp最近方法如何选择最近的中心？

编程问答 2022-05-28

问题描述

我有一个二维 xarray 数据集，我想在 lon 和 lot 坐标上进行插值，以便获得更高的分辨率，但这些值与每个坐标处的原始值完全对应。 我认为优秀的 xr.interp 函数能够做到这一点，但在 the example 之后，我发现原始值和内插值之间存在一些差异。我将经度和纬度分辨率增加 4，因此除了在原始数据集中出现一次的所有 air 值之外，在内插数据集中出现 16 次，但事实并非如此。

有谁知道原始数据集和插值数据集不对齐的原因是什么，我该如何解决？

ds = xr.tutorial.open_dataset("air_temperature").isel(time=0)
fig,axes = plt.subplots(ncols=2,figsize=(10,4))
ds_sel=ds.sel(lon=slice(250,260),lat=slice(40,30))
ds.air.plot(ax=axes[0],xlim=(250,ylim=(30,40))

axes[0].set_title("Raw data")

# Interpolated data
new_lon = np.linspace(ds.lon[0],ds.lon[-1],ds.dims["lon"] * 4)
new_lat = np.linspace(ds.lat[0],ds.lat[-1],ds.dims["lat"] * 4)

dsi = ds.interp(lat=new_lat,lon=new_lon,method="nearest")
dsi_sel=dsi.sel(lon=slice(250,30))
dsi.air.plot(ax=axes[1],40))


axes[1].set_title("Interpolated data")

Link

用

显示唯一值 
unique,counts = np.unique(ds_sel.air.values,return_counts=True)
print("original values",dict(zip(unique,counts)))
unique,counts = np.unique(dsi_sel.air.values,return_counts=True)
print("interpolated values",counts)))

我明白

original values {262.1: 1,263.1: 1,263.9: 1,264.4: 1,265.19998: 1,266.6: 1,266.79: 1,266.9: 2,268.29: 1,269.79: 1,270.4: 1,273.0: 1,273.6: 1,275.19998: 1,276.29: 1,278.0: 1,278.5: 1,278.6: 1,281.5: 1,282.1: 1,282.29: 1,284.6: 1,286.79: 1,288.0: 1}

interpolated values {262.1: 4,263.1: 8,263.9: 8,264.4: 8,265.19998: 4,266.6: 16,266.79: 16,266.9: 24,268.29: 8,269.79: 20,270.4: 10,273.0: 20,273.6: 16,275.19998: 8,276.29: 20,278.0: 16,278.5: 10,278.6: 8,281.5: 4,282.1: 16,282.29: 8,284.6: 8,286.79: 8,288.0: 4}

解决方法

我认为您在概念上遇到了围栏错误（请参阅本页上的部分：https://en.wikipedia.org/wiki/Off-by-one_error）

您应该将 xarray 坐标解释为“中点”，而不是单元格边界。

您的 new_lon 没有很好地分为 1/2、1/4、1/8 等：

print(new_lon)
[200.         200.61611374 201.23222749 201.84834123 202.46445498
 203.08056872 203.69668246 204.31279621 204.92890995]

并且它不包括所有原始坐标。

考虑到“逐个”：

new_lon = np.linspace(ds.lon[0],ds.lon[-1],(ds.dims["lon"] - 1) * 4 + 1)
new_lat = np.linspace(ds.lat[0],ds.lat[-1],(ds.dims["lat"] - 1) * 4 + 1)
print(new_lon)
[200.    200.625 201.25  201.875 202.5   203.125 203.75  204.375 205.   ]

然后你可以例如检查原始第一行和插值行的部分：

selection = ds["air"][0,:3]
selection_i = dsi["air"][0,:9]
print(selection["lon"])
print(selection.values)
print(selection_i["lon"])
print(selection_i.values)

这对我来说很好：

[200.  202.5 205. ]
[241.2 242.5 243.5]
[200.    200.625 201.25  201.875 202.5   203.125 203.75  204.375 205.   ]
[241.2   241.2   241.2   242.5   242.5   242.5   242.5   243.5   243.5]

当然，在进行最近插值时，您可能会得到平局： 0.5 与 0.0 的距离与 1.0 的距离同样远——因此您不得不在不经意间偏向“向上”或“向下”以获得一个最接近的值。

另请注意，绘制 Matplotlib QuadMesh 的 .plot() 命令必须以某种方式从中点推断边界。这有时会导致绘制的边界与您可能想到的略有不同（尤其是在坐标间隔不均匀的情况下）。

interpolationpythonpython-xarray