diff --git a/iris_grib/_load_convert.py b/iris_grib/_load_convert.py index f27c551d..927dc753 100644 --- a/iris_grib/_load_convert.py +++ b/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 diff --git a/iris_grib/_save_rules.py b/iris_grib/_save_rules.py index 0f0e4df3..cf0be617 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 @@ -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) diff --git a/iris_grib/tests/unit/load_convert/test_product_definition_6.py b/iris_grib/tests/unit/load_convert/test_product_definition_6.py new file mode 100644 index 00000000..de0f9a3d --- /dev/null +++ b/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() 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..c7d79cd8 --- /dev/null +++ b/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() diff --git a/iris_grib/tests/unit/save_rules/test_product_definition_template_6.py b/iris_grib/tests/unit/save_rules/test_product_definition_template_6.py new file mode 100644 index 00000000..b41bbe5a --- /dev/null +++ b/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()