Cordex domains

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!

Working with domain information

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

{'short_name': 'EUR-11',
 'region': 4,
 'long_name': 'Europe',
 'nlon': 424,
 'nlat': 412,
 'll_lon': -28.375,
 'ur_lon': 18.155,
 'll_lat': -23.375,
 'ur_lat': 21.835,
 'dlon': 0.11,
 'dlat': 0.11,
 'pollon': -162.0,
 'pollat': 39.25}

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

dict_keys(['cordex', 'cordex-high-res', 'cordex-fps', 'cordex-core', 'cordex-regular'])

All available cordex domains are in those tables:

cx.domains.table

	region	long_name	nlon	nlat	ll_lon	ur_lon	ll_lat	ur_lat	dlon	dlat	pollon	pollat
short_name
SAM-44	1	South America	146	167	143.92000	207.72000	-38.28000	34.76000	0.4400	0.4400	-56.06	70.60
CAM-44	2	Central America	210	113	-52.80000	39.16000	-28.60000	20.68000	0.4400	0.4400	113.98	75.74
NAM-44	3	North America	155	130	-33.88000	33.88000	-28.40000	28.36000	0.4400	0.4400	83.00	42.50
EUR-44	4	Europe	106	103	-28.21000	17.99000	-23.21000	21.67000	0.4400	0.4400	-162.00	39.25
AFR-44	5	Africa	194	201	-24.64000	60.28000	-45.76000	42.24000	0.4400	0.4400	180.00	90.00
WAS-44	6	South Asia	193	130	-32.12000	52.36000	-21.56000	35.20000	0.4400	0.4400	-123.34	79.95
EAS-44	7	East Asia	203	167	-40.92000	47.96000	-26.84000	46.20000	0.4400	0.4400	-64.78	77.61
CAS-44	8	Central Asia	153	100	-34.32000	32.56000	-20.68000	22.88000	0.4400	0.4400	-103.39	43.48
AUS-44	9	Australasia	200	129	142.16000	229.72000	-22.88000	33.44000	0.4400	0.4400	141.38	60.31
ANT-44	10	Antarctica	125	97	152.72000	207.28000	-27.72000	14.52000	0.4400	0.4400	-166.92	6.08
ARC-44	11	Arctic	116	133	-22.88000	27.72000	-24.20000	33.88000	0.4400	0.4400	0.00	6.55
MED-44	12	Mediterranean	98	63	-23.22000	19.46000	-21.34000	5.94000	0.4400	0.4400	198.00	39.25
MNA-44	13	Middle East and North Africa	232	118	-26.40000	75.24000	-6.60000	44.88000	0.4400	0.4400	180.00	90.00
MNA-22	13	Middle East and North Africa	464	236	-26.51000	75.35000	-6.71000	44.99000	0.2200	0.2200	180.00	90.00
SAM-11	1	South America	584	668	143.75500	207.88500	-38.44500	34.92500	0.1100	0.1100	-56.06	70.60
CAM-11	2	Central America	840	452	-52.96500	39.32500	-28.76500	20.84500	0.1100	0.1100	113.98	75.74
NAM-11	3	North America	620	520	-34.04500	34.04500	-28.56500	28.52500	0.1100	0.1100	83.00	42.50
EUR-11	4	Europe	424	412	-28.37500	18.15500	-23.37500	21.83500	0.1100	0.1100	-162.00	39.25
AFR-11	5	Africa	776	804	-24.80500	60.44500	-45.92500	42.40500	0.1100	0.1100	180.00	90.00
WAS-11	6	South Asia	772	520	-32.28500	52.52500	-21.72500	35.36500	0.1100	0.1100	-123.34	79.95
EAS-11	7	East Asia	812	668	-41.08500	48.12500	-27.00500	46.36500	0.1100	0.1100	-64.78	77.61
CAS-11	8	Central Asia	612	400	-34.48500	32.72500	-20.84500	23.04500	0.1100	0.1100	-103.39	43.48
AUS-11	9	Australasia	800	516	141.99500	229.88500	-23.04500	33.60500	0.1100	0.1100	141.38	60.31
ANT-11	10	Antarctica	500	388	152.55500	207.44500	-27.88500	14.68500	0.1100	0.1100	-166.92	6.08
ARC-11	11	Arctic	464	532	-23.04500	27.88500	-24.36500	34.04500	0.1100	0.1100	0.00	6.55
MED-11	12	Mediterranean	392	252	-23.38500	19.62500	-21.50500	6.10500	0.1100	0.1100	198.00	39.25
MNA-11	13	Middle East and North Africa	928	472	-26.56500	75.40500	-6.76500	45.04500	0.1100	0.1100	180.00	90.00
GAR-0275	4	Greater Alpine Region	476	444	-13.51125	-0.44875	-10.93125	1.25125	0.0275	0.0275	-162.00	39.25
SAM-22	1	South America	292	334	143.81000	207.83000	-38.39000	34.87000	0.2200	0.2200	-56.06	70.60
CAM-22	2	Central America	420	226	-52.91000	39.27000	-28.71000	20.79000	0.2200	0.2200	113.98	75.74
NAM-22	3	North America	310	260	-33.99000	33.99000	-28.51000	28.47000	0.2200	0.2200	83.00	42.50
EUR-22	4	Europe	212	206	-28.32000	18.10000	-23.32000	21.78000	0.2200	0.2200	-162.00	39.25
AFR-22	5	Africa	388	402	-24.75000	60.39000	-45.87000	42.35000	0.2200	0.2200	180.00	90.00
WAS-22	6	South Asia	386	260	-32.23000	52.47000	-21.67000	35.31000	0.2200	0.2200	-123.34	79.95
EAS-22	7	East Asia	396	251	-43.23000	43.67000	-22.10000	32.90000	0.2200	0.2200	296.30	61.00
CAS-22	8	Central Asia	306	200	-34.43000	32.67000	-20.79000	22.99000	0.2200	0.2200	-103.39	43.48
AUS-22	9	Australasia	400	258	142.05000	229.83000	-22.99000	33.55000	0.2200	0.2200	141.38	60.31
SEA-22	14	South East Asia	264	194	89.26000	147.12000	-15.18000	27.28000	0.2200	0.2200	180.00	90.00
SAM-44i	1	South America	181	155	-106.25000	-16.25000	-58.25000	18.75000	0.5000	0.5000	NaN	NaN
CAM-44i	2	Central America	207	111	-124.75000	-21.75000	-19.75000	35.25000	0.5000	0.5000	NaN	NaN
NAM-44i	3	North America	300	129	-171.75000	-22.25000	12.25000	76.25000	0.5000	0.5000	NaN	NaN
EUR-44i	4	Europe	221	103	-44.75000	65.25000	21.75000	72.75000	0.5000	0.5000	NaN	NaN
AFR-44i	5	Africa	173	179	-25.25000	60.75000	-46.25000	42.75000	0.5000	0.5000	NaN	NaN
WAS-44i	6	South Asia	195	124	19.25000	116.25000	-15.75000	45.75000	0.5000	0.5000	NaN	NaN
EAS-44i	7	East Asia	227	157	62.75000	175.75000	-18.75000	59.25000	0.5000	0.5000	NaN	NaN
CAS-44i	8	Central Asia	260	133	10.75000	140.25000	17.75000	69.75000	0.5000	0.5000	NaN	NaN
AUS-44i	9	Australasia	238	133	88.75000	207.25000	-53.25000	12.75000	0.5000	0.5000	NaN	NaN
ANT-44i	10	Antarctica	720	70	-179.75000	179.75000	-89.75000	-55.25000	0.5000	0.5000	NaN	NaN
ARC-44i	11	Arctic	720	83	-179.75000	179.75000	48.75000	89.75000	0.5000	0.5000	NaN	NaN
MED-44i	12	Mediterranean	144	65	-20.75000	51.75000	25.25000	57.25000	0.5000	0.5000	NaN	NaN
MNA-44i	13	Middle East and North Africa	206	106	-26.75000	75.75000	-7.25000	45.25000	0.5000	0.5000	NaN	NaN
MNA-22i	13	Middle East and North Africa high res.	410	209	-26.62500	75.62500	-6.87500	45.12500	0.2500	0.2500	NaN	NaN
EUR-11i	4	Europe high res.	881	408	-44.81250	65.18750	21.81250	72.68750	0.1250	0.1250	NaN	NaN
SEA-22i	14	South East Asia	233	172	89.12500	147.12500	-15.37500	27.37500	0.2500	0.2500	NaN	NaN

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

