Managing CRS to and from CF

http://cfconventions.org/cf-conventions/cf-conventions.html

Exporting CRS to CF

When exporting a CRS to the Climate and Forecast (CF) conventions, you need both the grid mapping as well as the coordinate system. If you don’t use the coordinate system, then you will lose the units of your projection.

In this example, this is the CRS we will use:

from pyproj import CRS

crs = CRS("EPSG:4326")

To get the grid mapping you use pyproj.crs.CRS.to_cf():

Added in version 2.2.0.

cf_grid_mapping = crs.to_cf()

Contents of cf_grid_mapping:

{'crs_wkt': 'GEOGCRS["WGS 84",DATUM["World Geodetic System '
            ....,ID["EPSG",4326]]',
'geographic_crs_name': 'WGS 84',
'grid_mapping_name': 'latitude_longitude',
'inverse_flattening': 298.257223563,
'longitude_of_prime_meridian': 0.0,
'prime_meridian_name': 'Greenwich',
'reference_ellipsoid_name': 'WGS 84',
'semi_major_axis': 6378137.0,
'semi_minor_axis': 6356752.314245179}

To get the coordinate system, you use pyproj.crs.CRS.cs_to_cf():

Added in version 3.0.0.

cf_coordinate_system = crs.cs_to_cf()

Contents of cf_coordinate_system:

[{'long_name': 'geodetic latitude coordinate',
'standard_name': 'latitude',
'units': 'degrees_north',
'axis': 'Y'},
{'long_name': 'geodetic longitude coordinate',
'standard_name': 'longitude',
'units': 'degrees_east',
'axis': 'X'}]

Importing CRS from CF

When importing a CRS from the Climate and Forecast (CF) conventions, you need both the grid mapping as well as the coordinate system. If you don’t use the coordinate system, then you will lose the units of your projection.

Note

If the CF crs_wkt attribute is available, the coordinate system is inside of the WKT and can be used to create the CRS in a single step.

Warning

If building from grid mapping, be mindful of the axis order. https://github.com/cf-convention/cf-conventions/pull/224

Build the CRS from CF grid mapping:

In this example, this is the grid mapping and coordinate system we will use:

variables:
  double x(x) ;
    x:standard_name = "projection_x_coordinate" ;
    x:long_name = "Easting" ;
    x:units = "m" ;
  double y(y) ;
    y:standard_name = "projection_y_coordinate" ;
    y:long_name = "Northing" ;
    y:units = "m" ;
  int crsOSGB ;
    crsOSGB:grid_mapping_name = "transverse_mercator";
    crsOSGB:semi_major_axis = 6377563.396 ;
    crsOSGB:inverse_flattening = 299.3249646 ;
    crsOSGB:longitude_of_prime_meridian = 0.0 ;
    crsOSGB:latitude_of_projection_origin = 49.0 ;
    crsOSGB:longitude_of_central_meridian = -2.0 ;
    crsOSGB:scale_factor_at_central_meridian = 0.9996012717 ;
    crsOSGB:false_easting = 400000.0 ;
    crsOSGB:false_northing = -100000.0 ;

Note

If the units are meters as in this example, then no further changes are necessary.

from pyproj import CRS

crs = CRS.from_cf(
    {
        "grid_mapping_name": "transverse_mercator",
        "semi_major_axis": 6377563.396,
        "inverse_flattening": 299.3249646,
        "longitude_of_prime_meridian": 0.0,
        "latitude_of_projection_origin": 49.0,
        "longitude_of_central_meridian": -2.0,
        "scale_factor_at_central_meridian": 0.9996012717,
        "false_easting": 400000.0,
        "false_northing": -100000.0,
    }
)

Modify the CRS with coordinate system:

Added in version 3.0.0.

Note

If the CF crs_wkt attribute is available, the coordinate system is inside of the WKT and can be used to create the CRS in a single step.

Warning

Be mindful of the axis order. https://github.com/cf-convention/cf-conventions/pull/224

In this example, assume everything is the same as above. However, the units are instead US_Survey_Foot:

variables:
  double x(x) ;
    x:standard_name = "projection_x_coordinate" ;
    x:long_name = "Easting" ;
    x:units = "US_Survey_Foot" ;
  double y(y) ;
    y:standard_name = "projection_y_coordinate" ;
    y:long_name = "Northing" ;
    y:units = "US_Survey_Foot" ;
  ...

In this case, you will need to get the unit conversion factor:

https://github.com/SciTools/cf-units

from cf_units import Unit
from pyproj import CRS

cf_unit = Unit("US_Survey_Foot")
unit = {
    "type": "LinearUnit",
    "name": "US Survey Foot",
    "conversion_factor": cf_unit.convert(1, "m"),
}
cartesian_cs = {
    "$schema": "https://proj.org/schemas/v0.2/projjson.schema.json",
    "type": "CoordinateSystem",
    "subtype": "Cartesian",
    "axis": [
        {"name": "Easting", "abbreviation": "E", "direction": "east", "unit": unit},
        {"name": "Northing", "abbreviation": "N", "direction": "north", "unit": unit},
    ],
}
crs = CRS.from_cf(
    {
        "grid_mapping_name": "transverse_mercator",
        "semi_major_axis": 6377563.396,
        "inverse_flattening": 299.3249646,
        "longitude_of_prime_meridian": 0.0,
        "latitude_of_projection_origin": 49.0,
        "longitude_of_central_meridian": -2.0,
        "scale_factor_at_central_meridian": 0.9996012717,
        "false_easting": 400000.0,
        "false_northing": -100000.0,
    },
    cartesian_cs=cartesian_cs,
)