diff --git a/docs/src/whatsnew/latest.rst b/docs/src/whatsnew/latest.rst index 1c64610f89..fcb0d05835 100644 --- a/docs/src/whatsnew/latest.rst +++ b/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 diff --git a/lib/iris/__init__.py b/lib/iris/__init__.py index 3e847acad7..7442ff1173 100644 --- a/lib/iris/__init__.py +++ b/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 diff --git a/lib/iris/coord_systems.py b/lib/iris/coord_systems.py index 510aafcb48..1f755998b1 100644 --- a/lib/iris/coord_systems.py +++ b/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. diff --git a/lib/iris/fileformats/_nc_load_rules/helpers.py b/lib/iris/fileformats/_nc_load_rules/helpers.py index 954250b4a2..65ab3d8b5b 100644 --- a/lib/iris/fileformats/_nc_load_rules/helpers.py +++ b/lib/iris/fileformats/_nc_load_rules/helpers.py @@ -19,7 +19,9 @@ import cf_units import numpy as np import numpy.ma as ma +import pyproj +import iris import iris.aux_factory from iris.common.mixin import _get_valid_standard_name import iris.coord_systems @@ -131,6 +133,8 @@ CF_ATTR_BOUNDS = "bounds" CF_ATTR_CALENDAR = "calendar" CF_ATTR_CLIMATOLOGY = "climatology" +CF_ATTR_GRID_CRS_WKT = "crs_wkt" +CF_ATTR_GRID_DATUM = "horizontal_datum_name" CF_ATTR_GRID_INVERSE_FLATTENING = "inverse_flattening" CF_ATTR_GRID_EARTH_RADIUS = "earth_radius" CF_ATTR_GRID_MAPPING_NAME = "grid_mapping_name" @@ -233,7 +237,10 @@ def build_cube_metadata(engine): ################################################################################ def _get_ellipsoid(cf_grid_var): - """Return the ellipsoid definition.""" + """ + Return a :class:`iris.coord_systems.GeogCS` using the relevant properties of + `cf_grid_var`. Returns None if no relevant properties are specified. + """ major = getattr(cf_grid_var, CF_ATTR_GRID_SEMI_MAJOR_AXIS, None) minor = getattr(cf_grid_var, CF_ATTR_GRID_SEMI_MINOR_AXIS, None) inverse_flattening = getattr( @@ -248,21 +255,51 @@ def _get_ellipsoid(cf_grid_var): if major is None and minor is None and inverse_flattening is None: major = getattr(cf_grid_var, CF_ATTR_GRID_EARTH_RADIUS, None) - return major, minor, inverse_flattening + datum = getattr(cf_grid_var, CF_ATTR_GRID_DATUM, None) + # Check crs_wkt if no datum + if datum is None: + crs_wkt = getattr(cf_grid_var, CF_ATTR_GRID_CRS_WKT, None) + if crs_wkt is not None: + proj_crs = pyproj.crs.CRS.from_wkt(crs_wkt) + if proj_crs.datum is not None: + datum = proj_crs.datum.name + + # An unknown crs datum will be treated as None + if datum == "unknown": + datum = None + + if not iris.FUTURE.datum_support: + wmsg = ( + "Ignoring a datum in netCDF load for consistency with existing " + "behaviour. In a future version of Iris, this datum will be " + "applied. To apply the datum when loading, use the " + "iris.FUTURE.datum_support flag." + ) + warnings.warn(wmsg, FutureWarning) + datum = None + + if datum is not None: + return iris.coord_systems.GeogCS.from_datum(datum) + elif major is None and minor is None and inverse_flattening is None: + return None + else: + return iris.coord_systems.GeogCS(major, minor, inverse_flattening) ################################################################################ def build_coordinate_system(engine, cf_grid_var): """Create a coordinate system from the CF-netCDF grid mapping variable.""" - major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) - - return iris.coord_systems.GeogCS(major, minor, inverse_flattening) + coord_system = _get_ellipsoid(cf_grid_var) + if coord_system is None: + raise ValueError("No ellipsoid specified") + else: + return coord_system ################################################################################ def build_rotated_coordinate_system(engine, cf_grid_var): """Create a rotated coordinate system from the CF-netCDF grid mapping variable.""" - major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + ellipsoid = _get_ellipsoid(cf_grid_var) north_pole_latitude = getattr( cf_grid_var, CF_ATTR_GRID_NORTH_POLE_LAT, 90.0 @@ -277,14 +314,6 @@ def build_rotated_coordinate_system(engine, cf_grid_var): cf_grid_var, CF_ATTR_GRID_NORTH_POLE_GRID_LON, 0.0 ) - ellipsoid = None - if ( - major is not None - or minor is not None - or inverse_flattening is not None - ): - ellipsoid = iris.coord_systems.GeogCS(major, minor, inverse_flattening) - rcs = iris.coord_systems.RotatedGeogCS( north_pole_latitude, north_pole_longitude, @@ -302,7 +331,7 @@ def build_transverse_mercator_coordinate_system(engine, cf_grid_var): grid mapping variable. """ - major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + ellipsoid = _get_ellipsoid(cf_grid_var) latitude_of_projection_origin = getattr( cf_grid_var, CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN, None @@ -327,14 +356,6 @@ def build_transverse_mercator_coordinate_system(engine, cf_grid_var): cf_grid_var, CF_ATTR_GRID_SCALE_FACTOR_AT_PROJ_ORIGIN, None ) - ellipsoid = None - if ( - major is not None - or minor is not None - or inverse_flattening is not None - ): - ellipsoid = iris.coord_systems.GeogCS(major, minor, inverse_flattening) - cs = iris.coord_systems.TransverseMercator( latitude_of_projection_origin, longitude_of_central_meridian, @@ -354,7 +375,7 @@ def build_lambert_conformal_coordinate_system(engine, cf_grid_var): grid mapping variable. """ - major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + ellipsoid = _get_ellipsoid(cf_grid_var) latitude_of_projection_origin = getattr( cf_grid_var, CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN, None @@ -368,14 +389,6 @@ def build_lambert_conformal_coordinate_system(engine, cf_grid_var): cf_grid_var, CF_ATTR_GRID_STANDARD_PARALLEL, None ) - ellipsoid = None - if ( - major is not None - or minor is not None - or inverse_flattening is not None - ): - ellipsoid = iris.coord_systems.GeogCS(major, minor, inverse_flattening) - cs = iris.coord_systems.LambertConformal( latitude_of_projection_origin, longitude_of_central_meridian, @@ -395,7 +408,7 @@ def build_stereographic_coordinate_system(engine, cf_grid_var): grid mapping variable. """ - major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + ellipsoid = _get_ellipsoid(cf_grid_var) latitude_of_projection_origin = getattr( cf_grid_var, CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN, None @@ -408,14 +421,6 @@ def build_stereographic_coordinate_system(engine, cf_grid_var): # Iris currently only supports Stereographic projections with a scale # factor of 1.0. This is checked elsewhere. - ellipsoid = None - if ( - major is not None - or minor is not None - or inverse_flattening is not None - ): - ellipsoid = iris.coord_systems.GeogCS(major, minor, inverse_flattening) - cs = iris.coord_systems.Stereographic( latitude_of_projection_origin, longitude_of_projection_origin, @@ -435,7 +440,7 @@ def build_mercator_coordinate_system(engine, cf_grid_var): grid mapping variable. """ - major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + ellipsoid = _get_ellipsoid(cf_grid_var) longitude_of_projection_origin = getattr( cf_grid_var, CF_ATTR_GRID_LON_OF_PROJ_ORIGIN, None @@ -449,14 +454,6 @@ def build_mercator_coordinate_system(engine, cf_grid_var): cf_grid_var, CF_ATTR_GRID_SCALE_FACTOR_AT_PROJ_ORIGIN, None ) - ellipsoid = None - if ( - major is not None - or minor is not None - or inverse_flattening is not None - ): - ellipsoid = iris.coord_systems.GeogCS(major, minor, inverse_flattening) - cs = iris.coord_systems.Mercator( longitude_of_projection_origin, ellipsoid=ellipsoid, @@ -476,7 +473,7 @@ def build_lambert_azimuthal_equal_area_coordinate_system(engine, cf_grid_var): grid mapping variable. """ - major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + ellipsoid = _get_ellipsoid(cf_grid_var) latitude_of_projection_origin = getattr( cf_grid_var, CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN, None @@ -487,14 +484,6 @@ def build_lambert_azimuthal_equal_area_coordinate_system(engine, cf_grid_var): false_easting = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_EASTING, None) false_northing = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_NORTHING, None) - ellipsoid = None - if ( - major is not None - or minor is not None - or inverse_flattening is not None - ): - ellipsoid = iris.coord_systems.GeogCS(major, minor, inverse_flattening) - cs = iris.coord_systems.LambertAzimuthalEqualArea( latitude_of_projection_origin, longitude_of_projection_origin, @@ -513,7 +502,7 @@ def build_albers_equal_area_coordinate_system(engine, cf_grid_var): grid mapping variable. """ - major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + ellipsoid = _get_ellipsoid(cf_grid_var) latitude_of_projection_origin = getattr( cf_grid_var, CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN, None @@ -527,14 +516,6 @@ def build_albers_equal_area_coordinate_system(engine, cf_grid_var): cf_grid_var, CF_ATTR_GRID_STANDARD_PARALLEL, None ) - ellipsoid = None - if ( - major is not None - or minor is not None - or inverse_flattening is not None - ): - ellipsoid = iris.coord_systems.GeogCS(major, minor, inverse_flattening) - cs = iris.coord_systems.AlbersEqualArea( latitude_of_projection_origin, longitude_of_central_meridian, @@ -554,7 +535,7 @@ def build_vertical_perspective_coordinate_system(engine, cf_grid_var): grid mapping variable. """ - major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + ellipsoid = _get_ellipsoid(cf_grid_var) latitude_of_projection_origin = getattr( cf_grid_var, CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN, None @@ -568,14 +549,6 @@ def build_vertical_perspective_coordinate_system(engine, cf_grid_var): false_easting = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_EASTING, None) false_northing = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_NORTHING, None) - ellipsoid = None - if ( - major is not None - or minor is not None - or inverse_flattening is not None - ): - ellipsoid = iris.coord_systems.GeogCS(major, minor, inverse_flattening) - cs = iris.coord_systems.VerticalPerspective( latitude_of_projection_origin, longitude_of_projection_origin, @@ -595,7 +568,7 @@ def build_geostationary_coordinate_system(engine, cf_grid_var): grid mapping variable. """ - major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + ellipsoid = _get_ellipsoid(cf_grid_var) latitude_of_projection_origin = getattr( cf_grid_var, CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN, None @@ -612,14 +585,6 @@ def build_geostationary_coordinate_system(engine, cf_grid_var): cf_grid_var, CF_ATTR_GRID_SWEEP_ANGLE_AXIS, None ) - ellipsoid = None - if ( - major is not None - or minor is not None - or inverse_flattening is not None - ): - ellipsoid = iris.coord_systems.GeogCS(major, minor, inverse_flattening) - cs = iris.coord_systems.Geostationary( latitude_of_projection_origin, longitude_of_projection_origin, diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py index 4818f621c8..2337403f4c 100644 --- a/lib/iris/fileformats/netcdf.py +++ b/lib/iris/fileformats/netcdf.py @@ -2519,6 +2519,8 @@ def add_ellipsoid(ellipsoid): else: cf_var_grid.semi_major_axis = semi_major cf_var_grid.semi_minor_axis = semi_minor + if ellipsoid.datum is not None: + cf_var_grid.horizontal_datum_name = ellipsoid.datum # latlon if isinstance(cs, iris.coord_systems.GeogCS): diff --git a/lib/iris/tests/integration/test_Datums.py b/lib/iris/tests/integration/test_Datums.py new file mode 100755 index 0000000000..77b7f28249 --- /dev/null +++ b/lib/iris/tests/integration/test_Datums.py @@ -0,0 +1,53 @@ +# Copyright Iris contributors +# +# This file is part of Iris and is released under the LGPL license. +# See COPYING and COPYING.LESSER in the root of the repository for full +# licensing details. +"""Integration tests for :class:`iris.coord_systems` datum suppport.""" + +# Import iris.tests first so that some things can be initialised before +# importing anything else. +import iris.tests as tests # isort:skip + +import cartopy.crs as ccrs +import numpy as np + +from iris.coord_systems import GeogCS, LambertConformal + + +class TestDatumTransformation(tests.IrisTest): + def setUp(self): + self.x_points = np.array([-1.5]) + self.y_points = np.array([50.5]) + + self.start_crs = ccrs.OSGB(False) + + def test_transform_points_datum(self): + + # Iris version + wgs84 = GeogCS.from_datum("WGS84") + iris_cs = LambertConformal( + central_lat=54, + central_lon=-4, + secant_latitudes=[52, 56], + ellipsoid=wgs84, + ) + iris_cs_as_cartopy = iris_cs.as_cartopy_crs() + + # Cartopy equivalent + cartopy_cs = ccrs.LambertConformal( + central_latitude=54, + central_longitude=-4, + standard_parallels=[52, 56], + globe=ccrs.Globe("WGS84"), + ) + + expected = cartopy_cs.transform_points( + self.start_crs, self.x_points, self.y_points + ) + + actual = iris_cs_as_cartopy.transform_points( + self.start_crs, self.x_points, self.y_points + ) + + self.assertArrayEqual(expected, actual) diff --git a/lib/iris/tests/integration/test_netcdf.py b/lib/iris/tests/integration/test_netcdf.py index 2a45561e17..8c913a1043 100644 --- a/lib/iris/tests/integration/test_netcdf.py +++ b/lib/iris/tests/integration/test_netcdf.py @@ -24,7 +24,8 @@ import numpy.ma as ma import iris -from iris.coords import CellMethod +import iris.coord_systems +from iris.coords import CellMethod, DimCoord from iris.cube import Cube, CubeList from iris.fileformats.netcdf import ( CF_CONVENTIONS_VERSION, @@ -32,6 +33,7 @@ UnknownCellMethodWarning, ) import iris.tests.stock as stock +import iris.tests.unit.fileformats.netcdf.test_load_cubes as tlc @tests.skip_data @@ -484,6 +486,12 @@ def test_unknown_method(self): @tests.skip_data class TestCoordSystem(tests.IrisTest): + def setUp(self): + tlc.setUpModule() + + def tearDown(self): + tlc.tearDownModule() + def test_load_laea_grid(self): cube = iris.load_cube( tests.get_data_path( @@ -492,6 +500,175 @@ def test_load_laea_grid(self): ) self.assertCML(cube, ("netcdf", "netcdf_laea.cml")) + datum_cf_var_cdl = """ + netcdf output { + dimensions: + y = 4 ; + x = 3 ; + variables: + float data(y, x) ; + data :standard_name = "toa_brightness_temperature" ; + data :units = "K" ; + data :grid_mapping = "mercator" ; + int mercator ; + mercator:grid_mapping_name = "mercator" ; + mercator:longitude_of_prime_meridian = 0. ; + mercator:earth_radius = 6378169. ; + mercator:horizontal_datum_name = "OSGB36" ; + float y(y) ; + y:axis = "Y" ; + y:units = "m" ; + y:standard_name = "projection_y_coordinate" ; + float x(x) ; + x:axis = "X" ; + x:units = "m" ; + x:standard_name = "projection_x_coordinate" ; + + // global attributes: + :Conventions = "CF-1.7" ; + :standard_name_vocabulary = "CF Standard Name Table v27" ; + + data: + + data = + 0, 1, 2, + 3, 4, 5, + 6, 7, 8, + 9, 10, 11 ; + + mercator = _ ; + + y = 1, 2, 3, 5 ; + + x = -6, -4, -2 ; + + } + """ + + datum_wkt_cdl = """ +netcdf output5 { +dimensions: + y = 4 ; + x = 3 ; +variables: + float data(y, x) ; + data :standard_name = "toa_brightness_temperature" ; + data :units = "K" ; + data :grid_mapping = "mercator" ; + int mercator ; + mercator:grid_mapping_name = "mercator" ; + mercator:longitude_of_prime_meridian = 0. ; + mercator:earth_radius = 6378169. ; + mercator:longitude_of_projection_origin = 0. ; + mercator:false_easting = 0. ; + mercator:false_northing = 0. ; + mercator:scale_factor_at_projection_origin = 1. ; + mercator:crs_wkt = "PROJCRS[\\"unknown\\",BASEGEOGCRS[\\"unknown\\",DATUM[\\"OSGB36\\",ELLIPSOID[\\"unknown\\",6378169,0,LENGTHUNIT[\\"metre\\",1,ID[\\"EPSG\\",9001]]]],PRIMEM[\\"Greenwich\\",0,ANGLEUNIT[\\"degree\\",0.0174532925199433],ID[\\"EPSG\\",8901]]],CONVERSION[\\"unknown\\",METHOD[\\"Mercator (variant B)\\",ID[\\"EPSG\\",9805]],PARAMETER[\\"Latitude of 1st standard parallel\\",0,ANGLEUNIT[\\"degree\\",0.0174532925199433],ID[\\"EPSG\\",8823]],PARAMETER[\\"Longitude of natural origin\\",0,ANGLEUNIT[\\"degree\\",0.0174532925199433],ID[\\"EPSG\\",8802]],PARAMETER[\\"False easting\\",0,LENGTHUNIT[\\"metre\\",1],ID[\\"EPSG\\",8806]],PARAMETER[\\"False northing\\",0,LENGTHUNIT[\\"metre\\",1],ID[\\"EPSG\\",8807]]],CS[Cartesian,2],AXIS[\\"(E)\\",east,ORDER[1],LENGTHUNIT[\\"metre\\",1,ID[\\"EPSG\\",9001]]],AXIS[\\"(N)\\",north,ORDER[2],LENGTHUNIT[\\"metre\\",1,ID[\\"EPSG\\",9001]]]]" ; + float y(y) ; + y:axis = "Y" ; + y:units = "m" ; + y:standard_name = "projection_y_coordinate" ; + float x(x) ; + x:axis = "X" ; + x:units = "m" ; + x:standard_name = "projection_x_coordinate" ; + +// global attributes: + :standard_name_vocabulary = "CF Standard Name Table v27" ; + :Conventions = "CF-1.7" ; +data: + + data = + 0, 1, 2, + 3, 4, 5, + 6, 7, 8, + 9, 10, 11 ; + + mercator = _ ; + + y = 1, 2, 3, 5 ; + + x = -6, -4, -2 ; +} + """ + + def test_load_datum_wkt(self): + expected = "OSGB 1936" + nc_path = tlc.cdl_to_nc(self.datum_wkt_cdl) + with iris.FUTURE.context(datum_support=True): + cube = iris.load_cube(nc_path) + test_crs = cube.coord("projection_y_coordinate").coord_system + actual = str(test_crs.as_cartopy_crs().datum) + self.assertStringEqual(expected, actual) + + def test_no_load_datum_wkt(self): + nc_path = tlc.cdl_to_nc(self.datum_wkt_cdl) + with self.assertWarnsRegexp("iris.FUTURE.datum_support"): + cube = iris.load_cube(nc_path) + test_crs = cube.coord("projection_y_coordinate").coord_system + actual = str(test_crs.as_cartopy_crs().datum) + self.assertStringEqual(actual, "unknown") + + def test_load_datum_cf_var(self): + expected = "OSGB 1936" + nc_path = tlc.cdl_to_nc(self.datum_cf_var_cdl) + with iris.FUTURE.context(datum_support=True): + cube = iris.load_cube(nc_path) + test_crs = cube.coord("projection_y_coordinate").coord_system + actual = str(test_crs.as_cartopy_crs().datum) + self.assertStringEqual(expected, actual) + + def test_no_load_datum_cf_var(self): + nc_path = tlc.cdl_to_nc(self.datum_cf_var_cdl) + with self.assertWarnsRegexp("iris.FUTURE.datum_support"): + cube = iris.load_cube(nc_path) + test_crs = cube.coord("projection_y_coordinate").coord_system + actual = str(test_crs.as_cartopy_crs().datum) + self.assertStringEqual(actual, "unknown") + + def test_save_datum(self): + expected = "OSGB 1936" + # saved_crs = iris.coord_systems.GeogCS.from_datum(datum="OSGB36") + saved_crs = iris.coord_systems.Mercator( + ellipsoid=iris.coord_systems.GeogCS.from_datum("OSGB36") + ) + + base_cube = stock.realistic_3d() + base_lat_coord = base_cube.coord("grid_latitude") + test_lat_coord = DimCoord( + base_lat_coord.points, + standard_name="projection_y_coordinate", + coord_system=saved_crs, + ) + base_lon_coord = base_cube.coord("grid_longitude") + test_lon_coord = DimCoord( + base_lon_coord.points, + standard_name="projection_x_coordinate", + coord_system=saved_crs, + ) + test_cube = Cube( + base_cube.data, + standard_name=base_cube.standard_name, + units=base_cube.units, + dim_coords_and_dims=( + (base_cube.coord("time"), 0), + (test_lat_coord, 1), + (test_lon_coord, 2), + ), + ) + + with self.temp_filename(suffix=".nc") as filename: + 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) + self.assertStringEqual(expected, actual) + def _get_scale_factor_add_offset(cube, datatype): """Utility function used by netCDF data packing tests.""" diff --git a/lib/iris/tests/test_coordsystem.py b/lib/iris/tests/test_coordsystem.py index 4229125969..046d76b15a 100644 --- a/lib/iris/tests/test_coordsystem.py +++ b/lib/iris/tests/test_coordsystem.py @@ -216,7 +216,9 @@ def test_as_cartopy_crs(self): cs = GeogCS(6543210, 6500000) res = cs.as_cartopy_crs() globe = ccrs.Globe( - semimajor_axis=6543210.0, semiminor_axis=6500000.0, ellipse=None + semimajor_axis=6543210.0, + semiminor_axis=6500000.0, + ellipse=None, ) expected = ccrs.Geodetic(globe) self.assertEqual(res, expected) @@ -243,7 +245,10 @@ def test_init(self): class Test_RotatedGeogCS_repr(tests.IrisTest): def test_repr(self): rcs = RotatedGeogCS( - 30, 40, north_pole_grid_longitude=50, ellipsoid=GeogCS(6371229) + 30, + 40, + north_pole_grid_longitude=50, + ellipsoid=GeogCS(6371229), ) expected = ( "RotatedGeogCS(30.0, 40.0, " @@ -263,7 +268,10 @@ def test_repr(self): class Test_RotatedGeogCS_str(tests.IrisTest): def test_str(self): rcs = RotatedGeogCS( - 30, 40, north_pole_grid_longitude=50, ellipsoid=GeogCS(6371229) + 30, + 40, + north_pole_grid_longitude=50, + ellipsoid=GeogCS(6371229), ) expected = ( "RotatedGeogCS(30.0, 40.0, " @@ -488,5 +496,23 @@ def test_south_cutoff(self): self.assertEqual(ccrs.cutoff, 30) +class Test_Datums(tests.IrisTest): + def test_default_none(self): + cs = GeogCS(6543210, 6500000) # Arbitrary radii + cartopy_crs = cs.as_cartopy_crs() + self.assertStringEqual(cartopy_crs.datum.name, "unknown") + + def test_set_persist(self): + cs = GeogCS.from_datum(datum="WGS84") + cartopy_crs = cs.as_cartopy_crs() + self.assertStringEqual( + cartopy_crs.datum.name, "World Geodetic System 1984" + ) + + cs = GeogCS.from_datum(datum="OSGB36") + cartopy_crs = cs.as_cartopy_crs() + self.assertStringEqual(cartopy_crs.datum.name, "OSGB 1936") + + if __name__ == "__main__": tests.main() 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 0d3ed932e8..f79aa494f3 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 @@ -11,6 +11,7 @@ import shutil import subprocess import tempfile +import warnings import iris.fileformats._nc_load_rules.engine from iris.fileformats.cf import CFReader @@ -95,10 +96,20 @@ def load_cube_from_cdl(self, cdl_string, cdl_path, nc_path): # Use 'patch' so it is restored after the test. self.patch("iris.fileformats.netcdf.DEBUG", self.debug) - # Call the main translation function to load a single cube. - # _load_cube establishes per-cube facts, activates rules and - # produces an actual cube. - cube = _load_cube(engine, cf, cf_var, nc_path) + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", + message="Ignoring a datum in netCDF load for consistency with existing " + "behaviour. In a future version of Iris, this datum will be " + "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 + # produces an actual cube. + cube = _load_cube(engine, cf, cf_var, nc_path) # Also Record, on the cubes, which hybrid coord elements were identified # by the rules operation. diff --git a/lib/iris/tests/unit/fileformats/netcdf/test_Saver.py b/lib/iris/tests/unit/fileformats/netcdf/test_Saver.py index 6b534b52b5..61b37fe477 100644 --- a/lib/iris/tests/unit/fileformats/netcdf/test_Saver.py +++ b/lib/iris/tests/unit/fileformats/netcdf/test_Saver.py @@ -769,7 +769,7 @@ def check_call(self, coord_name, coord_system, units, expected_units): self.assertEqual(result, expected_units) def test_geogcs_latitude(self): - crs = iris.coord_systems.GeogCS(60, 0) + crs = iris.coord_systems.GeogCS(60, 30) self.check_call( "latitude", coord_system=crs, @@ -778,7 +778,7 @@ def test_geogcs_latitude(self): ) def test_geogcs_longitude(self): - crs = iris.coord_systems.GeogCS(60, 0) + crs = iris.coord_systems.GeogCS(60, 30) self.check_call( "longitude", coord_system=crs, diff --git a/requirements/ci/py38.yml b/requirements/ci/py38.yml index 91cd9d3f5f..4382db1f5e 100644 --- a/requirements/ci/py38.yml +++ b/requirements/ci/py38.yml @@ -18,6 +18,7 @@ dependencies: - netcdf4 - numpy >=1.19 - python-xxhash + - pyproj - scipy # Optional dependencies.