This notebook is supposed to support a question concerning the order of vertices (lat lon bounds) in curvilinear grids and how it is handled in cdo
and the CF-conventions.
import numpy as np
import xarray as xr
from cartopy import crs as ccrs
from matplotlib import pyplot as plt
We want to plot cell vertices of a cmorized file created using the cdo cmor
operator. During the cmorization, cdo cmor
adds cell vertices for curvilinear grids (e.g. rotated pole):
ds = xr.open_dataset(
"/mnt/lustre02/work/ch0636/g300046/ESGF/euro-cordex/CDO_CMOR_RESULTS/CORDEX/output/EUR-11/GERICS/ECMWF-ERAINT/evaluation/r0i0p0/GERICS-REMO2015/v1/fx/orog/orog_EUR-11_ECMWF-ERAINT_evaluation_r0i0p0_GERICS-REMO2015_v1_fx_r0i0p0.nc"
)
pole = (
ds.rotated_latitude_longitude.grid_north_pole_longitude,
ds.rotated_latitude_longitude.grid_north_pole_latitude,
)
subset = (ds).isel(rlon=slice(0, 3), rlat=slice(0, 3))
subset
<xarray.Dataset> Dimensions: (rlat: 3, rlon: 3, vertices: 4) Coordinates: * rlat (rlat) float64 -23.38 -23.27 -23.16 * rlon (rlon) float64 -28.38 -28.27 -28.16 lat (rlat, rlon) float32 21.99 22.03 ... 22.23 22.27 lon (rlat, rlon) float32 349.9 350.0 ... 349.9 350.0 Dimensions without coordinates: vertices Data variables: rotated_latitude_longitude int32 -2147483647 lat_vertices (rlat, rlon, vertices) float32 22.02 ... 22.34 lon_vertices (rlat, rlon, vertices) float32 349.9 ... 350.1 orog (rlat, rlon) float32 265.9 258.9 ... 269.4 263.2 Attributes: (12/34) institution: Helmholtz-Zentrum Geesthacht, Climate Ser... institute_id: GERICS experiment_id: evaluation source: GERICS-REMO2015 model_id: GERICS-REMO2015 forcing: N/A ... ... title: GERICS-REMO2015 model output prepared for... parent_experiment: N/A modeling_realm: atmos realization: 0 cmor_version: 2.9.1 tracking_id: hdl:21.14103/c25ebf26-1647-49c8-b3fd-4e4e...
array([-23.375, -23.265, -23.155])
array([-28.375, -28.265, -28.155])
array([[21.987913, 22.027924, 22.0678 ], [22.088873, 22.128944, 22.168882], [22.189817, 22.22995 , 22.269953]], dtype=float32)
array([[349.9361 , 350.03607, 350.13617], [349.88898, 349.9891 , 350.08932], [349.8418 , 349.94205, 350.04242]], dtype=float32)
array(-2147483647, dtype=int32)
array([[[22.018372, 21.917442, 21.957424, 22.058414], [22.058414, 21.957424, 21.997337, 22.098389], [22.098389, 21.997337, 22.037117, 22.13823 ]], [[22.119293, 22.018372, 22.058414, 22.159397], [22.159397, 22.058414, 22.098389, 22.199432], [22.199432, 22.098389, 22.13823 , 22.239338]], [[22.220201, 22.119293, 22.159397, 22.260365], [22.260365, 22.159397, 22.199432, 22.300465], [22.300465, 22.199432, 22.239338, 22.340431]]], dtype=float32)
array([[[349.8625 , 349.90967, 350.0096 , 349.96255], [349.96255, 350.0096 , 350.10956, 350.06265], [350.06265, 350.10956, 350.20963, 350.16287]], [[349.81528, 349.8625 , 349.96255, 349.91547], [349.91547, 349.96255, 350.06265, 350.01572], [350.01572, 350.06265, 350.16287, 350.11606]], [[349.76797, 349.81528, 349.91547, 349.86832], [349.86832, 349.91547, 350.01572, 349.96872], [349.96872, 350.01572, 350.11606, 350.0692 ]]], dtype=float32)
array([[265.9244 , 258.8555 , 252.97688], [269.90173, 263.28323, 257.1503 ], [275.1792 , 269.3506 , 263.21967]], dtype=float32)
Let's plot lon_vertices
and lat_vertices
with a rotated pole projection and annotate cell and vertices indices.
rotated = ccrs.RotatedPole(*pole)
ax = plt.subplot(projection=rotated)
# plot cell centers
subset.orog.plot.pcolormesh(x="rlon", y="rlat", ax=ax)
# plot cell centers
# xr.plot.scatter(subset, x='lon', y='lat', transform=ccrs.PlateCarree())
xr.plot.scatter(subset, x="rlon", y="rlat", transform=rotated)
# plot cell vertices
xr.plot.scatter(
subset, x="lon_vertices", y="lat_vertices", transform=ccrs.PlateCarree()
)
# annotate cell centers with logical grid index
for (i, j), element in np.ndenumerate(subset.lon.values):
ax.annotate(
f"({i},{j})", (subset.rlon.isel(rlon=i), subset.rlat.isel(rlat=j))
) # xycoords=transform_plate)
transform_plate = ccrs.PlateCarree()._as_mpl_transform(ax)
i = 1
j = 1
for k in range(subset.vertices.size):
ax.annotate(
f" {k}",
(
subset.lon_vertices.isel(rlon=i, rlat=j, vertices=k),
subset.lat_vertices.isel(rlon=i, rlat=j, vertices=k),
),
xycoords=transform_plate,
color="red",
)
Note that the ordering of vertices starts at the top left (facing cell (i-1
, j+1
)). It seems to be different from what's expected in the CF-conventions which states that the 0-index should define the vertex facing cell (i-1
, j-1
). This causes problems in other CF compliant tools like, e.g., cf_xarray
, see also cf_xarray.bounds_to_vertices. This seems to be the same for a lot of cmorized datasets i looked at.
Just to demonstrate here, i use cf_xarray
to expand vertices into bounds of dimensions (rlon
+1, rlat
+1) and again plot those.
import cf_xarray as cfxr
def roll_vertices(ds):
"""roll vertices into an order beginning at lower left"""
ds = ds.copy()
ds["lon_vertices"] = ds.lon_vertices.roll(vertices=-1)
ds["lat_vertices"] = ds.lat_vertices.roll(vertices=-1)
return ds
def bounds_to_vertices(ds):
lon_bounds = cfxr.bounds_to_vertices(ds.lon_vertices, bounds_dim="vertices")
lat_bounds = cfxr.bounds_to_vertices(ds.lat_vertices, bounds_dim="vertices")
lon_bounds.name = "lon_bounds"
lat_bounds.name = "lat_bounds"
return xr.merge([lon_bounds, lat_bounds])
bounds = bounds_to_vertices(subset)
ax = plt.subplot(projection=rotated)
# plot cell centers
subset.orog.plot.pcolormesh(x="rlon", y="rlat", ax=ax)
xr.plot.scatter(subset.isel(rlon=1, rlat=1), x="rlon", y="rlat", transform=rotated)
# the bounds for cell (1,1) should be ((1,2), (1,2))
xr.plot.scatter(
bounds.isel(rlon_vertices=slice(1, 3), rlat_vertices=slice(1, 3)),
x="lon_bounds",
y="lat_bounds",
transform=ccrs.PlateCarree(),
)
# annotate cell centers with logical grid index
for (i, j), element in np.ndenumerate(subset.lon.values):
ax.annotate(f"({i},{j})", (subset.rlon.isel(rlon=i), subset.rlat.isel(rlat=j)))
So this shifts vertices to an incorrect position since cf_xarray
assumes the 0th vertices at the lower left. If we roll the vertices, we again get the correct positions of the vertices.
bounds = bounds_to_vertices(roll_vertices(subset))
ax = plt.subplot(projection=rotated)
# plot cell centers
subset.orog.plot.pcolormesh(x="rlon", y="rlat", ax=ax)
xr.plot.scatter(subset.isel(rlon=1, rlat=1), x="rlon", y="rlat", transform=rotated)
# the bounds for cell (1,1) should be ((1,2), (1,2))
xr.plot.scatter(
bounds.isel(rlon_vertices=slice(1, 3), rlat_vertices=slice(1, 3)),
x="lon_bounds",
y="lat_bounds",
transform=ccrs.PlateCarree(),
)
# annotate cell centers with logical grid index
for (i, j), element in np.ndenumerate(subset.lon.values):
ax.annotate(f"({i},{j})", (subset.rlon.isel(rlon=i), subset.rlat.isel(rlat=j)))
The problem is basically when it comes to conservative remapping. E.g., xesmf
relies on cf_xarray
to decode vertices correctly while cf_xarray
relies on the CF-convention. The question is how cdo
manages the vertices order internally, not only for the cdo cmor
operator but also for remapcon
if a curvilinear grid is present. At least, i experienced that cdo remapcon
is sensitive to the vertices order, e.g., i get different results with and without rolling the vertices before remapping.