regridding strategies for comparison with observations

[JavierDiezSierra]

We are going to need to regrid CORDEX data (only the climatologies, not the raw data) to compare the results with observations/reanalysis. We implemented some Python functions for "conservative" remapping using the xESMF package within the framework of the Copernicus Atlas.

@larsbuntemeyer I mention this because I saw that you use the same Python package, and might be interesting.

@jesusff

[larsbuntemeyer]

Yes, with xESMF regridding becomes basically a one liner (if data complies with CF conventions). I think there is this rule of thumb also mentioned in https://doi.org/10.5194/gmd-7-1297-2014

In order to ensure a fair evaluation,
our strategy was to always use the coarser grid as reference,
except for mean sea-level pressure

which would mean for us:

• eobs (0.1°): regrid eobs to EURO-CORDEX grid (which is slightly coarser)
• ERA5: regrid EURO-CORDEX grid to ERA5

I think, the paper does not mention this, but in case of flux variables like pr we should use conservative remapping (i think in xesmf terms, it's "conservative_normed") and for grid point values like tas we can use bilinear interpolation.

and also:

Additionally, an elevation correction was carried
out for temperature assuming a uniform temperature lapse
rate of 0.0064 K m−1

see, e.g., https://github.com/euro-cordex/evaltools/blob/b42eb5a5228872f467b6530fd70fdf8684e66a4e/evaltools/eval.py#L89-L100

[larsbuntemeyer]

Yes, this is my prototype for comparison with eobs seasonal mean.

[JavierDiezSierra]

@larsbuntemeyer If we decide to use bilinear interpolation, it is a one-liner with xESMF, as you mentioned. However, for the conservative remapping method, it is necessary to calculate the lon and lat boundaries of the cells. For that, we use this approximation estimate_boundaries

[larsbuntemeyer]

Yes, we need bounds but they should be provided in the datasets as lon_vertices and lat_vertices which xESMF should understand. So we shouldnt' have to estimate them!

However, your code looks interesting, did you check cf.add_bounds? Should be interesting to compare!

For CORDEX rotated pole grid mappings, i have a function transform_bounds that can compute bounds exactly if they are missing using pyproj.

One problem i see is that ALADIN datasets don't have lon and lat bounds, but it should be possible also to compute them exactly using transform_bounds.

[larsbuntemeyer]

Here is a quick prototype how to handle bounds when regridding eobs conservatively in xesmf:

import xarray as xr
import numpy as np
import cf_xarray as cfxr
from xesmf import Regridder

store = "https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/pangeo-forge/EOBS-feedstock/eobs-tg-tn-tx-rr-hu-pp.zarr"
eobs = xr.open_dataset(store, engine="zarr", chunks={})
#eobs["mask"] = xr.where(~np.isnan(eobs.pp.isel(time=0)), 1, 0)

ds = xr.open_dataset("/mnt/CORDEX_CMIP6_tmp/sim_data/CORDEX/CMIP6/DD/EUR-12/GERICS/ERA5/evaluation/r1i1p1f1/REMO2020/v1-r1/mon/pr/v20241120/pr_EUR-12_ERA5_evaluation_r1i1p1f1_GERICS_REMO2020_v1-r1_mon_197901-198812.nc", decode_coords="all")

ds['lon_b'] = cfxr.bounds_to_vertices(
    ds.vertices_lon, bounds_dim="vertices", order="counterclockwise"
)
ds['lat_b'] = cfxr.bounds_to_vertices(
    ds.vertices_lat, bounds_dim="vertices", order="counterclockwise"
)

regridder = Regridder(eobs, ds, method="conservative_normed")
eobs_on_cordex_grid = regridder(eobs)
eobs_on_cordex_grid.pp.where(eobs_on_cordex_grid.pp>0).sel(time="2000-01-01T00:00:00").plot()

Note, that eobs has missing values over water so we should probably handle this using a mask for regridding.

[JavierDiezSierra]

Really nice, @larsbuntemeyer!

We tried to generalize the regridding process by calculating longitude and latitude boundaries on our own through an approximation because some curvilinear simulations didn't provide this information. However, I didn't know about the transform_bounds function. It is really interesting, and I will try it!

[larsbuntemeyer]

Yes, it's quite convenient! The core code ist based on pyproj.Transformer. So if you have native coordinates (e.g., in a rotated crs) and the grid mapping definition, you can use that to transform coordinates in any other CRS.

[larsbuntemeyer]

Note, that eobs has missing values over water so we should probably handle this using a mask for regridding.

Actually, we should use unmapped_to_nan=True to set eobs values that have no data on the cordex grid to nan. That's cleaner.

[gnikulin]

In general, just to keep in mind, remapping/regridding depends on direction, interpolation or aggregation. If we are going from coarser to finer grids (ERA5 -> EUR-12) or remapping between 2 grids with about the same resolution (EUR-12 and E-OBS), the bilinear interpolation should work fine. If we are going from finer to coarser grids (EUR-12 -> ERA5) we need to aggregate data. In this case, simple area-weighted remapping should work fine. The conservative remapping also works, although I'm not sure that we need to preserve fluxes for such kind of analysis.