<xarray.Dataset>
Dimensions:                     (rlon: 424, rlat: 412)
Coordinates:
  * rlon                        (rlon) float64 -28.38 -28.27 ... 18.05 18.16
  * rlat                        (rlat) float64 -23.38 -23.27 ... 21.73 21.84
    lon                         (rlat, rlon) float64 -10.06 -9.964 ... 64.96
    lat                         (rlat, rlon) float64 21.99 22.03 ... 66.75 66.69
Data variables:
    rotated_latitude_longitude  int32 0
    topo                        (rlat, rlon) float32 284.0 246.0 ... 509.0 509.0
Attributes:
    CORDEX_domain:  EUR-11

xarray.Dataset

• Dimensions:
  • rlon: 424
  • rlat: 412
• Coordinates: (4)
  • rlon
    (rlon)
    float64
    -28.38 -28.27 ... 18.05 18.16

    axis :

    X

    standard_name :

    grid_longitude

    long_name :

    longitude in rotated pole grid

    units :

    degrees

    array([-28.375, -28.265, -28.155, ...,  17.935,  18.045,  18.155])

  • rlat
    (rlat)
    float64
    -23.38 -23.27 ... 21.73 21.84

    axis :

    Y

    standard_name :

    grid_latitude

    long_name :

    latitude in rotated pole grid

    units :

    degrees

    array([-23.375, -23.265, -23.155, ...,  21.615,  21.725,  21.835])

  • lon
    (rlat, rlon)
    float64
    -10.06 -9.964 ... 64.76 64.96

    standard_name :

    longitude

    long_name :

    longitude

    units :

    degrees_east

    array([[-10.06387966,  -9.96388628,  -9.86380363, ...,  36.19872654,
             36.30631294,  36.41382968],
           [-10.11099589, -10.01086154,  -9.91063756, ...,  36.23077703,
             36.33853934,  36.44623165],
           [-10.15817955, -10.05790422,  -9.95753886, ...,  36.26288362,
             36.37082205,  36.47869015],
           ...,
           [-44.26744553, -44.1287729 , -43.98957131, ...,  64.21746364,
             64.42305922,  64.62774455],
           [-44.43028713, -44.29178077, -44.15274426, ...,  64.38488381,
             64.59063409,  64.79546746],
           [-44.59386389, -44.45552794, -44.31666067, ...,  64.55349992,
             64.75940009,  64.96437667]])

  • lat
    (rlat, rlon)
    float64
    21.99 22.03 22.07 ... 66.75 66.69

    standard_name :

    latitude

    long_name :

    latitude

    units :

    degrees_north

    array([[21.98782876, 22.02783807, 22.06771602, ..., 25.1677739 ,
            25.14109405, 25.11426239],
           [22.08878725, 22.1288583 , 22.16879786, ..., 25.27388366,
            25.24715843, 25.22028117],
           [22.18973223, 22.22986507, 22.2698663 , ..., 25.37998647,
            25.3532158 , 25.32629289],
           ...,
           [60.05553554, 60.13088567, 60.20613923, ..., 66.63852915,
            66.57692365, 66.51510791],
           [60.12974984, 60.20521265, 60.28057952, ..., 66.72632651,
            66.66454929, 66.60256233],
           [60.20376337, 60.27933874, 60.35481882, ..., 66.81394504,
            66.75199541, 66.68983654]])

