Projection lat-lon coordinates including new GridMapping objects

cf/__init__.py

@@ -254,6 +254,9 @@
 from .domainaxis import DomainAxis
 from .fieldancillary import FieldAncillary
 from .field import Field
+
+from .gridmappings import *  # noqa: F403
+
 from .data import Data
 from .data.array import (
     CFANetCDFArray,

cf/constants.py

@@ -72,26 +72,123 @@
 
 _stash2standard_name = {}
 
+
 # ---------------------------------------------------------------------
-# Coordinate reference constants TODO: turn these into functions
+# Coordinate reference and grid mapping related constants
 # ---------------------------------------------------------------------
+# Grid Mapping valid attribute names. Taken from the list given in
+# 'Table F.1. Grid Mapping Attributes' in Appendix F: Grid Mappings'
+# of the Conventions document.
+#
+# Values indicate if attribute values are expected to be numeric (instead
+# of string) for the given attribute key (the table defines this via N, S).
+# For any value where this is True, canonical units must be provided
+# in the 'cr_canonical_units' dict below.
+cr_gm_valid_attr_names_are_numeric = {
+    "grid_mapping_name": False,
+    # *Those which describe the ellipsoid and prime meridian:*
+    "earth_radius": True,
+    "inverse_flattening": True,
+    "longitude_of_prime_meridian": True,
+    "prime_meridian_name": False,
+    "reference_ellipsoid_name": False,
+    "semi_major_axis": True,
+    "semi_minor_axis": True,
+    # *Specific/applicable to only given grid mapping(s):*
+    # ...projection origin related:
+    "longitude_of_projection_origin": True,
+    "latitude_of_projection_origin": True,
+    "scale_factor_at_projection_origin": True,
+    # ...false-Xings:
+    "false_easting": True,
+    "false_northing": True,
+    # ...angle axis related:
+    "sweep_angle_axis": False,
+    "fixed_angle_axis": False,
+    # ...central meridian related:
+    "longitude_of_central_meridian": True,
+    "scale_factor_at_central_meridian": True,
+    # ...pole coordinates related:
+    "grid_north_pole_latitude": True,
+    "grid_north_pole_longitude": True,
+    "north_pole_grid_longitude": True,
+    # ...other:
+    "standard_parallel": True,
+    "perspective_point_height": True,
+    "azimuth_of_central_line": True,
+    "straight_vertical_longitude_from_pole": True,
+    # *Other, not needed for a specific grid mapping but also listed
+    # in 'Table F.1. Grid Mapping Attributes':*
+    "crs_wkt": False,
+    "geographic_crs_name": False,
+    "geoid_name": False,
+    "geopotential_datum_name": False,
+    "horizontal_datum_name": False,
+    "projected_crs_name": False,
+    "towgs84": True,
+}
 cr_canonical_units = {
+    # *Those which describe the ellipsoid and prime meridian:*
     "earth_radius": Units("m"),
-    "grid_north_pole_latitude": Units("degrees_north"),
-    "grid_north_pole_longitude": Units("degrees_east"),
     "inverse_flattening": Units("1"),
-    "latitude_of_projection_origin": Units("degrees_north"),
-    "longitude_of_central_meridian": Units("degrees_east"),
     "longitude_of_prime_meridian": Units("degrees_east"),
-    "longitude_of_projection_origin": Units("degrees_east"),
-    "north_pole_grid_longitude": Units("degrees"),
-    "perspective_point_height": Units("m"),
-    "scale_factor_at_central_meridian": Units("1"),
-    "scale_factor_at_projection_origin": Units("1"),
     "semi_major_axis": Units("m"),
     "semi_minor_axis": Units("m"),
+    # *Specific/applicable to only given grid mapping(s):*
+    # ...projection origin related:
+    "longitude_of_projection_origin": Units("degrees_east"),
+    "latitude_of_projection_origin": Units("degrees_north"),
+    "scale_factor_at_projection_origin": Units("1"),
+    # ...false-Xings:
+    "false_easting": Units("m"),
+    "false_northing": Units("m"),
+    # ...central meridian related:
+    "longitude_of_central_meridian": Units("degrees_east"),
+    "scale_factor_at_central_meridian": Units("1"),
+    # ...pole coordinates related:
+    "grid_north_pole_latitude": Units("degrees_north"),
+    "grid_north_pole_longitude": Units("degrees_east"),
+    "north_pole_grid_longitude": Units("degrees"),
+    # ...other:
     "standard_parallel": Units("degrees_north"),
-    "straight_vertical_longitude_from_pole": Units("degrees_north"),
+    "perspective_point_height": Units("m"),
+    "azimuth_of_central_line": Units("degrees"),
+    # TODOGM check the below is correct, was 'degrees_north' before but it
+    # is a longitude so surely that was wrong?:
+    "straight_vertical_longitude_from_pole": Units("degrees_east"),
+    # *Other, not needed for a specific grid mapping but also listed
+    # in 'Table F.1. Grid Mapping Attributes':*
+    # TODOGM towgs84 is very complicated, with multiple components
+    # of translations, roations and scaling each with different units
+    # requirements, so we should probably sort this out later, see:
+    # https://proj.org/en/9.3/operations/transformations/
+    # helmert.html#helmert-transform
+    # where translations are in meters, rotations in arc-seconds, and
+    # scaling is in parts per million.
+    "towgs84": None,
+}
+
+# Indicates what PROJ string ID component(s) refer(s) to the grid mapping
+# attribute in question when '+' is added as a suffix, for example
+# 'earth_radius' is '+R' and 'inverse_flattening' is '+rf'
+cr_gm_attr_to_proj_string_comps = {
+    "earth_radius": "R",
+    "inverse_flattening": "rf",
+    "prime_meridian_name": "pm",
+    "reference_ellipsoid_name": "ellps",
+    "semi_major_axis": "a",
+    "semi_minor_axis": "b",
+    "longitude_of_projection_origin": "lon_0",
+    "latitude_of_projection_origin": "lat_0",
+    "scale_factor_at_projection_origin": "k_0",
+    "false_easting": "x_0",
+    "false_northing": "y_0",
+    "sweep_angle_axis": "sweep",
+    "standard_parallel": ("lat_1", "lat_2"),
+    "perspective_point_height": "h",
+    "azimuth_of_central_line": ("alpha", "gamma"),
+    "crs_wkt": "crs_wkt",
+    "towgs84": "towgs84",
 }
 
 cr_coordinates = {

cf/domain.py

@@ -18,6 +18,7 @@
     indices_shape,
     parse_indices,
 )
+from .gridmappings import GM, InvalidGridMapping
 
 _empty_set = set()
 
@@ -585,6 +586,59 @@ def get_data(self, default=ValueError(), _units=None, _fill_value=True):
             default, message=f"{self.__class__.__name__} has no data"
         )
 
