diff --git a/iris_grib/_save_rules.py b/iris_grib/_save_rules.py index 116304c6..da0818a2 100644 --- a/iris_grib/_save_rules.py +++ b/iris_grib/_save_rules.py @@ -22,7 +22,8 @@ import iris from iris.aux_factory import HybridHeightFactory, HybridPressureFactory from iris.coord_systems import (GeogCS, RotatedGeogCS, Mercator, - TransverseMercator, LambertConformal) + TransverseMercator, LambertConformal, + LambertAzimuthalEqualArea) from iris.exceptions import TranslationError @@ -649,6 +650,72 @@ def grid_definition_template_30(cube, grib): eccodes.codes_set(grib, "longitudeOfSouthernPole", 0) +def grid_definition_template_140(cube, grib): + """ + Set keys within the provided grib message based on + Grid Definition Template 3.140. + + Template 3.140 is used to represent a Lambert Azimuthual Equal Area grid. + + """ + + eccodes.codes_set(grib, "gridDefinitionTemplateNumber", 140) + + # Retrieve some information from the cube. + y_coord = cube.coord(dimensions=[0]) + x_coord = cube.coord(dimensions=[1]) + cs = y_coord.coord_system + + # Normalise the coordinate values to millimetres - the resolution + # used in the GRIB message. + y_mm = points_in_unit(y_coord, 'mm') + x_mm = points_in_unit(x_coord, 'mm') + + # Encode the horizontal points. + + # NB. Since we're already in millimetres, our tolerance for + # discrepancy in the differences is 1. + try: + x_step = step(x_mm, atol=1) + y_step = step(y_mm, atol=1) + except ValueError: + msg = ('Irregular coordinates not supported for Lambert ' + 'Azimuthal Equal Area.') + raise TranslationError(msg) + eccodes.codes_set(grib, 'Dx', abs(x_step)) + eccodes.codes_set(grib, 'Dy', abs(y_step)) + + horizontal_grid_common(cube, grib, xy=True) + + # Transform first point into geographic CS + geog = cs.ellipsoid if cs.ellipsoid is not None else GeogCS(1) + first_x, first_y = geog.as_cartopy_crs().transform_point( + x_coord.points[0], + y_coord.points[0], + cs.as_cartopy_crs()) + first_x = first_x % 360 + central_lon = cs.longitude_of_projection_origin % 360 + central_lat = cs.latitude_of_projection_origin + + eccodes.codes_set(grib, "latitudeOfFirstGridPoint", + int(np.round(first_y / _DEFAULT_DEGREES_UNITS))) + eccodes.codes_set(grib, "longitudeOfFirstGridPoint", + int(np.round(first_x / _DEFAULT_DEGREES_UNITS))) + eccodes.codes_set(grib, 'standardParallelInMicrodegrees', + central_lat / _DEFAULT_DEGREES_UNITS) + eccodes.codes_set(grib, 'centralLongitudeInMicrodegrees', + central_lon / _DEFAULT_DEGREES_UNITS) + eccodes.codes_set(grib, 'resolutionAndComponentFlags', + 0x1 << _RESOLUTION_AND_COMPONENTS_GRID_WINDS_BIT) + if (not (np.isclose(cs.false_easting, 0.0, atol=1e-6)) or + not (np.isclose(cs.false_northing, 0.0, atol=1e-6))): + msg = ('non zero false easting ({:.2f}) or ' + 'non zero false northing ({:.2f})' + '; unsupported by GRIB Template 3.140' + '').format(cs.false_easting, cs.false_northing) + raise TranslationError(msg) + + def grid_definition_section(cube, grib): """ Set keys within the grid definition section of the provided grib message, @@ -685,6 +752,11 @@ def grid_definition_section(cube, grib): elif isinstance(cs, LambertConformal): # Lambert Conformal coordinate system (template 3.30). grid_definition_template_30(cube, grib) + + elif isinstance(cs, LambertAzimuthalEqualArea): + # Lambert Azimuthal Equal Area coordinate system (template 3.140). + grid_definition_template_140(cube, grib) + else: name = cs.grid_mapping_name.replace('_', ' ').title() emsg = 'Grib saving is not supported for coordinate system {!r}' diff --git a/iris_grib/tests/integration/save_rules/test_grid_definition_section.py b/iris_grib/tests/integration/save_rules/test_grid_definition_section.py index 64de2775..c1005cc5 100644 --- a/iris_grib/tests/integration/save_rules/test_grid_definition_section.py +++ b/iris_grib/tests/integration/save_rules/test_grid_definition_section.py @@ -18,7 +18,8 @@ Mercator, TransverseMercator, LambertConformal, - AlbersEqualArea) + AlbersEqualArea, + LambertAzimuthalEqualArea) import numpy as np from iris_grib._save_rules import grid_definition_section @@ -100,6 +101,16 @@ def test_grid_definition_template_30(self): grid_definition_section(test_cube, self.mock_grib) self._check_key('gridDefinitionTemplateNumber', 30) + def test_grid_definition_template_140(self): + # Lambert Conformal grid. + x_points = np.arange(3) + y_points = np.arange(3) + coord_units = 'm' + cs = LambertAzimuthalEqualArea(ellipsoid=self.ellipsoid) + test_cube = self._make_test_cube(cs, x_points, y_points, coord_units) + grid_definition_section(test_cube, self.mock_grib) + self._check_key('gridDefinitionTemplateNumber', 140) + def test_coord_system_not_supported(self): # Test an unsupported grid - let's choose Albers Equal Area. x_points = np.arange(3) diff --git a/iris_grib/tests/unit/save_rules/test_grid_definition_template_140.py b/iris_grib/tests/unit/save_rules/test_grid_definition_template_140.py new file mode 100644 index 00000000..04b69117 --- /dev/null +++ b/iris_grib/tests/unit/save_rules/test_grid_definition_template_140.py @@ -0,0 +1,133 @@ +# Copyright iris-grib contributors +# +# This file is part of iris-grib and is released under the LGPL license. +# See COPYING and COPYING.LESSER in the root of the repository for full +# licensing details. +""" +Unit tests for :meth:`iris_grib._save_rules.grid_definition_template_140`. + +""" + +# Import iris_grib.tests first so that some things can be initialised before +# importing anything else. +import iris_grib.tests as tests + +import numpy as np + +import iris.coords +from iris.coord_systems import GeogCS, LambertAzimuthalEqualArea +from iris.exceptions import TranslationError + +from iris_grib._save_rules import ( + grid_definition_template_140 as grid_definition_template, +) +from iris_grib.tests.unit.save_rules import GdtTestMixin + + +class FakeGribError(Exception): + pass + + +class Test(tests.IrisGribTest, GdtTestMixin): + def setUp(self): + self.default_ellipsoid = GeogCS(semi_major_axis=6377563.396, + semi_minor_axis=6356256.909) + self.test_cube = self._make_test_cube() + + GdtTestMixin.setUp(self) + + def _make_test_cube(self, cs=None, x_points=None, y_points=None): + # Create a cube with given properties, or minimal defaults. + if cs is None: + cs = self._default_coord_system() + if x_points is None: + x_points = self._default_x_points() + if y_points is None: + y_points = self._default_y_points() + + x_coord = iris.coords.DimCoord(x_points, 'projection_x_coordinate', + units='m', coord_system=cs) + y_coord = iris.coords.DimCoord(y_points, 'projection_y_coordinate', + units='m', coord_system=cs) + test_cube = iris.cube.Cube(np.zeros((len(y_points), len(x_points)))) + test_cube.add_dim_coord(y_coord, 0) + test_cube.add_dim_coord(x_coord, 1) + return test_cube + + def _default_coord_system(self, false_easting=0, false_northing=0): + return LambertAzimuthalEqualArea(latitude_of_projection_origin=54.9, + longitude_of_projection_origin=-2.5, + false_easting=false_easting, + false_northing=false_northing, + ellipsoid=self.default_ellipsoid) + + def test__template_number(self): + grid_definition_template(self.test_cube, self.mock_grib) + self._check_key('gridDefinitionTemplateNumber', 140) + + def test__shape_of_earth(self): + grid_definition_template(self.test_cube, self.mock_grib) + self._check_key('shapeOfTheEarth', 7) + self._check_key('scaleFactorOfEarthMajorAxis', 0) + self._check_key('scaledValueOfEarthMajorAxis', 6377563.396) + self._check_key('scaleFactorOfEarthMinorAxis', 0) + self._check_key('scaledValueOfEarthMinorAxis', 6356256.909) + + def test__grid_shape(self): + test_cube = self._make_test_cube(x_points=np.arange(13), + y_points=np.arange(6)) + grid_definition_template(test_cube, self.mock_grib) + self._check_key('Nx', 13) + self._check_key('Ny', 6) + + def test__grid_points(self): + test_cube = self._make_test_cube(x_points=[1e6, 3e6, 5e6, 7e6], + y_points=[4e6, 9e6]) + grid_definition_template(test_cube, self.mock_grib) + self._check_key("latitudeOfFirstGridPoint", 81330008) + self._check_key("longitudeOfFirstGridPoint", 98799008) + self._check_key("Dx", 2e9) + self._check_key("Dy", 5e9) + + # specific to grid_definition_template_140 + def test__template_specifics(self): + grid_definition_template(self.test_cube, self.mock_grib) + self._check_key("standardParallelInMicrodegrees", 54900000) + self._check_key("centralLongitudeInMicrodegrees", 357500000) + + def test__scanmode(self): + grid_definition_template(self.test_cube, self.mock_grib) + self._check_key('iScansPositively', 1) + self._check_key('jScansPositively', 1) + + def test__scanmode_reverse(self): + test_cube = self._make_test_cube(x_points=np.arange(7e6, 0, -1e6)) + grid_definition_template(test_cube, self.mock_grib) + self._check_key('iScansPositively', 0) + self._check_key('jScansPositively', 1) + + def __fail_false_easting_northing(self, false_easting, false_northing): + cs = self._default_coord_system(false_easting=false_easting, + false_northing=false_northing) + test_cube = self._make_test_cube(cs=cs) + msg = (r'non zero false easting \(\d*\.\d{2}\) or ' + r'non zero false northing \(\d*\.\d{2}\)' + r'; unsupported by GRIB Template 3\.140' + r'') + with self.assertRaisesRegex( + TranslationError, + msg): + grid_definition_template(test_cube, self.mock_grib) + + def test__fail_false_easting(self): + self.__fail_false_easting_northing(10.0, 0.0) + + def test__fail_false_northing(self): + self.__fail_false_easting_northing(0.0, 10.0) + + def test__fail_false_easting_northing(self): + self.__fail_false_easting_northing(10.0, 10.0) + + +if __name__ == "__main__": + tests.main()