• Data variables: (2)
  • rotated_latitude_longitude
    ()
    int32
    0

    grid_mapping_name :

    rotated_latitude_longitude

    grid_north_pole_latitude :

    39.25

    grid_north_pole_longitude :

    -162.0

    north_pole_grid_longitude :

    0.0

    array(0, dtype=int32)

  • topo
    (rlat, rlon)
    float32
    284.0 246.0 246.0 ... 509.0 509.0

    units :

    m

    grid_mapping :

    rotated_latitude_longitude

    array([[ 2.8400000e+02,  2.4600000e+02,  2.4600000e+02, ...,
            -1.6953334e+03, -1.6953334e+03, -1.6953334e+03],
           [ 2.8066666e+02,  2.8066666e+02,  2.4600000e+02, ...,
            -1.6953334e+03, -1.6953334e+03, -1.6953334e+03],
           [ 2.8066666e+02,  2.8066666e+02,  2.4600000e+02, ...,
            -1.6953334e+03, -1.6953334e+03, -1.6953334e+03],
           ...,
           [-6.6666669e-01, -6.6666669e-01,  4.6666665e+00, ...,
             5.3600000e+02,  5.3600000e+02,  5.0900000e+02],
           [-6.6666669e-01, -6.6666669e-01, -6.6666669e-01, ...,
             5.3600000e+02,  5.0900000e+02,  5.0900000e+02],
           [ 2.7966666e+02, -6.6666669e-01, -6.6666669e-01, ...,
             5.0900000e+02,  5.0900000e+02,  5.0900000e+02]], dtype=float32)

• Attributes: (1)

  CORDEX_domain :

  EUR-11

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, e.g.,

eur11.topo.plot(cmap="terrain")

<matplotlib.collections.QuadMesh at 0x1456ee850>

eur11.topo.plot(x="lon", y="lat", cmap="terrain")

