Using xarray interp to reproject a dataarray?

Question

I have looked a lot at the xarray documentation of the interp function and I cannot really make sense of it. I see it is a reprojection but it doesn't really fit a real case example. Is their someone that could make sense of it for example by reprojecting this dataset on a webmercator datum?

Something like the example:

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"))

And here, I am stuck

Should be something like ds.interp(x=X,y=Y) but the array is indexed on lat lon

It is all a bit confusing to me...

Answer 1

You could also use reproject from rioxarray as suggested. Here is the code:

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])