diff --git a/lib/iris/coord_systems.py b/lib/iris/coord_systems.py index 1f755998b1..9544585535 100644 --- a/lib/iris/coord_systems.py +++ b/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 diff --git a/lib/iris/fileformats/_nc_load_rules/helpers.py b/lib/iris/fileformats/_nc_load_rules/helpers.py index 65ab3d8b5b..34eecdd310 100644 --- a/lib/iris/fileformats/_nc_load_rules/helpers.py +++ b/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: diff --git a/lib/iris/tests/integration/test_netcdf.py b/lib/iris/tests/integration/test_netcdf.py index 8c913a1043..14f1ee59e1 100644 --- a/lib/iris/tests/integration/test_netcdf.py +++ b/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) diff --git a/lib/iris/tests/test_coordsystem.py b/lib/iris/tests/test_coordsystem.py index 046d76b15a..e4c776a063 100644 --- a/lib/iris/tests/test_coordsystem.py +++ b/lib/iris/tests/test_coordsystem.py @@ -224,6 +224,144 @@ def test_as_cartopy_crs(self): self.assertEqual(res, expected) +class Test_GeogCS_equality(tests.IrisTest): + """Test cached values don't break GeogCS equality""" + + def test_as_cartopy_globe(self): + cs_const = GeogCS(6543210, 6500000) + cs_mut = GeogCS(6543210, 6500000) + initial_globe = cs_mut.as_cartopy_globe() + new_globe = cs_mut.as_cartopy_globe() + + self.assertIs(new_globe, initial_globe) + self.assertEqual(cs_const, cs_mut) + + def test_as_cartopy_projection(self): + cs_const = GeogCS(6543210, 6500000) + cs_mut = GeogCS(6543210, 6500000) + initial_projection = cs_mut.as_cartopy_projection() + initial_globe = initial_projection.globe + new_projection = cs_mut.as_cartopy_projection() + new_globe = new_projection.globe + + self.assertIs(new_globe, initial_globe) + self.assertEqual(cs_const, cs_mut) + + def test_as_cartopy_crs(self): + cs_const = GeogCS(6543210, 6500000) + cs_mut = GeogCS(6543210, 6500000) + initial_crs = cs_mut.as_cartopy_crs() + initial_globe = initial_crs.globe + new_crs = cs_mut.as_cartopy_crs() + new_globe = new_crs.globe + + self.assertIs(new_crs, initial_crs) + self.assertIs(new_globe, initial_globe) + self.assertEqual(cs_const, cs_mut) + + def test_update_to_equivalent(self): + cs_const = GeogCS(6500000, 6000000) + # Cause caching + _ = cs_const.as_cartopy_crs() + + cs_mut = GeogCS(6543210, 6000000) + # Cause caching + _ = cs_mut.as_cartopy_crs() + # Set value + cs_mut.semi_major_axis = 6500000 + cs_mut.inverse_flattening = 13 + + self.assertEqual(cs_const.semi_major_axis, 6500000) + self.assertEqual(cs_mut.semi_major_axis, 6500000) + self.assertEqual(cs_const, cs_mut) + + +class Test_GeogCS_mutation(tests.IrisTest): + "Test that altering attributes of a GeogCS instance behaves as expected" + + def test_semi_major_axis_change(self): + # Clear datum + # Clear caches + cs = GeogCS.from_datum("OSGB 1936") + _ = cs.as_cartopy_crs() + self.assertEqual(cs.datum, "OSGB 1936") + cs.semi_major_axis = 6000000 + self.assertIsNone(cs.datum) + self.assertEqual(cs.as_cartopy_globe().semimajor_axis, 6000000) + + def test_semi_major_axis_no_change(self): + # Datum untouched + # Caches untouched + cs = GeogCS.from_datum("OSGB 1936") + initial_crs = cs.as_cartopy_crs() + self.assertEqual(cs.datum, "OSGB 1936") + cs.semi_major_axis = 6377563.396 + self.assertEqual(cs.datum, "OSGB 1936") + new_crs = cs.as_cartopy_crs() + self.assertIs(new_crs, initial_crs) + + def test_semi_minor_axis_change(self): + # Clear datum + # Clear caches + cs = GeogCS.from_datum("OSGB 1936") + _ = cs.as_cartopy_crs() + self.assertEqual(cs.datum, "OSGB 1936") + cs.semi_minor_axis = 6000000 + self.assertIsNone(cs.datum) + self.assertEqual(cs.as_cartopy_globe().semiminor_axis, 6000000) + + def test_semi_minor_axis_no_change(self): + # Datum untouched + # Caches untouched + cs = GeogCS.from_datum("OSGB 1936") + initial_crs = cs.as_cartopy_crs() + self.assertEqual(cs.datum, "OSGB 1936") + cs.semi_minor_axis = 6356256.909237285 + self.assertEqual(cs.datum, "OSGB 1936") + new_crs = cs.as_cartopy_crs() + self.assertIs(new_crs, initial_crs) + + def test_datum_change(self): + # Semi-major axis changes + # All internal ellipoid values set to None + # CRS changes + cs = GeogCS(6543210, 6500000) + _ = cs.as_cartopy_crs() + self.assertTrue("_globe" in cs.__dict__) + self.assertTrue("_crs" in cs.__dict__) + self.assertEqual(cs.semi_major_axis, 6543210) + cs.datum = "OSGB 1936" + self.assertEqual(cs.as_cartopy_crs().datum, "OSGB 1936") + self.assertIsNone(cs.__dict__["_semi_major_axis"]) + self.assertIsNone(cs.__dict__["_semi_minor_axis"]) + self.assertIsNone(cs.__dict__["_inverse_flattening"]) + self.assertEqual(cs.semi_major_axis, 6377563.396) + + def test_datum_no_change(self): + # Caches untouched + cs = GeogCS.from_datum("OSGB 1936") + initial_crs = cs.as_cartopy_crs() + cs.datum = "OSGB 1936" + new_crs = cs.as_cartopy_crs() + self.assertIs(new_crs, initial_crs) + + def test_inverse_flattening_change(self): + # Caches untouched + # Axes unchanged (this behaviour is odd, but matches existing behaviour) + # Warning about lack of effect on other aspects + cs = GeogCS(6543210, 6500000) + initial_crs = cs.as_cartopy_crs() + with self.assertWarnsRegex( + UserWarning, + "Setting inverse_flattening does not affect other properties of the GeogCS object.", + ): + cs.inverse_flattening = cs.inverse_flattening + 1 + new_crs = cs.as_cartopy_crs() + self.assertIs(new_crs, initial_crs) + self.assertEqual(cs.semi_major_axis, 6543210) + self.assertEqual(cs.semi_minor_axis, 6500000) + + class Test_RotatedGeogCS_construction(tests.IrisTest): def test_init(self): rcs = RotatedGeogCS( diff --git a/lib/iris/tests/unit/fileformats/nc_load_rules/actions/__init__.py b/lib/iris/tests/unit/fileformats/nc_load_rules/actions/__init__.py index f79aa494f3..2f627f0a2c 100644 --- a/lib/iris/tests/unit/fileformats/nc_load_rules/actions/__init__.py +++ b/lib/iris/tests/unit/fileformats/nc_load_rules/actions/__init__.py @@ -104,7 +104,6 @@ def load_cube_from_cdl(self, cdl_string, cdl_path, nc_path): "applied. To apply the datum when loading, use the " "iris.FUTURE.datum_support flag.", category=FutureWarning, - module="iris.fileformats._nc_load_rules.helpers", ) # Call the main translation function to load a single cube. # _load_cube establishes per-cube facts, activates rules and