diff --git a/lib/iris/fileformats/grib/_load_convert.py b/lib/iris/fileformats/grib/_load_convert.py index 03b226aaee..f84b1d49cf 100644 --- a/lib/iris/fileformats/grib/_load_convert.py +++ b/lib/iris/fileformats/grib/_load_convert.py @@ -214,6 +214,40 @@ def _hindcast_fix(forecast_time): return forecast_time +def fixup_float32_from_int32(value): + """ + Workaround for use when reading an IEEE 32-bit floating-point value + which the ECMWF GRIB API has erroneously treated as a 4-byte signed + integer. + + """ + # Convert from two's complement to sign-and-magnitude. + # NB. The bit patterns 0x00000000 and 0x80000000 will both be + # returned by the ECMWF GRIB API as an integer 0. Because they + # correspond to positive and negative zero respectively it is safe + # to treat an integer 0 as a positive zero. + if value < 0: + value = 0x80000000 - value + value_as_uint32 = np.array(value, dtype='u4') + value_as_float32 = value_as_uint32.view(dtype='f4') + return float(value_as_float32) + + +def fixup_int32_from_uint32(value): + """ + Workaround for use when reading a signed, 4-byte integer which the + ECMWF GRIB API has erroneously treated as an unsigned, 4-byte + integer. + + NB. This workaround is safe to use with values which are already + treated as signed, 4-byte integers. + + """ + if value >= 0x80000000: + value = 0x80000000 - value + return value + + ############################################################################### # # Identification Section 1 @@ -644,6 +678,83 @@ def grid_definition_template_5(section, metadata): 'grid_latitude', 'grid_longitude', cs) +def grid_definition_template_12(section, metadata): + """ + Translate template representing transverse Mercator. + + Updates the metadata in-place with the translations. + + Args: + + * section: + Dictionary of coded key/value pairs from section 3 of the message. + + * metadata: + :class:`collections.OrderedDict` of metadata. + + """ + major, minor, radius = ellipsoid_geometry(section) + geog_cs = ellipsoid(section['shapeOfTheEarth'], major, minor, radius) + + lat = section['latitudeOfReferencePoint'] * _GRID_ACCURACY_IN_DEGREES + lon = section['longitudeOfReferencePoint'] * _GRID_ACCURACY_IN_DEGREES + scale = section['scaleFactorAtReferencePoint'] + # Catch bug in ECMWF GRIB API (present at 1.12.1) where the scale + # is treated as a signed, 4-byte integer. + if isinstance(scale, int): + scale = fixup_float32_from_int32(scale) + CM_TO_M = 0.01 + easting = section['XR'] * CM_TO_M + northing = section['YR'] * CM_TO_M + cs = icoord_systems.TransverseMercator(lat, lon, easting, northing, + scale, geog_cs) + + # Deal with bug in ECMWF GRIB API (present at 1.12.1) where these + # values are treated as unsigned, 4-byte integers. + x1 = fixup_int32_from_uint32(section['x1']) + y1 = fixup_int32_from_uint32(section['y1']) + x2 = fixup_int32_from_uint32(section['x2']) + y2 = fixup_int32_from_uint32(section['y2']) + + # Rather unhelpfully this grid definition template seems to be + # overspecified, and thus open to inconsistency. + last_x = x1 + (section['Ni'] - 1) * section['Di'] + last_y = y1 + (section['Nj'] - 1) * section['Dj'] + if (last_x != x2 or last_y != y2): + raise TranslationError('Inconsistent grid definition') + + x1 = x1 * CM_TO_M + dx = section['Di'] * CM_TO_M + x_points = x1 + np.arange(section['Ni']) * dx + y1 = y1 * CM_TO_M + dy = section['Dj'] * CM_TO_M + y_points = y1 + np.arange(section['Nj']) * dy + + # This has only been tested with +x/+y scanning, so raise an error + # for other permutations. + scan = scanning_mode(section['scanningMode']) + if scan.i_negative: + raise TranslationError('Unsupported -x scanning') + if not scan.j_positive: + raise TranslationError('Unsupported -y scanning') + + # Create the X and Y coordinates. + y_coord = DimCoord(y_points, 'projection_y_coordinate', units='m', + coord_system=cs) + x_coord = DimCoord(x_points, 'projection_x_coordinate', units='m', + coord_system=cs) + + # Determine the lat/lon dimensions. + y_dim, x_dim = 0, 1 + scan = scanning_mode(section['scanningMode']) + if scan.j_consecutive: + y_dim, x_dim = 1, 0 + + # Add the X and Y coordinates to the metadata dim coords. + metadata['dim_coords_and_dims'].append((y_coord, y_dim)) + metadata['dim_coords_and_dims'].append((x_coord, x_dim)) + + def grid_definition_template_90(section, metadata): """ Translate template representing space view. @@ -792,6 +903,9 @@ def grid_definition_section(section, metadata): elif template == 5: # Process variable resolution rotated latitude/longitude. grid_definition_template_5(section, metadata) + elif template == 12: + # Process transverse Mercator. + grid_definition_template_12(section, metadata) elif template == 90: # Process space view. grid_definition_template_90(section, metadata) diff --git a/lib/iris/fileformats/grib/_message.py b/lib/iris/fileformats/grib/_message.py index 2e85f23d91..bb0d8b573e 100644 --- a/lib/iris/fileformats/grib/_message.py +++ b/lib/iris/fileformats/grib/_message.py @@ -101,7 +101,7 @@ def data(self): 'unsupported quasi-regular grid.') template = grid_section['gridDefinitionTemplateNumber'] - if template in (0, 1, 5, 90): + if template in (0, 1, 5, 12, 90): # We can ignore the first two bits (i-neg, j-pos) because # that is already captured in the coordinate values. if grid_section['scanningMode'] & 0x3f: diff --git a/lib/iris/tests/unit/fileformats/grib/load_convert/__init__.py b/lib/iris/tests/unit/fileformats/grib/load_convert/__init__.py index 45b0743731..16ae7153c2 100644 --- a/lib/iris/tests/unit/fileformats/grib/load_convert/__init__.py +++ b/lib/iris/tests/unit/fileformats/grib/load_convert/__init__.py @@ -17,3 +17,19 @@ """Unit tests for the :mod:`iris.fileformats.grib._load_convert` package.""" from __future__ import (absolute_import, division, print_function) + +from collections import OrderedDict + + +def empty_metadata(): + metadata = OrderedDict() + metadata['factories'] = [] + metadata['references'] = [] + metadata['standard_name'] = None + metadata['long_name'] = None + metadata['units'] = None + metadata['attributes'] = {} + metadata['cell_methods'] = [] + metadata['dim_coords_and_dims'] = [] + metadata['aux_coords_and_dims'] = [] + return metadata diff --git a/lib/iris/tests/unit/fileformats/grib/load_convert/test_fixup_float32_from_int32.py b/lib/iris/tests/unit/fileformats/grib/load_convert/test_fixup_float32_from_int32.py new file mode 100644 index 0000000000..ee2b14e06b --- /dev/null +++ b/lib/iris/tests/unit/fileformats/grib/load_convert/test_fixup_float32_from_int32.py @@ -0,0 +1,46 @@ +# (C) British Crown Copyright 2014, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +""" +Unit tests for `iris.fileformats.grib._load_convert.fixup_float32_from_int32`. + +""" + +from __future__ import (absolute_import, division, print_function) + +# Import iris.tests first so that some things can be initialised before +# importing anything else. +import iris.tests as tests + +from iris.fileformats.grib._load_convert import fixup_float32_from_int32 + + +class Test(tests.IrisTest): + def test_negative(self): + result = fixup_float32_from_int32(-0x3f000000) + self.assertEqual(result, -0.5) + + def test_zero(self): + result = fixup_float32_from_int32(0) + self.assertEqual(result, 0) + + def test_positive(self): + result = fixup_float32_from_int32(0x3f000000) + self.assertEqual(result, 0.5) + + +if __name__ == '__main__': + tests.main() diff --git a/lib/iris/tests/unit/fileformats/grib/load_convert/test_fixup_int32_from_uint32.py b/lib/iris/tests/unit/fileformats/grib/load_convert/test_fixup_int32_from_uint32.py new file mode 100644 index 0000000000..c603b6b9fb --- /dev/null +++ b/lib/iris/tests/unit/fileformats/grib/load_convert/test_fixup_int32_from_uint32.py @@ -0,0 +1,56 @@ +# (C) British Crown Copyright 2014, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +""" +Unit tests for `iris.fileformats.grib._load_convert.fixup_int32_from_uint32`. + +""" + +from __future__ import (absolute_import, division, print_function) + +# Import iris.tests first so that some things can be initialised before +# importing anything else. +import iris.tests as tests + +from iris.fileformats.grib._load_convert import fixup_int32_from_uint32 + + +class Test(tests.IrisTest): + def test_negative(self): + result = fixup_int32_from_uint32(0x80000005) + self.assertEqual(result, -5) + + def test_negative_zero(self): + result = fixup_int32_from_uint32(0x80000000) + self.assertEqual(result, 0) + + def test_zero(self): + result = fixup_int32_from_uint32(0) + self.assertEqual(result, 0) + + def test_positive(self): + result = fixup_int32_from_uint32(200000) + self.assertEqual(result, 200000) + + def test_already_negative(self): + # If we *already* have a negative value the fixup routine should + # leave it alone. + result = fixup_int32_from_uint32(-7) + self.assertEqual(result, -7) + + +if __name__ == '__main__': + tests.main() diff --git a/lib/iris/tests/unit/fileformats/grib/load_convert/test_grid_definition_template_12.py b/lib/iris/tests/unit/fileformats/grib/load_convert/test_grid_definition_template_12.py new file mode 100644 index 0000000000..15d2f6a453 --- /dev/null +++ b/lib/iris/tests/unit/fileformats/grib/load_convert/test_grid_definition_template_12.py @@ -0,0 +1,153 @@ +# (C) British Crown Copyright 2014, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +""" +Unit tests for +:func:`iris.fileformats.grib._load_convert.grid_definition_template_12`. + +""" + +from __future__ import (absolute_import, division, print_function) + +# import iris tests first so that some things can be initialised +# before importing anything else. +import iris.tests as tests + +import numpy as np + +import iris.coord_systems +import iris.coords +import iris.exceptions +from iris.tests.unit.fileformats.grib.load_convert import empty_metadata +from iris.fileformats.grib._load_convert import grid_definition_template_12 + + +MDI = 2 ** 32 - 1 + + +class Test(tests.IrisTest): + def section_3(self): + section = { + 'shapeOfTheEarth': 7, + 'scaleFactorOfRadiusOfSphericalEarth': MDI, + 'scaledValueOfRadiusOfSphericalEarth': MDI, + 'scaleFactorOfEarthMajorAxis': 3, + 'scaledValueOfEarthMajorAxis': 6377563396, + 'scaleFactorOfEarthMinorAxis': 3, + 'scaledValueOfEarthMinorAxis': 6356256909, + 'Ni': 4, + 'Nj': 3, + 'latitudeOfReferencePoint': 49000000, + 'longitudeOfReferencePoint': -2000000, + 'resolutionAndComponentFlags': 0, + 'scaleFactorAtReferencePoint': 0.9996012717, + 'XR': 40000000, + 'YR': -10000000, + 'scanningMode': 64, + 'Di': 200000, + 'Dj': 100000, + 'x1': 29300000, + 'y1': 9200000, + 'x2': 29900000, + 'y2': 9400000 + } + return section + + def expected(self, y_dim, x_dim): + # Prepare the expectation. + expected = empty_metadata() + ellipsoid = iris.coord_systems.GeogCS(6377563.396, 6356256.909) + cs = iris.coord_systems.TransverseMercator(49, -2, 400000, -100000, + 0.9996012717, ellipsoid) + nx = 4 + x_origin = 293000 + dx = 2000 + x = iris.coords.DimCoord(np.arange(nx) * dx + x_origin, + 'projection_x_coordinate', units='m', + coord_system=cs) + ny = 3 + y_origin = 92000 + dy = 1000 + y = iris.coords.DimCoord(np.arange(ny) * dy + y_origin, + 'projection_y_coordinate', units='m', + coord_system=cs) + expected['dim_coords_and_dims'].append((y, y_dim)) + expected['dim_coords_and_dims'].append((x, x_dim)) + return expected + + def test(self): + section = self.section_3() + metadata = empty_metadata() + grid_definition_template_12(section, metadata) + expected = self.expected(0, 1) + self.assertEqual(metadata, expected) + + def test_spherical(self): + section = self.section_3() + section['shapeOfTheEarth'] = 0 + metadata = empty_metadata() + grid_definition_template_12(section, metadata) + expected = self.expected(0, 1) + cs = expected['dim_coords_and_dims'][0][0].coord_system + cs.ellipsoid = iris.coord_systems.GeogCS(6367470) + self.assertEqual(metadata, expected) + + def test_negative_x(self): + section = self.section_3() + section['scanningMode'] = 0b11000000 + metadata = empty_metadata() + with self.assertRaisesRegexp(iris.exceptions.TranslationError, + '-x scanning'): + grid_definition_template_12(section, metadata) + + def test_negative_y(self): + section = self.section_3() + section['scanningMode'] = 0b00000000 + metadata = empty_metadata() + with self.assertRaisesRegexp(iris.exceptions.TranslationError, + '-y scanning'): + grid_definition_template_12(section, metadata) + + def test_transposed(self): + section = self.section_3() + section['scanningMode'] = 0b01100000 + metadata = empty_metadata() + grid_definition_template_12(section, metadata) + expected = self.expected(1, 0) + self.assertEqual(metadata, expected) + + def test_incompatible_grid_extent(self): + section = self.section_3() + section['x2'] += 100 + metadata = empty_metadata() + with self.assertRaisesRegexp(iris.exceptions.TranslationError, + 'grid'): + grid_definition_template_12(section, metadata) + + def test_scale_workaround(self): + section = self.section_3() + section['scaleFactorAtReferencePoint'] = 1065346526 + metadata = empty_metadata() + grid_definition_template_12(section, metadata) + expected = self.expected(0, 1) + # A float32 can't hold exactly the same value. + cs = expected['dim_coords_and_dims'][0][0].coord_system + cs.scale_factor_at_central_meridian = 0.9996012449264526 + self.assertEqual(metadata, expected) + + +if __name__ == '__main__': + tests.main() diff --git a/lib/iris/tests/unit/fileformats/grib/load_convert/test_grid_definition_template_90.py b/lib/iris/tests/unit/fileformats/grib/load_convert/test_grid_definition_template_90.py index e385c4f06b..48256ad04f 100644 --- a/lib/iris/tests/unit/fileformats/grib/load_convert/test_grid_definition_template_90.py +++ b/lib/iris/tests/unit/fileformats/grib/load_convert/test_grid_definition_template_90.py @@ -26,13 +26,12 @@ # before importing anything else. import iris.tests as tests -from collections import OrderedDict - import numpy as np import iris.coord_systems import iris.coords import iris.exceptions +from iris.tests.unit.fileformats.grib.load_convert import empty_metadata from iris.fileformats.grib._load_convert import grid_definition_template_90 @@ -40,19 +39,6 @@ class Test(tests.IrisTest): - def empty_metadata(self): - metadata = OrderedDict() - metadata['factories'] = [] - metadata['references'] = [] - metadata['standard_name'] = None - metadata['long_name'] = None - metadata['units'] = None - metadata['attributes'] = {} - metadata['cell_methods'] = [] - metadata['dim_coords_and_dims'] = [] - metadata['aux_coords_and_dims'] = [] - return metadata - def uk(self): section = { 'shapeOfTheEarth': 3, @@ -81,7 +67,7 @@ def uk(self): def expected_uk(self, y_dim, x_dim): # Prepare the expectation. - expected = self.empty_metadata() + expected = empty_metadata() major = 6378168.8 ellipsoid = iris.coord_systems.GeogCS(major, 6356584.0) height = (6610674e-6 - 1) * major @@ -147,7 +133,7 @@ def compare(self, metadata, expected): def test_uk(self): section = self.uk() - metadata = self.empty_metadata() + metadata = empty_metadata() grid_definition_template_90(section, metadata) expected = self.expected_uk(0, 1) self.compare(metadata, expected) @@ -155,7 +141,7 @@ def test_uk(self): def test_uk_transposed(self): section = self.uk() section['scanningMode'] = 0b11100000 - metadata = self.empty_metadata() + metadata = empty_metadata() grid_definition_template_90(section, metadata) expected = self.expected_uk(1, 0) self.compare(metadata, expected) @@ -163,7 +149,7 @@ def test_uk_transposed(self): def test_non_zero_latitude(self): section = self.uk() section['latitudeOfSubSatellitePoint'] = 1 - metadata = self.empty_metadata() + metadata = empty_metadata() with self.assertRaisesRegexp(iris.exceptions.TranslationError, 'non-zero latitude'): grid_definition_template_90(section, metadata) @@ -171,7 +157,7 @@ def test_non_zero_latitude(self): def test_rotated_meridian(self): section = self.uk() section['orientationOfTheGrid'] = 1 - metadata = self.empty_metadata() + metadata = empty_metadata() with self.assertRaisesRegexp(iris.exceptions.TranslationError, 'orientation'): grid_definition_template_90(section, metadata) @@ -179,7 +165,7 @@ def test_rotated_meridian(self): def test_zero_height(self): section = self.uk() section['Nr'] = 0 - metadata = self.empty_metadata() + metadata = empty_metadata() with self.assertRaisesRegexp(iris.exceptions.TranslationError, 'zero'): grid_definition_template_90(section, metadata) @@ -187,7 +173,7 @@ def test_zero_height(self): def test_orthographic(self): section = self.uk() section['Nr'] = MDI - metadata = self.empty_metadata() + metadata = empty_metadata() with self.assertRaisesRegexp(iris.exceptions.TranslationError, 'orthographic'): grid_definition_template_90(section, metadata) @@ -195,14 +181,14 @@ def test_orthographic(self): def test_scanning_mode_positive_x(self): section = self.uk() section['scanningMode'] = 0b01000000 - metadata = self.empty_metadata() + metadata = empty_metadata() with self.assertRaisesRegexp(iris.exceptions.TranslationError, r'\+x'): grid_definition_template_90(section, metadata) def test_scanning_mode_negative_y(self): section = self.uk() section['scanningMode'] = 0b10000000 - metadata = self.empty_metadata() + metadata = empty_metadata() with self.assertRaisesRegexp(iris.exceptions.TranslationError, '-y'): grid_definition_template_90(section, metadata) diff --git a/lib/iris/tests/unit/fileformats/grib/message/test__GribMessage.py b/lib/iris/tests/unit/fileformats/grib/message/test__GribMessage.py index c10af1e12d..0c42b96bdc 100644 --- a/lib/iris/tests/unit/fileformats/grib/message/test__GribMessage.py +++ b/lib/iris/tests/unit/fileformats/grib/message/test__GribMessage.py @@ -214,6 +214,11 @@ def section_3(self, scanning_mode): return _example_section_3(5, scanning_mode) +class Test_data__grid_template_12(tests.IrisTest, Mixin_data__grid_template): + def section_3(self, scanning_mode): + return _example_section_3(12, scanning_mode) + + class Test_data__grid_template_90(tests.IrisTest, Mixin_data__grid_template): def section_3(self, scanning_mode): section_3 = _example_section_3(90, scanning_mode)