pdt6 & gdt140

iris_grib/_load_convert.py

@@ -2066,6 +2066,37 @@ def product_definition_template_1(section, metadata, frt_coord):
     metadata['aux_coords_and_dims'].append((realization, None))
 
 
+def product_definition_template_6(section, metadata, frt_coord):
+    """
+    Translate template representing percentile forecast,
+    at a horizontal level or in a horizontal layer at a
+    point in time.
+
+    Updates the metadata in-place with the translations.
+
+    Args:
+
+    * section:
+        Dictionary of coded key/value pairs from section 4 of the message.
+
+    * metadata:
+        :class:`collectins.OrderedDict` of metadata.
+
+    * frt_coord:
+        The scalar forecast reference time :class:`iris.coords.DimCoord`.
+
+    """
+    # Perform identical message processing.
+    product_definition_template_0(section, metadata, frt_coord)
+
+    percentile = DimCoord(section['percentileValue'],
+                          long_name='percentile',
+                          units='%')
+
+    # Add the realization coordinate to the metadata aux coords.
+    metadata['aux_coords_and_dims'].append((percentile, None))
+
+
 def product_definition_template_8(section, metadata, frt_coord):
     """
     Translate template representing average, accumulation and/or extreme values

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
 
 
@@ -642,6 +643,72 @@ def grid_definition_template_30(cube, grib):
     gribapi.grib_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 azimuthal equal area
+    projection grid.
+
+    """
+
+    gribapi.grib_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 '
+               'Conformal.')
+        raise TranslationError(msg)
+    gribapi.grib_set(grib, 'Dx', abs(x_step))
+    gribapi.grib_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
+    proj_origin_lon = cs.longitude_of_projection_origin % 360
+
+    gribapi.grib_set(grib, "latitudeOfFirstGridPoint",
+                     int(np.round(first_y / _DEFAULT_DEGREES_UNITS)))
+    gribapi.grib_set(grib, "longitudeOfFirstGridPoint",
+                     int(np.round(first_x / _DEFAULT_DEGREES_UNITS)))
+    gribapi.grib_set(grib, 'resolutionAndComponentFlags',
+                     0x1 << _RESOLUTION_AND_COMPONENTS_GRID_WINDS_BIT)
+    gribapi.grib_set(grib, "standardParallel",
+                     int(np.round(proj_origin_lon
+                                  / _DEFAULT_DEGREES_UNITS)))
+    gribapi.grib_set(grib, "centralLongitude",
+                     int(np.round(cs.longitude_of_projection_origin
+                                  / _DEFAULT_DEGREES_UNITS)))
+    if cs.false_easting != 0.0 or cs.false_northing != 0.0:
+        msg = ('false easting non zero ({}) or false northing non zero ({}) '
+               '\n; 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,
@@ -678,6 +745,8 @@ 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):
+        grid_definition_template_140(cube, grib)
     else:
         name = cs.grid_mapping_name.replace('_', ' ').title()
         emsg = 'Grib saving is not supported for coordinate system {!r}'
@@ -1226,6 +1295,24 @@ def product_definition_template_1(cube, grib, full3d_cube=None):
     set_ensemble(cube, grib)
 
 