<matplotlib.collections.QuadMesh at 0x1458587f0>

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"""
    %matplotlib inline
    import cartopy.crs as ccrs
    import cartopy.feature as cf
    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=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

(-162.0, 39.25)

plot(eur11.topo, pole)

User defined domain

The domains are actually created from csv tables. To checkout the tables you can have a look at dm.TABLES. This is a dictionary of dataframes created during the import of the model from a number of csv tables that define standard cordex domains. E.g., available tables are:

cx.domains.tables["cordex-high-res"]

	region	long_name	nlon	nlat	ll_lon	ur_lon	ll_lat	ur_lat	dlon	dlat	pollon	pollat
short_name
SAM-11	1	South America	584	668	143.755	207.885	-38.445	34.925	0.11	0.11	-56.06	70.60
CAM-11	2	Central America	840	452	-52.965	39.325	-28.765	20.845	0.11	0.11	113.98	75.74
NAM-11	3	North America	620	520	-34.045	34.045	-28.565	28.525	0.11	0.11	83.00	42.50
EUR-11	4	Europe	424	412	-28.375	18.155	-23.375	21.835	0.11	0.11	-162.00	39.25
AFR-11	5	Africa	776	804	-24.805	60.445	-45.925	42.405	0.11	0.11	180.00	90.00
WAS-11	6	South Asia	772	520	-32.285	52.525	-21.725	35.365	0.11	0.11	-123.34	79.95
EAS-11	7	East Asia	812	668	-41.085	48.125	-27.005	46.365	0.11	0.11	-64.78	77.61
CAS-11	8	Central Asia	612	400	-34.485	32.725	-20.845	23.045	0.11	0.11	-103.39	43.48
AUS-11	9	Australasia	800	516	141.995	229.885	-23.045	33.605	0.11	0.11	141.38	60.31
ANT-11	10	Antarctica	500	388	152.555	207.445	-27.885	14.685	0.11	0.11	-166.92	6.08
ARC-11	11	Arctic	464	532	-23.045	27.885	-24.365	34.045	0.11	0.11	0.00	6.55
MED-11	12	Mediterranean	392	252	-23.385	19.625	-21.505	6.105	0.11	0.11	198.00	39.25
MNA-11	13	Middle East and North Africa	928	472	-26.565	75.405	-6.765	45.045	0.11	0.11	180.00	90.00

The domains are created using the create_dataset function, e.g.:

help(cx.create_dataset)

Help on function create_dataset in module cordex.core.domain:

create_dataset(nlon, nlat, dlon, dlat, ll_lon, ll_lat, pollon=None, pollat=None, name=None, dummy=False, add_vertices=False, attrs=None, mapping_name=None, bounds=None, **kwargs)
    Create domain dataset from grid information.

    Parameters
    ----------
    nlon : int
        longitudal number of grid boxes
    nlat : int
        latitudal number of grid boxes
    dlon : float
        longitudal resolution (degrees)
    dlat : float
        latitudal resolution (degrees)
    ll_lon : float
        lower left longitude (degrees)
    ll_lat : float
        lower left latitude (degrees)
    pollon : float
        pol longitude (degrees)
    pollat : float
        pol latitude (degrees)
    dummy : str or logical
        Name of dummy field, if dummy=topo, the cdo topo operator will be
        used to create some dummy topography data. dummy data is useful for
        looking at the domain with ncview.
    add_vertices : bool
        Add grid boundaries in the global coordinates (lon_vertices and lat_vertices).
    attrs: str or dict
        Global attributes that should be added to the dataset. If `attrs='CORDEX'`
        a set of standard CF global attributes.
    mapping_name: str
        Variable name of the grid mapping, if mapping_name is `None`, the CF standard
        variable name is used.

Let’s create the EUR-11 domain manually from the numbers in the table:

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)

True

You can now use the create_dataset function to create any domain as an xarray dataset.

Check out the Africa domain!

afr11 = cx.cordex_domain("AFR-11", dummy="topo")
afr11

<xarray.Dataset>
Dimensions:                     (rlon: 776, rlat: 804)
Coordinates:
  * rlon                        (rlon) float64 -24.8 -24.7 ... 60.34 60.45
  * rlat                        (rlat) float64 -45.92 -45.81 ... 42.3 42.41
    lon                         (rlat, rlon) float64 -24.81 -24.7 ... 60.44
    lat                         (rlat, rlon) float64 -45.92 -45.92 ... 42.41
Data variables:
    rotated_latitude_longitude  int32 0
    topo                        (rlat, rlon) float32 -4.429e+03 ... 161.0
Attributes:
    CORDEX_domain:  AFR-11

xarray.Dataset

• Dimensions:
  • rlon: 776
  • rlat: 804
• Coordinates: (4)
  • rlon
    (rlon)
    float64
    -24.8 -24.7 -24.59 ... 60.34 60.45

    axis :

    X

    standard_name :

    grid_longitude

    long_name :

    longitude in rotated pole grid

    units :

    degrees

    array([-24.805, -24.695, -24.585, ...,  60.225,  60.335,  60.445])

  • rlat
    (rlat)
    float64
    -45.92 -45.81 -45.7 ... 42.3 42.41

    axis :

    Y

    standard_name :

    grid_latitude

    long_name :

    latitude in rotated pole grid

    units :

    degrees

    array([-45.925, -45.815, -45.705, ...,  42.185,  42.295,  42.405])

  • lon
    (rlat, rlon)
    float64
    -24.81 -24.7 -24.59 ... 60.33 60.44

    standard_name :

    longitude

    long_name :

    longitude

    units :

    degrees_east

    array([[-24.805, -24.695, -24.585, ...,  60.225,  60.335,  60.445],
           [-24.805, -24.695, -24.585, ...,  60.225,  60.335,  60.445],
           [-24.805, -24.695, -24.585, ...,  60.225,  60.335,  60.445],
           ...,
           [-24.805, -24.695, -24.585, ...,  60.225,  60.335,  60.445],
           [-24.805, -24.695, -24.585, ...,  60.225,  60.335,  60.445],
           [-24.805, -24.695, -24.585, ...,  60.225,  60.335,  60.445]])

  • lat
    (rlat, rlon)
    float64
    -45.92 -45.92 ... 42.41 42.41

    standard_name :

    latitude

    long_name :

    latitude

    units :

    degrees_north

    array([[-45.925, -45.925, -45.925, ..., -45.925, -45.925, -45.925],
           [-45.815, -45.815, -45.815, ..., -45.815, -45.815, -45.815],
           [-45.705, -45.705, -45.705, ..., -45.705, -45.705, -45.705],
           ...,
           [ 42.185,  42.185,  42.185, ...,  42.185,  42.185,  42.185],
           [ 42.295,  42.295,  42.295, ...,  42.295,  42.295,  42.295],
           [ 42.405,  42.405,  42.405, ...,  42.405,  42.405,  42.405]])

• Data variables: (2)
  • rotated_latitude_longitude
    ()
    int32
    0

    grid_mapping_name :

    rotated_latitude_longitude

    grid_north_pole_latitude :

    90.0

    grid_north_pole_longitude :

    180.0

    north_pole_grid_longitude :

    0.0

    array(0, dtype=int32)

  • topo
    (rlat, rlon)
    float32
    -4.429e+03 -4.429e+03 ... 161.0

    units :

    m

    grid_mapping :

    rotated_latitude_longitude

    array([[-4429.    , -4429.    , -4429.    , ..., -4703.6665, -4703.6665,
            -4703.6665],
           [-4429.    , -4429.    , -4429.    , ..., -4703.6665, -4703.6665,
            -4703.6665],
           [-4429.    , -4429.    , -4429.    , ..., -4703.6665, -4703.6665,
            -4703.6665],
           ...,
           [-3270.    , -3270.    , -3270.    , ...,   161.    ,   161.    ,
              161.    ],
           [-3270.    , -3270.    , -3270.    , ...,   161.    ,   161.    ,
              161.    ],
           [-3270.    , -3270.    , -3270.    , ...,   161.    ,   161.    ,
              161.    ]], dtype=float32)

• Attributes: (1)

  CORDEX_domain :

  AFR-11

pole = (
    afr11.rotated_latitude_longitude.grid_north_pole_longitude,
    afr11.rotated_latitude_longitude.grid_north_pole_latitude,
)
pole

(180.0, 90.0)

plot(afr11.topo, pole)

Plot all cordex-core domains

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"""
    %matplotlib inline
    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:

