Datum support in Iris

docs/src/whatsnew/latest.rst

@@ -37,7 +37,7 @@ This document explains the changes made to Iris for this release
    :pull:`4589`)
 
 #. `@wjbenfold`_ added support for ``false_easting`` and ``false_northing`` to
-   :class:`~iris.coord_system.Mercator`. (:issue:`3107`, :pull:`4524`)
+   :class:`~iris.coord_systems.Mercator`. (:issue:`3107`, :pull:`4524`)
 
 #. `@rcomer`_ and `@wjbenfold`_ (reviewer) implemented lazy aggregation for the
    :obj:`iris.analysis.PERCENTILE` aggregator. (:pull:`3901`)
@@ -47,7 +47,13 @@ This document explains the changes made to Iris for this release
 
 #. `@wjbenfold`_ added support for CF-compliant treatment of
    ``standard_parallel`` and ``scale_factor_at_projection_origin`` to
-   :class:`~iris.coord_system.Mercator`. (:issue:`3844`, :pull:`4609`)
+   :class:`~iris.coord_systems.Mercator`. (:issue:`3844`, :pull:`4609`)
+
+#. `@wjbenfold`_ added support datums associated with coordinate systems (e.g.
+   :class:`~iris.coord_systems.GeogCS` other subclasses of
+   :class:`~iris.coord_systems.CoordSystem`). Loading of datum information from
+   a netCDF file only happens when the :obj:`iris.FUTURE.datum_support` flag is
+   set. (:issue:`4619`, :pull:`4704`)
 
 
 🐛 Bugs Fixed

lib/iris/__init__.py

@@ -136,7 +136,7 @@ def callback(cube, field, filename):
 class Future(threading.local):
     """Run-time configuration controller."""
 
-    def __init__(self):
+    def __init__(self, datum_support=False):
         """
         A container for run-time options controls.
 
@@ -151,22 +151,24 @@ def __init__(self):
         .. note::
 
             iris.FUTURE.example_future_flag does not exist. It is provided
-            as an example because there are currently no flags in
-            iris.Future.
+            as an example.
 
         """
-        # The flag 'example_future_flag' is provided as a future reference
-        # for the structure of this class.
+        # The flag 'example_future_flag' is provided as a reference for the
+        # structure of this class.
+        #
+        # Note that self.__dict__ is used explicitly due to the manner in which
+        # __setattr__ is overridden.
         #
         # self.__dict__['example_future_flag'] = example_future_flag
-        pass
+        self.__dict__["datum_support"] = datum_support
 
     def __repr__(self):
 
         # msg = ('Future(example_future_flag={})')
         # return msg.format(self.example_future_flag)
-        msg = "Future()"
-        return msg.format()
+        msg = "Future(datum_support={})"
+        return msg.format(self.datum_support)
 
     # deprecated_options = {'example_future_flag': 'warning',}
     deprecated_options = {}
@@ -211,8 +213,7 @@ def context(self, **kwargs):
         .. note::
 
             iris.FUTURE.example_future_flag does not exist and is
-            provided only as an example since there are currently no
-            flags in Future.
+            provided only as an example.
 
         """
         # Save the current context

lib/iris/coord_systems.py

@@ -122,6 +122,12 @@ def as_cartopy_projection(self):
         pass
 
 
+_short_datum_names = {
+    "OSGB 1936": "OSGB36",
+    "WGS 84": "WGS84",
+}
+
+
 class GeogCS(CoordSystem):
     """
     A geographic (ellipsoidal) coordinate system, defined by the shape of