+    def get_grid_mappings(self, as_class=False):
+        """Returns coordinate conversions with their grid mappings.
+
+        .. versionadded:: GMVER
+
+        :Parameters:
+
+            as_class: `bool`, optional
+                If `True`, return the grid mapping as the equivalent
+                CF Grid Mapping class, for example
+                cf.RotatedLatitudeLongitude, rather than as a string
+                corresponding to the value of the 'grid_mapping_name'
+                attribute, for example 'rotated_latitude_longitude'.
+                By default the 'grid_mapping_name' value is returned.
+
+        :Returns:
+
+                `dict`
+                     CoordinateConversion construct identifiers with
+                     values of their 'grid_mapping_name' attribute,
+                     or corresponding CF Grid Mapping class if
+                     as_class is `True`, for all CoordinateConversions
+                     of the domain that have a 'grid_mapping_name'
+                     parameter defined.
+
+        **Examples**
+
+        >>> f.get_grid_mappings()
+        {'coordinatereference1': "rotated_latitude_longitude"}
+        >>> f.get_grid_mappings(as_class=True)
+        {'coordinatereference1': cf.gridmappings.RotatedLatitudeLongitude}
+        >>> g.get_grid_mappings()
+        {}
+
+        """
+        gms = {}
+        for cref_name, cref in self.coordinate_references().items():
+            # If there is no coordinate conversion or parameters set on one,
+            # this will give None, so is safe
+            gm = cref.coordinate_conversion.get_parameter(
+                "grid_mapping_name", default=None
+            )
+            if gm:
+                if as_class:
+                    try:
+                        gm_class = GM(cref)
+                        gms[cref_name] = gm_class
+                    except InvalidGridMapping:
+                        pass  # not a supported GM so don't add
+                else:
+                    gms[cref_name] = gm
+        return gms
+
     def identity(self, default="", strict=False, relaxed=False, nc_only=False):
         """Return the canonical identity.
 

cf/field.py

@@ -15667,3 +15667,40 @@ def unlimited(self, *args):
             version="3.0.0",
             removed_at="4.0.0",
         )  # pragma: no cover
+
+    def get_grid_mappings(self, as_class=False):
+        """Returns coordinate conversions with their grid mappings.
+
+        .. versionadded:: GMVER
+
+        :Parameters:
+
+            as_class: `bool`, optional
+                If `True`, return the grid mapping as the equivalent
+                CF Grid Mapping class, for example
+                cf.RotatedLatitudeLongitude, rather than as a string
+                corresponding to the value of the 'grid_mapping_name'
+                attribute, for example 'rotated_latitude_longitude'.
+                By default the 'grid_mapping_name' value is returned.
+
+        :Returns:
+
+                `dict`
+                     CoordinateConversion construct identifiers with
+                     values of their 'grid_mapping_name' attribute,
+                     or corresponding CF Grid Mapping class if
+                     as_class is `True`, for all CoordinateConversions
+                     of the domain that have a 'grid_mapping_name'
+                     parameter defined.
+
+        **Examples**
+
+        >>> f.get_grid_mappings()
+        {'coordinatereference1': "rotated_latitude_longitude"}
+        >>> f.get_grid_mappings(as_class=True)
+        {'coordinatereference1': cf.gridmappings.RotatedLatitudeLongitude}
+        >>> g.get_grid_mappings()
+        {}
+
+        """
+        return self.domain.get_grid_mappings(as_class=as_class)

cf/gridmappings/__init__.py

@@ -0,0 +1,40 @@
+"""
+Module for Grid Mappings supported by the CF Conventions.
+
+For the full list of supported Grid Mappings and details, see Appendix F
+of the canonical document:
+
+https://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/
+cf-conventions.html#appendix-grid-mappings
+
+This module should be kept up to date with the Appendix, by adding or
+amending appropriate classes.
+
+Note that abstract classes to support the creation of concrete classes
+for concrete Grid Mappinds are defined in the 'abstract' module, so
+not included in the listing below.
+
+"""
+
+
+from .abstract import *
+from .gridmapping import (
+    GM,
+    InvalidGridMapping,
+    AlbersEqualArea,
+    AzimuthalEquidistant,
+    Geostationary,
+    LambertAzimuthalEqualArea,
+    LambertConformalConic,
+    LambertCylindricalEqualArea,
+    Mercator,
+    ObliqueMercator,
+    Orthographic,
+    PolarStereographic,
+    RotatedLatitudeLongitude,
+    LatitudeLongitude,
+    Sinusoidal,
+    Stereographic,
+    TransverseMercator,
+    VerticalPerspective,
+)

cf/gridmappings/abstract/__init__.py

@@ -0,0 +1,25 @@
+"""
+Abstract classes representing categories of Grid Mapping.
+
+Categories correspond to the type of projection a Grid Mapping is
+based upon, for example these are often based upon a geometric shape
+that forms the developable surface that is used to flatten the map
+of Earth (such as conic for a cone or cylindrical for a cylinder).
+
+The LatLonGridMapping case is special in that it (from the CF Conventions,
+Appendix F) 'defines the canonical 2D geographical coordinate system
+based upon latitude and longitude coordinates on a spherical Earth'.
+
+"""
+
+from .gridmappingbase import (
+    GridMapping,
+    convert_proj_angular_data_to_cf,
+    convert_cf_angular_data_to_proj,
+    _validate_map_parameter,
+)
+from .azimuthal import AzimuthalGridMapping
+from .conic import ConicGridMapping
+from .cylindrical import CylindricalGridMapping
+from .latlon import LatLonGridMapping
+from .perspective import PerspectiveGridMapping

cf/gridmappings/abstract/azimuthal.py

@@ -0,0 +1,65 @@
+from .gridmappingbase import (
+    GridMapping,
+    _validate_map_parameter,
+)
+
+
+class AzimuthalGridMapping(GridMapping):
+    """A Grid Mapping with Azimuthal classification.
+
+    .. versionadded:: GMVER
+
+    :Parameters:
+
+        longitude_of_projection_origin: number or scalar `Data`, optional
+            The longitude of projection center (PROJ 'lon_0' value).
+            If provided as a number or `Data` without units, the units
+            are taken as 'degrees_east', else the `Data` units are
+            taken and must be angular and compatible with longitude.
+            The default is 0.0 degrees_east.
+
+        latitude_of_projection_origin: number or scalar `Data`, optional
+            The latitude of projection center (PROJ 'lat_0' value).
+            If provided as a number or `Data` without units, the units
+            are taken as 'degrees_north', else the `Data` units are
+            taken and must be angular and compatible with latitude.
+            The default is 0.0 degrees_north.
+
+        false_easting: number or scalar `Data`, optional
+            The false easting (PROJ 'x_0') value.
+            If provided as a number or `Data` without units, the units
+            are taken as metres 'm', else the `Data` units are
+            taken and must be compatible with distance.
+            The default is 0.0 metres.
+
+        false_northing: number or scalar `Data`, optional
+            The false northing (PROJ 'y_0') value.
+            If provided as a number or `Data` without units, the units
+            are taken as metres 'm', else the `Data` units are
+            taken and must be compatible with distance.
+            The default is 0.0 metres.
+
+    """
+
+    def __init__(
+        self,
+        longitude_of_projection_origin=0.0,
+        latitude_of_projection_origin=0.0,
+        false_easting=0.0,
+        false_northing=0.0,
+        **kwargs,
+    ):
+        super().__init__(**kwargs)
+
+        self.longitude_of_projection_origin = _validate_map_parameter(
+            "longitude_of_projection_origin", longitude_of_projection_origin
+        )
+        self.latitude_of_projection_origin = _validate_map_parameter(
+            "latitude_of_projection_origin", latitude_of_projection_origin
+        )
+        self.false_easting = _validate_map_parameter(
+            "false_easting", false_easting
+        )
+        self.false_northing = _validate_map_parameter(
+            "false_northing", false_northing
+        )