plots(
    [
        cx.cordex_domain(name, dummy="topo")
        for name in cx.domains.tables["cordex-core"].index
    ],
    borders=False,
)

Using the grid information for interpolation

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:

import xarray as xr

# 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

<xarray.Dataset>
Dimensions:    (time: 1200, bnds: 2, lat: 144, lon: 192)
Coordinates:
  * time       (time) object 1850-01-16 00:00:00 ... 1949-12-16 00:00:00
  * lat        (lat) float64 -89.38 -88.12 -86.88 -85.62 ... 86.88 88.12 89.38
  * lon        (lon) float64 0.9375 2.812 4.688 6.562 ... 355.3 357.2 359.1
    height     float64 1.5
Dimensions without coordinates: bnds
Data variables:
    time_bnds  (time, bnds) object 1850-01-01 00:00:00 ... 1950-01-01 00:00:00
    lat_bnds   (lat, bnds) float64 -90.0 -88.75 -88.75 ... 88.75 88.75 90.0
    lon_bnds   (lon, bnds) float64 0.0 1.875 1.875 3.75 ... 358.1 358.1 360.0
    tas        (time, lat, lon) float32 ...
Attributes: (12/47)
    Conventions:                     CF-1.7 CMIP-6.2
    activity_id:                     CMIP
    branch_method:                   standard
    branch_time_in_child:            0.0
    branch_time_in_parent:           144000.0
    creation_date:                   2019-04-05T16:02:56Z
    ...                              ...
    variable_id:                     tas
    variant_label:                   r1i1p1f2
    license:                         CMIP6 model data produced by the Met Off...
    cmor_version:                    3.4.0
    tracking_id:                     hdl:21.14100/255d149c-12fc-41f1-878d-034...
    DODS_EXTRA.Unlimited_Dimension:  time

xarray.Dataset

• Dimensions:
  • time: 1200
  • bnds: 2
  • lat: 144
  • lon: 192
