The domain module should give some tools to work with preconfigured or user defined domains. Domains are defined as xarray datasets that will contain dimensions and coodinates according to CF-conventions.
NOTE: The domain module mostly focuses on working with rotated cordex domains and how they are defined in the cordex archive specifications. However, there are some regional models that use different mappings instead of rotated_pole or rotated_latitude_longitude which we focus on. Any expertise working with those different mappings is highly welcome!
import cordex as cx
The domain module contains some useful functions to work with cordex meta data, e.g., you can get some domain grid information using
cx.domain_info("EUR-11")
The domain information is stored in a number of csv tables. The module contains a tables dictionary that sorts the tables by resolution or project, e.g.
cx.domains.tables.keys()
All available cordex domains are in those tables:
cx.domains.table
EUR-11 example¶The heart of the module are some functions that create a dataset from the grid information, e.g.
eur11 = cx.cordex_domain("EUR-11", dummy="topo")
eur11
The dummy='topo' argument means, we want a dummy variable in the dataset to see how the domain looks like. For the dummy topography, we use the cdo topo operator in the background. So maybe you have to install python-cdo, e.g., conda install -c conda-forge python-cdo. Working with xarray datasets means, that we can use all the nice functions of xarray including plotting. The latest versions of py-cordex also come with an xarray accessor that allows you, e.g., to get a quick overview of the domain, using, e.g.
# from matplotlib import pyplot as plt
# plt.figure(figsize=(8,8))
eur11.cx.map()
We can also use the dummy topography to plot an overview:
eur11.topo.plot(cmap="terrain")
eur11.topo.plot(x="lon", y="lat", cmap="terrain")
Let's define a slightly more sophisticated plotting function that uses cartopy for the right projection with a rotated pole:
def plot(da, pole, vmin=None, vmax=None, borders=True, title=None):
"""plot a domain using the right projection with cartopy"""
import cartopy.crs as ccrs
import cartopy.feature as cf
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 10))
projection = ccrs.PlateCarree()
transform = ccrs.RotatedPole(pole_latitude=pole[1], pole_longitude=pole[0])
# ax = plt.axes(projection=projection)
ax = plt.axes(projection=transform)
# ax.set_extent([ds_sub.rlon.min(), ds_sub.rlon.max(), ds_sub.rlat.min(), ds_sub.rlat.max()], crs=transform)
ax.gridlines(
draw_labels=True,
linewidth=0.5,
color="gray",
xlocs=range(-180, 180, 10),
ylocs=range(-90, 90, 5),
)
da.plot(
ax=ax,
cmap="terrain",
transform=transform,
vmin=vmin,
vmax=vmax,
x="rlon",
y="rlat",
)
ax.coastlines(resolution="50m", color="black", linewidth=1)
if borders:
ax.add_feature(cf.BORDERS)
if title is not None:
ax.set_title(title)
pole = (
eur11.rotated_latitude_longitude.grid_north_pole_longitude,
eur11.rotated_latitude_longitude.grid_north_pole_latitude,
)
pole
plot(eur11.topo, pole)
The domains are actually created from a csv table. The domains are created using the create_dataset function, in which the data from the table is fed, e.g.:
cx.create_dataset?
Let's create the EUR-11 domain manually from the numbers in the table:
cx.domains.table.loc["EUR-11"]
eur11_user = cx.create_dataset(
nlon=424,
nlat=412,
dlon=0.11,
dlat=0.11,
ll_lon=-28.375,
ll_lat=-23.375,
pollon=-162.00,
pollat=39.25,
dummy="topo",
)
We can check that this gives the same result as our preconfigured domain.
eur11_user.equals(eur11)
You can now use the create_dataset function to create any domain as an xarray dataset.
We need a slightly modified plotting routine for this:
def plots(dsets, vmin=None, vmax=None, borders=True, title=None):
"""plot a domain using the right projection with cartopy"""
import cartopy.crs as ccrs
import cartopy.feature as cf
import matplotlib.patheffects as pe
import matplotlib.pyplot as plt
plt.figure(figsize=(20, 10))
projection = ccrs.PlateCarree()
# transform = ccrs.RotatedPole(pole_latitude=pole[1], pole_longitude=pole[0])
# ax = plt.axes(projection=projection)
ax = plt.axes(projection=projection)
# ax.set_extent([ds_sub.rlon.min(), ds_sub.rlon.max(), ds_sub.rlat.min(), ds_sub.rlat.max()], crs=transform)
ax.gridlines(
draw_labels=True,
linewidth=0.5,
color="gray",
xlocs=range(-180, 180, 15),
ylocs=range(-90, 90, 10),
)
# path_effects = [pe.Stroke(linewidth=50, foreground='g'), pe.Normal()]
for ds in dsets:
pole = (
ds.rotated_latitude_longitude.grid_north_pole_longitude,
ds.rotated_latitude_longitude.grid_north_pole_latitude,
)
transform = ccrs.RotatedPole(pole_latitude=pole[1], pole_longitude=pole[0])
ds.topo.plot(
ax=ax,
cmap="terrain",
transform=transform,
vmin=vmin,
vmax=vmax,
x="rlon",
y="rlat",
add_colorbar=False,
)
ax.coastlines(resolution="50m", color="black", linewidth=1)
if borders:
ax.add_feature(cf.BORDERS)
if title is not None:
ax.set_title("")
Now, let's plot all cordex core domains into one overview. We will identify them by their resolution:
table = cx.domains.table
plots(
[
cx.cordex_domain(name, dummy="topo")
for name in table.index
if table.loc[name].dlon == 0.22
],
borders=False,
)
The gridded Cordex datasets are particular usefule for regridding either with CDO or other interpolation packages.
We will use some CMIP6 model data from the ESGF to show how we can do the regridding:
!wget https://raw.githubusercontent.com/remo-rcm/pyremo-data/main/orog_fx_MPI-ESM1-2-HR_historical_r1i1p1f1_gn.nc
import xarray as xr
# https://raw.githubusercontent.com/remo-rcm/pyremo-data/main/orog_fx_MPI-ESM1-2-HR_historical_r1i1p1f1_gn.nc
ds = xr.open_dataset("orog_fx_MPI-ESM1-2-HR_historical_r1i1p1f1_gn.nc")
# ds = xr.open_dataset(
# "https://esgf3.dkrz.de/thredds/dodsC/cmip6/CMIP/MPI-M/MPI-ESM1-2-HR/historical/r1i1p1f1/Amon/tas/gn/v20190710/tas_Amon_MPI-ESM1-2-HR_historical_r1i1p1f1_gn_199501-199912.nc"
# )
# ds = xr.open_dataset(
# "http://esgf-data3.ceda.ac.uk/thredds/dodsC/esg_cmip6/CMIP6/CMIP/MOHC/UKESM1-0-LL/historical/r1i1p1f2/Amon/tas/gn/v20190406/tas_Amon_UKESM1-0-LL_historical_r1i1p1f2_gn_185001-194912.nc"
# )
ds
CDO¶We will use CDO's python bindings to control the cdo operators. Please note, that the python bindings in the background actually execute the cdo commands on the command line. CDO does have a huge IO overhead since it will always write a file between each operator step and will always need data from a file on the filesystem as input. If you give an xarray dataset as input (see below), the python binding will actually trigger a to_netcdf call to write the input as a temporary file to disk. You should be aware of this if you use huge xarray datasets as input.
We will first write the EUR-11 grid into a file on the disk so that we can use it as input to CDO
from cdo import Cdo
eur11.to_netcdf("EUR-11_grid.nc")
Now we will remap the first timestep of the CMIP6 modeldata to the EUR-11 grid:
remap_cdo = Cdo().remapbil("EUR-11_grid.nc", input=ds, returnXArray="orog")
remap_cdo
remap_cdo.plot()
A nice alternative is xesmf since it is based on xarray and will also very nicely work with dask.
import xesmf as xe
regridder = xe.Regridder(ds, eur11, "bilinear", periodic=True)
regridder
remap_xe = regridder(ds.orog)
remap_xe.plot()
We can easily compare both approaches
(remap_cdo - remap_xe).plot(x="lon", y="lat")
The nice thing about xesmf is that it works together with xarray and will keep all meta information. Another consequence is that xesmf works well with dask and it's vectorization. That means, if we have a long time axis along we want to regrid, this can easily be parallelized using, e.g., dask.distributed and will also work nicely with large datasets that don't fit into memory. The xesmf regridder can also store and reuse regridding weights for faster interpolation.