Grid Mappings#

See

  1. Dataset.cf.grid_mapping_names,

  2. DataArray.cf.grid_mapping_name

cf_xarray understands the concept of coordinate projections using the grid_mapping attribute convention. For example, the dataset might contain two sets of coordinates:

  • native coordinates in which the data is defined, e.g., regular 1D coordinates

  • projected coordinates which probably denote some “real” coordinates in latitude and longitude

Due to the projection, those real coordinates are probably 2D data variables. The grid_mapping attribute of a data variable makes a connection to another data variable defining the coordinate reference system (CRS) of those native coordinates. It should enable you to project the native coordinates into any other CRS, including the real 2D latitude and longitude coordinates. This is often useful for plotting, e.g., you can tell cartopy how to correctly plot coastlines for the CRS your data is defined in.

Extracting grid mapping info#

Dataset#

To access grid_mapping attributes, consider this example:

from cf_xarray.datasets import rotds

rotds

<xarray.Dataset> Size: 844B
Dimensions:       (rlon: 3, rlat: 3, bounds: 4)
Coordinates:
  * rlon          (rlon) float64 24B 17.93 18.05 18.16
  * rlat          (rlat) float64 24B 21.61 21.73 21.84
    lon           (rlon, rlat) float64 72B 64.22 64.42 64.63 ... 64.76 64.96
    lat           (rlon, rlat) float64 72B 66.64 66.58 66.52 ... 66.75 66.69
Dimensions without coordinates: bounds
Data variables:
    lon_bounds    (bounds, rlon, rlat) float64 288B 64.03 64.24 ... 64.74 64.95
    lat_bounds    (bounds, rlon, rlat) float64 288B 66.63 66.56 ... 66.83 66.76
    rotated_pole  int32 4B 0
    temp          (rlat, rlon) float64 72B 0.6337 0.6367 ... 0.06773 0.2431

The related grid mappings can be discovered using Dataset.cf.grid_mapping_names which maps a “grid mapping name” to the appropriate variable name:

rotds.cf.grid_mapping_names

{'rotated_latitude_longitude': ['rotated_pole']}

Access the grid_mapping variable as

rotds.cf["grid_mapping"]

<xarray.DataArray 'rotated_pole' ()> Size: 4B
array(0, dtype=int32)
Attributes:
    grid_mapping_name:          rotated_latitude_longitude
    grid_north_pole_latitude:   39.25
    grid_north_pole_longitude:  -162.0

DataArrays#

Grid mapping variables are propagated when extracting DataArrays:

da = rotds.cf["temp"]
da

<xarray.DataArray 'temp' (rlat: 3, rlon: 3)> Size: 72B
array([[0.63372339, 0.63671979, 0.51343171],
       [0.35921895, 0.17884925, 0.30960602],
       [0.43538307, 0.06773382, 0.24308926]])
Coordinates:
  * rlon          (rlon) float64 24B 17.93 18.05 18.16
  * rlat          (rlat) float64 24B 21.61 21.73 21.84
    rotated_pole  int32 4B 0
Attributes:
    standard_name:  air_temperature
    grid_mapping:   rotated_pole

To find the grid mapping name use the singular DataArray.cf.grid_mapping_name

da.cf.grid_mapping_name

'rotated_latitude_longitude'

And to get the grid mapping variable

da.cf["grid_mapping"]

<xarray.DataArray 'rotated_pole' ()> Size: 4B
array(0, dtype=int32)
Attributes:
    grid_mapping_name:          rotated_latitude_longitude
    grid_north_pole_latitude:   39.25
    grid_north_pole_longitude:  -162.0

Use grid_mapping in projections#

The grid mapping information use very useful in projections, e.g., for plotting. pyproj understands CF conventions right away, e.g.

from pyproj import CRS

CRS.from_cf(rotds.cf["grid_mapping"].attrs)

gives you more details on the projection:

<Derived Geographic 2D CRS: {"$schema": "https://proj.org/schemas/v0.2/projjso ...>
Name: undefined
Axis Info [ellipsoidal]:
- lon[east]: Longitude (degree)
- lat[north]: Latitude (degree)
Area of Use:
- undefined
Coordinate Operation:
- name: Pole rotation (netCDF CF convention)
- method: Pole rotation (netCDF CF convention)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

For use in cartopy, there is some more overhead due to this issue. So you should select the right cartopy CRS and just feed in the grid mapping info:

from cartopy import crs as ccrs

grid_mapping = rotds.cf["grid_mapping"]
pole_latitude = grid_mapping.grid_north_pole_latitude
pole_longitude = grid_mapping.grid_north_pole_longitude

ccrs.RotatedPole(pole_longitude, pole_latitude)

cartopy rotated pole projection