• Coordinates: (4)
  • time
    (time)
    object
    1850-01-16 00:00:00 ... 1949-12-...

    bounds :

    time_bnds

    axis :

    T

    long_name :

    time

    standard_name :

    time

    _ChunkSizes :

    524288

    array([cftime.Datetime360Day(1850, 1, 16, 0, 0, 0, 0, has_year_zero=True),
           cftime.Datetime360Day(1850, 2, 16, 0, 0, 0, 0, has_year_zero=True),
           cftime.Datetime360Day(1850, 3, 16, 0, 0, 0, 0, has_year_zero=True), ...,
           cftime.Datetime360Day(1949, 10, 16, 0, 0, 0, 0, has_year_zero=True),
           cftime.Datetime360Day(1949, 11, 16, 0, 0, 0, 0, has_year_zero=True),
           cftime.Datetime360Day(1949, 12, 16, 0, 0, 0, 0, has_year_zero=True)],
          dtype=object)

  • lat
    (lat)
    float64
    -89.38 -88.12 ... 88.12 89.38

    bounds :

    lat_bnds

    units :

    degrees_north

    axis :

    Y

    long_name :

    Latitude

    standard_name :

    latitude

    array([-89.375, -88.125, -86.875, -85.625, -84.375, -83.125, -81.875, -80.625,
           -79.375, -78.125, -76.875, -75.625, -74.375, -73.125, -71.875, -70.625,
           -69.375, -68.125, -66.875, -65.625, -64.375, -63.125, -61.875, -60.625,
           -59.375, -58.125, -56.875, -55.625, -54.375, -53.125, -51.875, -50.625,
           -49.375, -48.125, -46.875, -45.625, -44.375, -43.125, -41.875, -40.625,
           -39.375, -38.125, -36.875, -35.625, -34.375, -33.125, -31.875, -30.625,
           -29.375, -28.125, -26.875, -25.625, -24.375, -23.125, -21.875, -20.625,
           -19.375, -18.125, -16.875, -15.625, -14.375, -13.125, -11.875, -10.625,
            -9.375,  -8.125,  -6.875,  -5.625,  -4.375,  -3.125,  -1.875,  -0.625,
             0.625,   1.875,   3.125,   4.375,   5.625,   6.875,   8.125,   9.375,
            10.625,  11.875,  13.125,  14.375,  15.625,  16.875,  18.125,  19.375,
            20.625,  21.875,  23.125,  24.375,  25.625,  26.875,  28.125,  29.375,
            30.625,  31.875,  33.125,  34.375,  35.625,  36.875,  38.125,  39.375,
            40.625,  41.875,  43.125,  44.375,  45.625,  46.875,  48.125,  49.375,
            50.625,  51.875,  53.125,  54.375,  55.625,  56.875,  58.125,  59.375,
            60.625,  61.875,  63.125,  64.375,  65.625,  66.875,  68.125,  69.375,
            70.625,  71.875,  73.125,  74.375,  75.625,  76.875,  78.125,  79.375,
            80.625,  81.875,  83.125,  84.375,  85.625,  86.875,  88.125,  89.375])

  • lon
    (lon)
    float64
    0.9375 2.812 4.688 ... 357.2 359.1

    bounds :

    lon_bnds

    units :

    degrees_east

    axis :

    X

    long_name :

    Longitude

    standard_name :

    longitude

    array([  0.9375,   2.8125,   4.6875,   6.5625,   8.4375,  10.3125,  12.1875,
            14.0625,  15.9375,  17.8125,  19.6875,  21.5625,  23.4375,  25.3125,
            27.1875,  29.0625,  30.9375,  32.8125,  34.6875,  36.5625,  38.4375,
            40.3125,  42.1875,  44.0625,  45.9375,  47.8125,  49.6875,  51.5625,
            53.4375,  55.3125,  57.1875,  59.0625,  60.9375,  62.8125,  64.6875,
            66.5625,  68.4375,  70.3125,  72.1875,  74.0625,  75.9375,  77.8125,
            79.6875,  81.5625,  83.4375,  85.3125,  87.1875,  89.0625,  90.9375,
            92.8125,  94.6875,  96.5625,  98.4375, 100.3125, 102.1875, 104.0625,
           105.9375, 107.8125, 109.6875, 111.5625, 113.4375, 115.3125, 117.1875,
           119.0625, 120.9375, 122.8125, 124.6875, 126.5625, 128.4375, 130.3125,
           132.1875, 134.0625, 135.9375, 137.8125, 139.6875, 141.5625, 143.4375,
           145.3125, 147.1875, 149.0625, 150.9375, 152.8125, 154.6875, 156.5625,
           158.4375, 160.3125, 162.1875, 164.0625, 165.9375, 167.8125, 169.6875,
           171.5625, 173.4375, 175.3125, 177.1875, 179.0625, 180.9375, 182.8125,
           184.6875, 186.5625, 188.4375, 190.3125, 192.1875, 194.0625, 195.9375,
           197.8125, 199.6875, 201.5625, 203.4375, 205.3125, 207.1875, 209.0625,
           210.9375, 212.8125, 214.6875, 216.5625, 218.4375, 220.3125, 222.1875,
           224.0625, 225.9375, 227.8125, 229.6875, 231.5625, 233.4375, 235.3125,
           237.1875, 239.0625, 240.9375, 242.8125, 244.6875, 246.5625, 248.4375,
           250.3125, 252.1875, 254.0625, 255.9375, 257.8125, 259.6875, 261.5625,
           263.4375, 265.3125, 267.1875, 269.0625, 270.9375, 272.8125, 274.6875,
           276.5625, 278.4375, 280.3125, 282.1875, 284.0625, 285.9375, 287.8125,
           289.6875, 291.5625, 293.4375, 295.3125, 297.1875, 299.0625, 300.9375,
           302.8125, 304.6875, 306.5625, 308.4375, 310.3125, 312.1875, 314.0625,
           315.9375, 317.8125, 319.6875, 321.5625, 323.4375, 325.3125, 327.1875,
           329.0625, 330.9375, 332.8125, 334.6875, 336.5625, 338.4375, 340.3125,
           342.1875, 344.0625, 345.9375, 347.8125, 349.6875, 351.5625, 353.4375,
           355.3125, 357.1875, 359.0625])

  • height
    ()
    float64
    ...

    units :

    m

    axis :

    Z

    positive :

    up

    long_name :

    height

    standard_name :

    height

    array(1.5)

