ENH: Add support for coordinate systems with CRS using CF conventions (…

…#660)

docs/build_crs_cf.rst

@@ -0,0 +1,193 @@
+.. _build_crs_cf:
+
+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:
+
+.. code-block:: python
+
+    from pyproj import CRS
+
+    crs = CRS("EPSG:4326")
+
+
+To get the grid mapping you use :meth:`pyproj.crs.CRS.to_cf`:
+
+.. versionadded:: 2.2.0
+
+.. code-block:: python
+
+    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 :meth:`pyproj.crs.CRS.cs_to_cf`:
+
+.. versionadded:: 3.0.0
+
+.. code-block:: python
+
+    cf_coordinate_system = crs.cs_to_cf()
+
+
+Contents of `cf_coordinate_system`::
+
+    [{'long_name': 'geodetic latitude coordinate',
+    'standard_name': 'geodetic latitude',
+    'unit': 'degrees_north'},
+    {'long_name': 'geodetic longitude coordinate',
+    'standard_name': 'geodetic longitude',
+    'unit': 'degrees_east'}]
+
+
+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.
+
+.. code-block:: python
+
+    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:
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. versionadded:: 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
+
+
+.. code-block:: python
+
+    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,
+    )

docs/history.rst

@@ -15,6 +15,7 @@ Change Log
 * ENH: Use from_user_input in __eq__ when comparing CRS sub-classes (i.e. PrimeMeridian, Datum, Ellipsoid, etc.) (issue #606)
 * ENH: Added :meth:`pyproj.transformer.TransformerGroup.download_grids` (pull #643)
 * ENH: Added transformation grid sync API/CLI (issue #572)
+* ENH: Add support for coordinate systems with CRS using CF conventions (issue #536)
 
 2.6.1
 ~~~~~

docs/index.rst

@@ -18,12 +18,13 @@ GitHub Repository: https://github.com/pyproj4/pyproj
    examples
    transformation_grids
    gotchas
-   advanced_examples
-   build_crs
    api/index
    cli
-   optimize_transformations
+   advanced_examples
+   build_crs
+   build_crs_cf
    crs_compatibility
+   optimize_transformations
    history
    past_versions
 

pyproj/_crs.pyi

@@ -115,6 +115,7 @@ class CoordinateSystem(_CRSParts):
     def from_json_dict(coordinate_system_dict: dict) -> "CoordinateSystem": ...
     @staticmethod
     def from_json(coordinate_system_json_str: str) -> "CoordinateSystem": ...
+    def to_cf(self, rotated_pole: bool = False) -> List[dict]: ...
 
 class Param:
     name: str

pyproj/_crs.pyx

@@ -710,6 +710,85 @@ cdef class CoordinateSystem(_CRSParts):
             _load_proj_json(coordinate_system_json_str)
         )
 
+    def to_cf(self, rotated_pole=False):
+        """
+        .. versionadded:: 3.0.0
+
+        This converts a :obj:`pyproj.crs.CoordinateSystem` axis
+        to a list of Climate and Forecast (CF) Version 1.8 dicts.
+
+        Parameters
+        ----------
+        rotated_pole: bool, optional
+            If True, the geographic coordinates are on a rotated pole grid.
+            This corresponds to the rotated_latitude_longitude grid_mapping_name.
+            Default is False.
+
+        Returns
+        -------
+        List[dict]:
+            CF-1.8 version of the CoordinateSystem.
+        """
+        axis_list = self.to_json_dict()["axis"]
+        cf_params = []
+        def get_linear_unit(axis):
+            try:
+                return f'{axis["unit"]["conversion_factor"]} metre'
+            except TypeError:
+                return axis["unit"]
+
+        if self.name == "cartesian":
+            for axis in axis_list:
+                if axis["name"].lower() == "easting":
+                    cf_axis = "X"
+                else:
+                    cf_axis = "Y"
+                cf_params.append(dict(
+                    axis=cf_axis,
+                    long_name=axis["name"],
+                    standard_name=f"projection_{cf_axis.lower()}_coordinate",
+                    unit=get_linear_unit(axis),
+                ))
+        elif self.name == "ellipsoidal":
+            for axis in axis_list:
+                if axis["abbreviation"].upper() == "H":
+                    cf_params.append(dict(
+                        standard_name="height_above_reference_ellipsoid",
+                        long_name=axis["name"],
+                        unit=axis["unit"],
+                        positive=axis["direction"],
+                        axis="Z",
+                    ))
+                else:
+                    name = axis["name"].lower()
+                    if rotated_pole:
+                        cf_params.append(dict(
+                            standard_name=f"grid_{name}",
+                            long_name=f"{name} in rotated pole grid",
+                            unit="degrees",
+                        ))
+                    else:
+                        cf_params.append(dict(
+                            standard_name=name,
+                            long_name=f"{name} coordinate",
+                            unit=f'degrees_{axis["direction"]}',
+                        ))
+        elif self.name == "vertical":
+            for axis in axis_list:
+                if axis["abbreviation"].upper() == "H":
+                    standard_name = "height"
+                else:
+                    standard_name = "depth"
+                cf_params.append(dict(
+                    standard_name=standard_name,
+                    long_name=axis["name"],
+                    unit=get_linear_unit(axis),
+                    positive=axis["direction"],
+                    axis="Z",
+                ))
+
+        return cf_params
+
 
 cdef class Ellipsoid(_CRSParts):
     """