@@ -139,29 +145,28 @@ def __init__(
         longitude_of_prime_meridian=None,
     ):
         """
-        Creates a new GeogCS.
-
-        Kwargs:
+        Create a new GeogCS.
 
+        Parameters
+        ----------
         * semi_major_axis, semi_minor_axis:
-            Axes of ellipsoid, in metres.  At least one must be given
-            (see note below).
-
+            Axes of ellipsoid, in metres.  At least one must be given (see note
+            below).
         * inverse_flattening:
-            Can be omitted if both axes given (see note below).
-            Defaults to 0.0 .
-
+            Can be omitted if both axes given (see note below). Default 0.0
         * longitude_of_prime_meridian:
-            Specifies the prime meridian on the ellipsoid, in degrees.
-            Defaults to 0.0 .
+            Specifies the prime meridian on the ellipsoid, in degrees. Default 0.0
+        * datum:
+            Name of the datum used to modify the ellipsoid. Default None
 
+        Notes
+        -----
         If just semi_major_axis is set, with no semi_minor_axis or
         inverse_flattening, then a perfect sphere is created from the given
         radius.
 
-        If just two of semi_major_axis, semi_minor_axis, and
-        inverse_flattening are given the missing element is calculated from the
-        formula:
+        If just two of semi_major_axis, semi_minor_axis, and inverse_flattening
+        are given the missing element is calculated from the formula:
         :math:`flattening = (major - minor) / major`
 
         Currently, Iris will not allow over-specification (all three ellipsoid
@@ -194,58 +199,34 @@ def __init__(
         ):
             raise ValueError("Ellipsoid is overspecified")
 
-        # Perfect sphere (semi_major_axis only)? (1 0 0)
-        elif semi_major_axis is not None and (
-            semi_minor_axis is None and not inverse_flattening
+        # We didn't get enough to specify an ellipse.
+        if semi_major_axis is None and (
+            semi_minor_axis is None or inverse_flattening is None
         ):
-            semi_minor_axis = semi_major_axis
-            inverse_flattening = 0.0
+            raise ValueError("Insufficient ellipsoid specification")
 
-        # Calculate semi_major_axis? (0 1 1)
-        elif semi_major_axis is None and (
-            semi_minor_axis is not None and inverse_flattening is not None
-        ):
+        # Making a globe needs a semi_major_axis and a semi_minor_axis
+        if semi_major_axis is None:
             semi_major_axis = -semi_minor_axis / (
                 (1.0 - inverse_flattening) / inverse_flattening
             )
-
-        # Calculate semi_minor_axis? (1 0 1)
-        elif semi_minor_axis is None and (
-            semi_major_axis is not None and inverse_flattening is not None
-        ):
+        if semi_minor_axis is None and inverse_flattening:
             semi_minor_axis = semi_major_axis - (
                 (1.0 / inverse_flattening) * semi_major_axis
             )
 
-        # Calculate inverse_flattening? (1 1 0)
-        elif inverse_flattening is None and (
-            semi_major_axis is not None and semi_minor_axis is not None
-        ):
-            if semi_major_axis == semi_minor_axis:
-                inverse_flattening = 0.0
-            else:
-                inverse_flattening = 1.0 / (
-                    (semi_major_axis - semi_minor_axis) / semi_major_axis
-                )
-
-        # We didn't get enough to specify an ellipse.
-        else:
-            raise ValueError("Insufficient ellipsoid specification")
-
-        #: Major radius of the ellipsoid in metres.
-        self.semi_major_axis = float(semi_major_axis)
-
-        #: Minor radius of the ellipsoid in metres.
-        self.semi_minor_axis = float(semi_minor_axis)
-
-        #: :math:`1/f` where :math:`f = (a-b)/a`.
-        self.inverse_flattening = float(inverse_flattening)
-
         #: Describes 'zero' on the ellipsoid in degrees.
         self.longitude_of_prime_meridian = _arg_default(
             longitude_of_prime_meridian, 0
         )
 
+        globe = ccrs.Globe(
+            ellipse=None,
+            semimajor_axis=semi_major_axis,
+            semiminor_axis=semi_minor_axis,
+        )
+        self._crs = ccrs.Geodetic(globe)
+
     def _pretty_attrs(self):
         attrs = [("semi_major_axis", self.semi_major_axis)]
         if self.semi_major_axis != self.semi_minor_axis:
@@ -257,6 +238,14 @@ def _pretty_attrs(self):
                     self.longitude_of_prime_meridian,
                 )
             )
+        # An unknown crs datum will be treated as None
+        if self.datum is not None and self.datum != "unknown":
+            attrs.append(
+                (
+                    "datum",
+                    self.datum,
+                )
+            )
         return attrs
 
     def __repr__(self):
@@ -294,7 +283,7 @@ def xml_element(self, doc):
         return CoordSystem.xml_element(self, doc, attrs)
 
     def as_cartopy_crs(self):
-        return ccrs.Geodetic(self.as_cartopy_globe())
+        return self._crs
 
     def as_cartopy_projection(self):
         return ccrs.PlateCarree(
@@ -303,13 +292,38 @@ def as_cartopy_projection(self):
         )
 
     def as_cartopy_globe(self):
-        # Explicitly set `ellipse` to None as a workaround for
-        # Cartopy setting WGS84 as the default.
-        return ccrs.Globe(
-            semimajor_axis=self.semi_major_axis,
-            semiminor_axis=self.semi_minor_axis,
-            ellipse=None,
+        return self._crs.globe
+
+    def __getattr__(self, name):
+        if name == "semi_major_axis":
+            return self._crs.ellipsoid.semi_major_metre
+        if name == "semi_minor_axis":
+            return self._crs.ellipsoid.semi_minor_metre
+        if name == "inverse_flattening":
+            return self._crs.ellipsoid.inverse_flattening
+        if name == "datum":
+            datum = self._crs.datum.name
+            # An unknown crs datum will be treated as None
+            if datum == "unknown":
+                return None
+            return datum
+        return getattr(super(), name)
+
+    @classmethod
+    def from_datum(cls, datum, longitude_of_prime_meridian=None):
+
+        short_datum = _short_datum_names.get(datum, datum)
+
+        # Cartopy doesn't actually enact datums unless they're provided without
+        # ellipsoid axes, so only provide the datum
+        crs = super().__new__(cls)
+        crs._crs = ccrs.Geodetic(ccrs.Globe(short_datum, ellipse=None))
+
+        #: Describes 'zero' on the ellipsoid in degrees.
+        crs.longitude_of_prime_meridian = _arg_default(
+            longitude_of_prime_meridian, 0
         )
+        return crs
 
 
 class RotatedGeogCS(CoordSystem):
@@ -1110,6 +1124,10 @@ def __init__(
         * false_northing:
             Y offset from the planar origin in metres. Defaults to 0.0.
 
+        * datum:
+            If given, specifies the datumof the coordinate system. Only
+            respected if iris.Future.daum_support is set.
+
         Note: Only one of ``standard_parallel`` and
         ``scale_factor_at_projection_origin`` should be included.