• Data variables: (4)
  • time_bnds
    (time, bnds)
    object
    ...

    _ChunkSizes :

    [1 2]

    array([[cftime.Datetime360Day(1850, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.Datetime360Day(1850, 2, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.Datetime360Day(1850, 2, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.Datetime360Day(1850, 3, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.Datetime360Day(1850, 3, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.Datetime360Day(1850, 4, 1, 0, 0, 0, 0, has_year_zero=True)],
           ...,
           [cftime.Datetime360Day(1949, 10, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.Datetime360Day(1949, 11, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.Datetime360Day(1949, 11, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.Datetime360Day(1949, 12, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.Datetime360Day(1949, 12, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.Datetime360Day(1950, 1, 1, 0, 0, 0, 0, has_year_zero=True)]],
          dtype=object)

  • lat_bnds
    (lat, bnds)
    float64
    ...

    _ChunkSizes :

    [144 2]

    array([[-90.  , -88.75],
           [-88.75, -87.5 ],
           [-87.5 , -86.25],
           ...,
           [ 86.25,  87.5 ],
           [ 87.5 ,  88.75],
           [ 88.75,  90.  ]])

  • lon_bnds
    (lon, bnds)
    float64
    ...

    _ChunkSizes :

    [192 2]

    array([[  0.   ,   1.875],
           [  1.875,   3.75 ],
           [  3.75 ,   5.625],
           ...,
           [354.375, 356.25 ],
           [356.25 , 358.125],
           [358.125, 360.   ]])

  • tas
    (time, lat, lon)
    float32
    ...

    standard_name :

    air_temperature

    long_name :

    Near-Surface Air Temperature

    comment :

    near-surface (usually, 2 meter) air temperature

    units :

    K

    original_name :

    mo: (stash: m01s03i236, lbproc: 128)

    cell_methods :

    area: time: mean

    cell_measures :

    area: areacella

    history :

    2019-04-05T16:02:56Z altered by CMOR: Treated scalar dimension: 'height'. 2019-04-05T16:02:56Z altered by CMOR: replaced missing value flag (-1.07374e+09) with standard missing value (1e+20).

    _ChunkSizes :

    [ 1 144 192]

    [33177600 values with dtype=float32]

• Attributes: (47)

  Conventions :

  CF-1.7 CMIP-6.2

  activity_id :

  CMIP

  branch_method :

  standard

  branch_time_in_child :

  0.0

  branch_time_in_parent :

  144000.0

  creation_date :

  2019-04-05T16:02:56Z

  cv_version :

  6.2.20.1

  data_specs_version :

  01.00.29

  experiment :

  all-forcing simulation of the recent past

  experiment_id :

  historical

  external_variables :

  areacella

  forcing_index :

  2

  frequency :

  mon

  further_info_url :

  https://furtherinfo.es-doc.org/CMIP6.MOHC.UKESM1-0-LL.historical.none.r1i1p1f2

  grid :

  Native N96 grid; 192 x 144 longitude/latitude

  grid_label :

  gn

  history :

  2019-04-05T15:50:03Z ; CMOR rewrote data to be consistent with CMIP6, CF-1.7 CMIP-6.2 and CF standards.; 2019-04-05T15:49:43Z MIP Convert v1.0.2, Python v2.7.12, Iris v1.13.0, Numpy v1.13.3, netcdftime v1.4.1.

  initialization_index :

  1

  institution :

  Met Office Hadley Centre, Fitzroy Road, Exeter, Devon, EX1 3PB, UK

  institution_id :

  MOHC

  mip_era :

  CMIP6

  mo_runid :

  u-bc179

  nominal_resolution :

  250 km

  parent_activity_id :

  CMIP

  parent_experiment_id :

  piControl

  parent_mip_era :

  CMIP6

  parent_source_id :

  UKESM1-0-LL

  parent_time_units :

  days since 1850-01-01-00-00-00

  parent_variant_label :

  r1i1p1f2

  physics_index :

  1

  product :

  model-output

  realization_index :

  1

  realm :

  atmos

  source :

  UKESM1.0-LL (2018): aerosol: UKCA-GLOMAP-mode atmos: MetUM-HadGEM3-GA7.1 (N96; 192 x 144 longitude/latitude; 85 levels; top level 85 km) atmosChem: UKCA-StratTrop land: JULES-ES-1.0 landIce: none ocean: NEMO-HadGEM3-GO6.0 (eORCA1 tripolar primarily 1 deg with meridional refinement down to 1/3 degree in the tropics; 360 x 330 longitude/latitude; 75 levels; top grid cell 0-1 m) ocnBgchem: MEDUSA2 seaIce: CICE-HadGEM3-GSI8 (eORCA1 tripolar primarily 1 deg; 360 x 330 longitude/latitude)

  source_id :

  UKESM1-0-LL

  source_type :

  AOGCM AER BGC CHEM

  sub_experiment :

  none

  sub_experiment_id :

  none

  table_id :

  Amon

  table_info :

  Creation Date:(13 December 2018) MD5:2b12b5db6db112aa8b8b0d6c1645b121

  title :

  UKESM1-0-LL output prepared for CMIP6

  variable_id :

  tas

  variant_label :

  r1i1p1f2

  license :

  CMIP6 model data produced by the Met Office Hadley Centre is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License (https://creativecommons.org/licenses). Consult https://pcmdi.llnl.gov/CMIP6/TermsOfUse for terms of use governing CMIP6 output, including citation requirements and proper acknowledgment. Further information about this data, including some limitations, can be found via the further_info_url (recorded as a global attribute in this file) and at https://ukesm.ac.uk/cmip6. The data producers and data providers make no warranty, either express or implied, including, but not limited to, warranties of merchantability and fitness for a particular purpose. All liabilities arising from the supply of the information (including any liability arising in negligence) are excluded to the fullest extent permitted by law.

  cmor_version :

  3.4.0

  tracking_id :

  hdl:21.14100/255d149c-12fc-41f1-878d-0344c8b797f4

  DODS_EXTRA.Unlimited_Dimension :

  time

Interpolation using 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.isel(time=0), returnXArray="tas")
remap_cdo

<xarray.DataArray 'tas' (rlat: 412, rlon: 424)>
[174688 values with dtype=float32]
Coordinates:
    lon      (rlat, rlon) float64 ...
    lat      (rlat, rlon) float64 ...
  * rlon     (rlon) float64 -28.38 -28.27 -28.16 -28.05 ... 17.94 18.05 18.16
  * rlat     (rlat) float64 -23.38 -23.27 -23.16 -23.05 ... 21.62 21.73 21.84
    height   float64 1.5
Attributes:
    standard_name:  air_temperature
    long_name:      Near-Surface Air Temperature
    units:          K
    grid_mapping:   rotated_latitude_longitude
    comment:        near-surface (usually, 2 meter) air temperature
    original_name:  mo: (stash: m01s03i236, lbproc: 128)
    cell_methods:   area: time: mean
    cell_measures:  area: areacella
    history:        2019-04-05T16:02:56Z altered by CMOR: Treated scalar dime...
    _ChunkSizes:    [  1 144 192]

xarray.DataArray

'tas'

• rlat: 412
• rlon: 424

• ...

  [174688 values with dtype=float32]

• Coordinates: (5)
  • lon
    (rlat, rlon)
    float64
    ...

    standard_name :

    longitude

    long_name :

    longitude

    units :

    degrees_east

    _CoordinateAxisType :

    Lon

    [174688 values with dtype=float64]

  • lat
    (rlat, rlon)
    float64
    ...

    standard_name :

    latitude

    long_name :

    latitude

    units :

    degrees_north

    _CoordinateAxisType :

    Lat

    [174688 values with dtype=float64]

  • rlon
    (rlon)
    float64
    -28.38 -28.27 ... 18.05 18.16

    standard_name :

    grid_longitude

    long_name :

    longitude in rotated pole grid

    units :

    degrees

    axis :

    X

    array([-28.375, -28.265, -28.155, ...,  17.935,  18.045,  18.155])

  • rlat
    (rlat)
    float64
    -23.38 -23.27 ... 21.73 21.84

    standard_name :

    grid_latitude

    long_name :

    latitude in rotated pole grid

    units :

    degrees

    axis :

    Y

    array([-23.375, -23.265, -23.155, ...,  21.615,  21.725,  21.835])

  • height
    ()
    float64
    ...

    standard_name :

    height

    long_name :

    height

    units :

    m

    positive :

    up

    axis :

    Z

    array(1.5)

• Attributes: (10)

  standard_name :

  air_temperature

  long_name :

  Near-Surface Air Temperature

  units :

  K

  grid_mapping :

  rotated_latitude_longitude

  comment :

  near-surface (usually, 2 meter) air temperature

  original_name :

  mo: (stash: m01s03i236, lbproc: 128)

  cell_methods :

  area: time: mean

  cell_measures :

  area: areacella

  history :

  2019-04-05T16:02:56Z altered by CMOR: Treated scalar dimension: 'height'. 2019-04-05T16:02:56Z altered by CMOR: replaced missing value flag (-1.07374e+09) with standard missing value (1e+20).

  _ChunkSizes :

  [ 1 144 192]

(remap_cdo - 273.5).plot()

<matplotlib.collections.QuadMesh at 0x14696f820>

Alternative interpolation methods

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

xESMF Regridder
Regridding algorithm:       bilinear
Weight filename:            bilinear_144x192_412x424_peri.nc
Reuse pre-computed weights? False
Input grid shape:           (144, 192)
Output grid shape:          (412, 424)
Periodic in longitude?      True

remap_xe = regridder(ds.tas.isel(time=0))

(remap_xe - 273.5).plot()

<matplotlib.collections.QuadMesh at 0x1184813d0>

We can easily compare both approaches

(remap_cdo - remap_xe).plot(x="lon", y="lat")

<matplotlib.collections.QuadMesh at 0x1460e6a30>

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.