Datum performance improvements

lib/iris/coord_systems.py

@@ -9,6 +9,7 @@
 """
 
 from abc import ABCMeta, abstractmethod
+from functools import cached_property
 import warnings
 
 import cartopy.crs as ccrs
@@ -53,10 +54,28 @@ class CoordSystem(metaclass=ABCMeta):
     grid_mapping_name = None
 
     def __eq__(self, other):
-        return (
-            self.__class__ == other.__class__
-            and self.__dict__ == other.__dict__
-        )
+        """
+        Override equality
+
+        The `_globe` and `_crs` attributes are not compared because they are
+        cached properties and completely derived from other attributes. The
+        nature of caching means that they can appear on one object and not on
+        another despite the objects being identical, and them being completely
+        derived from other attributes means they will only differ if other
+        attributes that are being tested for equality differ.
+        """
+        if self.__class__ != other.__class__:
+            return False
+        self_keys = set(self.__dict__.keys())
+        other_keys = set(other.__dict__.keys())
+        check_keys = (self_keys | other_keys) - {"_globe", "_crs"}
+        for key in check_keys:
+            try:
+                if self.__dict__[key] != other.__dict__[key]:
+                    return False
+            except KeyError:
+                return False
+        return True
 
     def __ne__(self, other):
         # Must supply __ne__, Python does not defer to __eq__ for
@@ -132,7 +151,6 @@ class GeogCS(CoordSystem):
     """
     A geographic (ellipsoidal) coordinate system, defined by the shape of
     the Earth and a prime meridian.
