Merge pull request SciTools#1475 from rhattersley/grib2-gdt12

Add load support for GRIB2 GDT 12.

Author: esc24

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)

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:

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

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 <http://www.gnu.org/licenses/>.
+"""
+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()

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 <http://www.gnu.org/licenses/>.
+"""
+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()