+def product_definition_template_6(cube, grib, full3d_cube=None):
+    """
+    Set keys within the provided grib message based on Product Definition
+    Template 4.6.
+
+    Template 4.6 is used to represent a percentile forecast at a point in time
+    interval.
+
+    """
+    gribapi.grib_set(grib, "productDefinitionTemplateNumber", 6)
+    if not (cube.coords('percentile') and
+            len(cube.coord('percentile').points) == 1):
+        raise ValueError("A cube 'percentile' coordinate with one "
+                         "point is required, but not present.")
+    gribapi.grib_set(grib, "percentileValue",
+                     int(cube.coord('percentile').points[0]))
+
+
 def product_definition_template_8(cube, grib, full3d_cube=None):
     """
     Set keys within the provided grib message based on Product
@@ -1420,6 +1507,8 @@ def product_definition_section(cube, grib, full3d_cube=None):
         elif 'spatial_processing_type' in cube.attributes:
             # spatial process (template 4.15)
             product_definition_template_15(cube, grib, full3d_cube)
+        elif cube.coords('percentile'):
+            product_definition_template_6(cube, grib, full3d_cube)
         else:
             # forecast (template 4.0)
             product_definition_template_0(cube, grib, full3d_cube)

iris_grib/tests/unit/load_convert/test_product_definition_6.py

@@ -0,0 +1,67 @@
+# 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.
+"""
+Test function
+:func:`iris_grib._load_convert.product_definition_template_1`.
+
+"""
+
+# import iris_grib.tests first so that some things can be initialised
+# before importing anything else.
+import iris_grib.tests as tests
+
+from copy import deepcopy
+from unittest import mock
+import warnings
+
+from iris.coords import DimCoord
+
+from iris_grib._load_convert import product_definition_template_6
+
+
+class Test(tests.IrisGribTest):
+    def setUp(self):
+        def func(s, m, f):
+            return m['cell_methods'].append(self.cell_method)
+
+        module = 'iris_grib._load_convert'
+        self.patch('warnings.warn')
+        this = '{}.product_definition_template_0'.format(module)
+        self.cell_method = mock.sentinel.cell_method
+        self.patch(this, side_effect=func)
+        self.metadata = {'factories': [], 'references': [],
+                         'standard_name': None,
+                         'long_name': None, 'units': None, 'attributes': {},
+                         'cell_methods': [], 'dim_coords_and_dims': [],
+                         'aux_coords_and_dims': []}
+
+    def _check(self, request_warning):
+        this = 'iris_grib._load_convert.options'
+        with mock.patch(this, warn_on_unsupported=request_warning):
+            metadata = deepcopy(self.metadata)
+            percentile = 50
+            section = {'percentileValue': percentile}
+            forecast_reference_time = mock.sentinel.forecast_reference_time
+            # The called being tested.
+            product_definition_template_6(section, metadata,
+                                          forecast_reference_time)
+            expected = deepcopy(self.metadata)
+            expected['cell_methods'].append(self.cell_method)
+            percentile = DimCoord(percentile,
+                                  long_name='percentile',
+                                  units='%')
+            expected['aux_coords_and_dims'].append((percentile, None))
+            self.assertEqual(metadata, expected)
+
+    def test_pdt_no_warn(self):
+        self._check(False)
+
+    def test_pdt_warn(self):
+        self._check(True)
+
+
+if __name__ == '__main__':
+    tests.main()

iris_grib/tests/unit/save_rules/test_grid_definition_template_140.py

@@ -0,0 +1,107 @@
+# 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
+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=6378137.0,
+                                        semi_minor_axis=6356752.314140356)
+        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):
+        return LambertAzimuthalEqualArea(latitude_of_projection_origin=54.9,
+                                         longitude_of_projection_origin=-2.5,
+                                         false_easting=0.0, false_northing=0.0,
+                                         ellipsoid=self.default_ellipsoid)
+
+    def test__template_number(self):
+        grid_definition_template_140(self.test_cube, self.mock_grib)
+        self._check_key('gridDefinitionTemplateNumber', 140)
+
+    def test__shape_of_earth(self):
+        grid_definition_template_140(self.test_cube, self.mock_grib)
+        self._check_key('shapeOfTheEarth', 7)
+        self._check_key('scaleFactorOfEarthMajorAxis', 0)
+        self._check_key('scaledValueOfEarthMajorAxis', 6378137.0)
+        self._check_key('scaleFactorOfEarthMinorAxis', 0)
+        self._check_key('scaledValueOfEarthMinorAxis', 6356752.314140356)
+
+    def test__grid_shape(self):
+        test_cube = self._make_test_cube(x_points=np.arange(13),
+                                         y_points=np.arange(6))
+        grid_definition_template_140(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_140(test_cube, self.mock_grib)
+        self._check_key("latitudeOfFirstGridPoint", 81331520)
+        self._check_key("longitudeOfFirstGridPoint", 98776029)
+        self._check_key("Dx", 2e9)
+        self._check_key("Dy", 5e9)
+
+    def test__template_specifics(self):
+        grid_definition_template_140(self.test_cube, self.mock_grib)
+        self._check_key("standardParallel", 357500000)
+        self._check_key("centralLongitude", -2500000)
+
+    def test__scanmode(self):
+        grid_definition_template_140(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_140(test_cube, self.mock_grib)
+        self._check_key('iScansPositively', 0)
+        self._check_key('jScansPositively', 1)
+
+
+if __name__ == "__main__":
+    tests.main()

iris_grib/tests/unit/save_rules/test_product_definition_template_6.py

@@ -0,0 +1,59 @@
+# 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 :func:`iris_grib._save_rules.product_definition_template_6`
+
+"""
+
+# Import iris_grib.tests first so that some things can be initialised before
+# importing anything else.
+import iris_grib.tests as tests
+
+from unittest import mock
+
+from cf_units import Unit
+import gribapi
+
+from iris.coords import CellMethod, DimCoord
+import iris.tests.stock as stock
+
+from iris_grib._save_rules import product_definition_template_6
+
+
+class TestRealizationIdentifier(tests.IrisGribTest):
+    def setUp(self):
+        self.cube = stock.lat_lon_cube()
+        # Rename cube to avoid warning about unknown discipline/parameter.
+        self.cube.rename('air_temperature')
+        coord = DimCoord([45], 'time',
+                         units=Unit('days since epoch', calendar='standard'))
+        self.cube.add_aux_coord(coord)
+
+    @mock.patch.object(gribapi, 'grib_set')
+    def test_percentile(self, mock_set):
+        cube = self.cube
+        coord = DimCoord(10, long_name='percentile', units='%')
+        cube.add_aux_coord(coord)
+
+        product_definition_template_6(cube, mock.sentinel.grib)
+        mock_set.assert_any_call(mock.sentinel.grib,
+                                 "productDefinitionTemplateNumber", 6)
+        mock_set.assert_any_call(mock.sentinel.grib,
+                                 "percentileValue", 10)
+
+    @mock.patch.object(gribapi, 'grib_set')
+    def test_multiple_percentile_values(self, mock_set):
+        cube = self.cube
+        coord = DimCoord([8, 9, 10], long_name='percentile', units='%')
+        cube.add_aux_coord(coord, 0)
+
+        msg = "'percentile' coordinate with one point is required"
+        with self.assertRaisesRegex(ValueError, msg):
+            product_definition_template_6(cube, mock.sentinel.grib)
+
+
+if __name__ == "__main__":
+    tests.main()