-
     """
 
     grid_mapping_name = "latitude_longitude"
@@ -156,8 +174,6 @@ def __init__(
             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. Default 0.0
-        * datum:
-            Name of the datum used to modify the ellipsoid. Default None
 
         Notes
         -----
@@ -172,6 +188,10 @@ def __init__(
         Currently, Iris will not allow over-specification (all three ellipsoid
         parameters).
 
+        After object creation, altering any of these properties will not update
+        the others. semi_major_axis and semi_minor_axis are used when creating
+        Cartopy objects.
+
         Examples::
 
             cs = GeogCS(6371229)
@@ -199,34 +219,60 @@ def __init__(
         ):
             raise ValueError("Ellipsoid is overspecified")
 
-        # 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
+        # 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
         ):
-            raise ValueError("Insufficient ellipsoid specification")
+            semi_minor_axis = semi_major_axis
+            inverse_flattening = 0.0
 
-        # Making a globe needs a semi_major_axis and a semi_minor_axis
-        if semi_major_axis is None:
+        # 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
+        ):
             semi_major_axis = -semi_minor_axis / (
                 (1.0 - inverse_flattening) / inverse_flattening
             )
-        if semi_minor_axis is None and 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
+        ):
             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)
+
+        self._datum = None
+
         #: 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:
@@ -292,37 +338,159 @@ def as_cartopy_projection(self):
         )
 
     def as_cartopy_globe(self):
-        return self._crs.globe
+        return self._globe
 
-    def __getattr__(self, name):
-        if name == "semi_major_axis":
+    @cached_property
+    def _globe(self):
+        """
+        A representation of this CRS as a Cartopy Globe.
+
+        Note
+        ----
+        This property is created when required and then cached for speed. That
+        cached value is cleared when an assignment is made to a property of the
+        class that invalidates the cache.
+        """
+        if self._datum is not None:
+            short_datum = _short_datum_names.get(self._datum, self._datum)
+            # Cartopy doesn't actually enact datums unless they're provided without
+            # ellipsoid axes, so only provide the datum
+            return ccrs.Globe(short_datum, ellipse=None)
+        else:
+            return ccrs.Globe(
+                ellipse=None,
+                semimajor_axis=self._semi_major_axis,
+                semiminor_axis=self._semi_minor_axis,
+            )
+
+    @cached_property
+    def _crs(self):
+        """
+        A representation of this CRS as a Cartopy CRS.
+
+        Note
+        ----
+        This property is created when required and then cached for speed. That
+        cached value is cleared when an assignment is made to a property of the
+        class that invalidates the cache.
+        """
+        return ccrs.Geodetic(self._globe)
+
+    def _wipe_cached_properties(self):
+        """
+        Wipes the cached properties on the object as part of any update to a
+        value that invalidates the cache.
+        """
+        try:
+            delattr(self, "_crs")
+        except AttributeError:
+            pass
+        try:
+            delattr(self, "_globe")
+        except AttributeError:
+            pass
+
+    @property
+    def semi_major_axis(self):
+        if self._semi_major_axis is not None:
+            return self._semi_major_axis
+        else:
             return self._crs.ellipsoid.semi_major_metre
-        if name == "semi_minor_axis":
+
+    @semi_major_axis.setter
+    def semi_major_axis(self, value):
+        """
+        Setting this property to a different value invalidates the current datum
+        (if any) because a datum encodes a specific semi-major axis. This also
+        invalidates the cached `cartopy.Globe` and `cartopy.CRS`.
+        """
+        value = float(value)
+        if not np.isclose(self.semi_major_axis, value):
+            self._datum = None
+            self._wipe_cached_properties()
+        self._semi_major_axis = value
+
+    @property
+    def semi_minor_axis(self):
+        if self._semi_minor_axis is not None:
+            return self._semi_minor_axis
+        else:
             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
+
+    @semi_minor_axis.setter
+    def semi_minor_axis(self, value):
+        """
+        Setting this property to a different value invalidates the current datum
+        (if any) because a datum encodes a specific semi-minor axis. This also
+        invalidates the cached `cartopy.Globe` and `cartopy.CRS`.
+        """
+        value = float(value)
+        if not np.isclose(self.semi_minor_axis, value):
+            self._datum = None
+            self._wipe_cached_properties()
+        self._semi_minor_axis = value
+
+    @property
+    def inverse_flattening(self):
+        if self._inverse_flattening is not None:
+            return self._inverse_flattening
+        else:
+            self._crs.ellipsoid.inverse_flattening
+
+    @inverse_flattening.setter
+    def inverse_flattening(self, value):
+        """
+        Setting this property to a different value does not affect the behaviour
+        of this object any further than the value of this property.
+        """
+        wmsg = (
+            "Setting inverse_flattening does not affect other properties of "
+            "the GeogCS object. To change other properties set them explicitly"
+            " or create a new GeogCS instance."
+        )
+        warnings.warn(wmsg, UserWarning)
+        value = float(value)
+        self._inverse_flattening = value
+
+    @property
+    def datum(self):
+        if self._datum is None:
+            return None
+        else:
+            datum = self._datum
             return datum
-        return getattr(super(), name)
+
+    @datum.setter
+    def datum(self, value):
+        """
+        Setting this property to a different value invalidates the current
+        values of the ellipsoid measurements because a datum encodes its own
+        ellipse. This also invalidates the cached `cartopy.Globe` and
+        `cartopy.CRS`.
+        """
+        if self._datum != value:
+            self._semi_major_axis = None
+            self._semi_minor_axis = None
+            self._inverse_flattening = None
+            self._wipe_cached_properties()
+        self._datum = value
 
     @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))
+
+        crs._semi_major_axis = None
+        crs._semi_minor_axis = None
+        crs._inverse_flattening = None
 
         #: Describes 'zero' on the ellipsoid in degrees.
         crs.longitude_of_prime_meridian = _arg_default(
             longitude_of_prime_meridian, 0
         )
+
+        crs._datum = datum
+
         return crs
 
 

lib/iris/fileformats/_nc_load_rules/helpers.py

@@ -275,7 +275,7 @@ def _get_ellipsoid(cf_grid_var):
             "applied. To apply the datum when loading, use the "
             "iris.FUTURE.datum_support flag."
         )
-        warnings.warn(wmsg, FutureWarning)
+        warnings.warn(wmsg, FutureWarning, stacklevel=14)
         datum = None
 
     if datum is not None:

lib/iris/tests/integration/test_netcdf.py

@@ -661,9 +661,6 @@ def test_save_datum(self):
             iris.save(test_cube, filename)
             with iris.FUTURE.context(datum_support=True):
                 cube = iris.load_cube(filename)
-        print(cube)
-        for coord in cube.coords():
-            print(coord)
 
         test_crs = cube.coord("projection_y_coordinate").coord_system
         actual = str(test_crs.as_cartopy_crs().datum)