Save LAEA grid

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}'

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)

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 = ('non zero false easting ({:.2f}) or '
+               'non zero false northing ({:.2f})'
+               '; unsupported by GRIB Template 3.140'
+               '.').format(false_easting, false_northing)
+        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()