使用 xarray interp 重新投影数据数组?

Boo*_*hin 1 raster map-projections python-xarray

我已经看了很多interp 函数的 xarray 文档,但我无法真正理解它。我看到这是一个重新投影,但它并不真正适合真实的案例。他们是否有人可以理解它,例如通过将这个数据集重新投影到 webmercator 数据上?

就像这个例子:

import xarray as xr
from pyproj import Transformer

ds = xr.tutorial.open_dataset("air_temperature").isel(time=0)
fig, axes = plt.subplots(ncols=2, figsize=(10, 4))
lon, lat = np.meshgrid(ds.lon, ds.lat)
shp = lon.shape
# reproject the grid
gcs_to_3857 = Transformer.from_crs(4326, 3857, always_xy=True)
x, y = gcs_to_3857.transform(lon.ravel(), lat.ravel())
# future index for a regular raster
X= np.linspace(x.min(), x.max(), shp[1])
Y= np.linspace(y.min(), y.max(), shp[0])   
data["x"] = xr.DataArray(np.reshape(x, shp), dims=("lat", "lon"))
data["y"] = xr.DataArray(np.reshape(y, shp), dims=("lat", "lon"))

Run Code Online (Sandbox Code Playgroud)

而在这里,我被困住了

应该是类似的东西ds.interp(x=X,y=Y),但数组是在经纬度上索引的

这对我来说有点令人困惑......

小智 5

您还可以按照建议使用重新投影rioxarray。这是代码:

import xarray as xr
import numpy as np
from rasterio.enums import Resampling
import matplotlib.pyplot as plt

ds = xr.tutorial.open_dataset('air_temperature').isel(time=0)
ds = ds.rio.write_crs('EPSG:4326')
dst = ds.rio.reproject('EPSG:3857', shape=(250, 250), resampling=Resampling.bilinear, nodata=np.nan)

fig, ax = plt.subplots(1, 2, figsize=(14, 5))
ds.air.plot(ax=ax[0])
dst.air.plot(ax=ax[1])
Run Code Online (Sandbox Code Playgroud)

在此输入图像描述