Merge pull request SciTools#1477 from dkillick/grib_save_gdt12

Add GRIB save GDT12 (Transverse Mercator)

Author: pp-mo

lib/iris/fileformats/grib/__init__.py

@@ -960,15 +960,13 @@ def save_grib2(cube, target, append=False, **kwargs):
     else:
         raise ValueError("Can only save grib to filename or writable")
 
-    # discover the lat and lon coords (this bit is common to the pp and grib savers...)
-    lat_coords = filter(lambda coord: "latitude" in coord.name(), cube.coords())
-    lon_coords = filter(lambda coord: "longitude" in coord.name(), cube.coords())
-    if len(lat_coords) != 1 or len(lon_coords) != 1:
-        raise TranslationError("Did not find one (and only one) "
-                               "latitude or longitude coord")
+    x_coords = cube.coords(axis='x', dim_coords=True)
+    y_coords = cube.coords(axis='y', dim_coords=True)
+    if len(x_coords) != 1 or len(y_coords) != 1:
+        raise TranslationError("Did not find one (and only one) x or y coord")
 
     # Save each latlon slice2D in the cube
-    for slice2D in cube.slices([lat_coords[0], lon_coords[0]]):
+    for slice2D in cube.slices([y_coords[0], x_coords[0]]):
 
         # Save this slice to the grib file
         grib_message = gribapi.grib_new_from_samples("GRIB2")

lib/iris/fileformats/grib/_save_rules.py

@@ -35,6 +35,7 @@
 import iris
 import iris.exceptions
 import iris.unit
+from iris.coord_systems import GeogCS, RotatedGeogCS, TransverseMercator
 from iris.fileformats.rules import is_regular, regular_step
 from iris.fileformats.grib import grib_phenom_translation as gptx
 from iris.fileformats.grib._load_convert import (_STATISTIC_TYPE_NAMES,
@@ -92,6 +93,19 @@ def fixup_int32_as_uint32(value):
     return value
 
 
+def ensure_set_int32_value(grib, key, value):
+    """
+    Ensure the workaround function :func:`fixup_int32_as_uint32` is applied as
+    necessary to problem keys.
+
+    """
+    try:
+        gribapi.grib_set(grib, key, value)
+    except gribapi.GribInternalError:
+        value = fixup_int32_as_uint32(value)
+        gribapi.grib_set(grib, key, value)
+
+
 ###############################################################################
 #
 # Constants
@@ -104,13 +118,16 @@ def fixup_int32_as_uint32(value):
 # Reference Regulation 92.1.6
 _DEFAULT_DEGREES_UNITS = 1.0e-6
 
+# Units scaling for Grid Definition Template 3.12
+_TRANSVERSE_MERCATOR_DEGREES_UNITS = 1.0e-2
 
 ###############################################################################
 #
 # Identification Section 1
 #
 ###############################################################################
 
+
 def centre(cube, grib):
     # TODO: read centre from cube
     gribapi.grib_set_long(grib, "centre", 74)  # UKMO
@@ -159,28 +176,35 @@ def identification(cube, grib):
 #
 ###############################################################################
 
-def shape_of_the_earth(cube, grib):
 
+def shape_of_the_earth(cube, grib):
     # assume latlon
     cs = cube.coord(dimensions=[0]).coord_system
 
-    # Turn them all missing to start with (255 for byte, -1 for long)
+    # Initially set shape_of_earth keys to missing (255 for byte, -1 for long).
     gribapi.grib_set_long(grib, "scaleFactorOfRadiusOfSphericalEarth", 255)
     gribapi.grib_set_long(grib, "scaledValueOfRadiusOfSphericalEarth", -1)
     gribapi.grib_set_long(grib, "scaleFactorOfEarthMajorAxis", 255)
     gribapi.grib_set_long(grib, "scaledValueOfEarthMajorAxis", -1)
     gribapi.grib_set_long(grib, "scaleFactorOfEarthMinorAxis", 255)
     gribapi.grib_set_long(grib, "scaledValueOfEarthMinorAxis", -1)
 
-    ellipsoid = cs
-    if isinstance(cs, iris.coord_systems.RotatedGeogCS):
+    if isinstance(cs, GeogCS):
+        ellipsoid = cs
+    else:
         ellipsoid = cs.ellipsoid
+        if ellipsoid is None:
+            msg = "Could not determine shape of the earth from coord system "\
+                  "of horizontal grid."
+            raise iris.exceptions.TranslationError(msg)
 
+    # Spherical earth.
     if ellipsoid.inverse_flattening == 0.0:
         gribapi.grib_set_long(grib, "shapeOfTheEarth", 1)
         gribapi.grib_set_long(grib, "scaleFactorOfRadiusOfSphericalEarth", 0)
         gribapi.grib_set_long(grib, "scaledValueOfRadiusOfSphericalEarth",
                               ellipsoid.semi_major_axis)
+    # Oblate spheroid earth.
     else:
         gribapi.grib_set_long(grib, "shapeOfTheEarth", 7)
         gribapi.grib_set_long(grib, "scaleFactorOfEarthMajorAxis", 0)
@@ -223,9 +247,8 @@ def latlon_first_last(x_coord, y_coord, grib):
 def dx_dy(x_coord, y_coord, grib):
     x_step = regular_step(x_coord)
     y_step = regular_step(y_coord)
-    # TODO: THIS USED BE "Dx" and "Dy"!!! DID THE API CHANGE AGAIN???
-    gribapi.grib_set_double(grib, "DxInDegrees", float(abs(x_step)))
-    gribapi.grib_set_double(grib, "DyInDegrees", float(abs(y_step)))
+    gribapi.grib_set(grib, "DxInDegrees", float(abs(x_step)))
+    gribapi.grib_set(grib, "DyInDegrees", float(abs(y_step)))
 
 
 def scanning_mode_flags(x_coord, y_coord, grib):
@@ -235,7 +258,7 @@ def scanning_mode_flags(x_coord, y_coord, grib):
                           int(y_coord.points[1] - y_coord.points[0] > 0))
 
 
-def latlon_grid_common(cube, grib):
+def horizontal_grid_common(cube, grib):
     # Grib encoding of the sequences of X and Y points.
     y_coord = cube.coord(dimensions=[0])
     x_coord = cube.coord(dimensions=[1])
@@ -312,7 +335,7 @@ def grid_definition_template_0(cube, grib):
     """
     # Constant resolution, aka 'regular' true lat-lon grid.
     gribapi.grib_set_long(grib, "gridDefinitionTemplateNumber", 0)
-    latlon_grid_common(cube, grib)
+    horizontal_grid_common(cube, grib)
     latlon_points_regular(cube, grib)
 
 
@@ -333,7 +356,7 @@ def grid_definition_template_1(cube, grib):
     rotated_pole(cube, grib)
 
     # Encode the lat/lon points.
-    latlon_grid_common(cube, grib)
+    horizontal_grid_common(cube, grib)
     latlon_points_regular(cube, grib)
 
 
@@ -351,19 +374,84 @@ def grid_definition_template_5(cube, grib):
     # Without this, setting "gridDefinitionTemplateNumber" = 5 causes an
     # immediate error.
     # See: https://software.ecmwf.int/issues/browse/SUP-1095
-    # This is acceptable, as the subsequent call to 'latlon_grid_common' will
-    # set these to the correct horizontal dimensions (by calling 'grid_dims').
+    # This is acceptable, as the subsequent call to 'horizontal_grid_common'
+    # will set these to the correct horizontal dimensions
+    # (by calling 'grid_dims').
     gribapi.grib_set(grib, "Ni", 1)
     gribapi.grib_set(grib, "Nj", 1)
     gribapi.grib_set(grib, "gridDefinitionTemplateNumber", 5)
 
     # Record details of the rotated coordinate system.
     rotated_pole(cube, grib)
     # Encode the lat/lon points.
-    latlon_grid_common(cube, grib)
+    horizontal_grid_common(cube, grib)
     latlon_points_irregular(cube, grib)
 
 
+def grid_definition_template_12(cube, grib):
+    """
+    Set keys within the provided grib message based on
+    Grid Definition Template 3.12.
+
+    Template 3.12 is used to represent a Transverse Mercator grid.
+
+    """
+    gribapi.grib_set(grib, "gridDefinitionTemplateNumber", 12)
+
+    # Retrieve some information from the cube.
+    y_coord = cube.coord(dimensions=[0])
+    x_coord = cube.coord(dimensions=[1])
+    cs = y_coord.coord_system
+
+    # Set some keys specific to GDT12.
+    # Encode the horizontal points.
+    x_step = regular_step(x_coord)
+    y_step = regular_step(y_coord)
+    gribapi.grib_set(grib, "Di",
+                     float(abs(x_step)) / _TRANSVERSE_MERCATOR_DEGREES_UNITS)
+    gribapi.grib_set(grib, "Dj",
+                     float(abs(y_step)) / _TRANSVERSE_MERCATOR_DEGREES_UNITS)
+    horizontal_grid_common(cube, grib)
+
+    # GRIBAPI expects unsigned ints in X1, X2, Y1, Y2 but it should accept
+    # signed ints, so work around this.
+    # See https://software.ecmwf.int/issues/browse/SUP-1101
+    ensure_set_int32_value(
+        grib, "Y1",
+        int(y_coord.points[0] / _TRANSVERSE_MERCATOR_DEGREES_UNITS))
+    ensure_set_int32_value(
+        grib, "Y2",
+        int(y_coord.points[-1] / _TRANSVERSE_MERCATOR_DEGREES_UNITS))
+    ensure_set_int32_value(
+        grib, "X1",
+        int(x_coord.points[0] / _TRANSVERSE_MERCATOR_DEGREES_UNITS))
+    ensure_set_int32_value(
+        grib, "X2",
+        int(x_coord.points[-1] / _TRANSVERSE_MERCATOR_DEGREES_UNITS))
+
+    # Lat and lon of reference point are measured in millionths of a degree.
+    gribapi.grib_set(grib, "latitudeOfReferencePoint",
+                     cs.latitude_of_projection_origin / _DEFAULT_DEGREES_UNITS)
+    gribapi.grib_set(grib, "longitudeOfReferencePoint",
+                     cs.longitude_of_central_meridian / _DEFAULT_DEGREES_UNITS)
+
+    # False easting and false northing are measured in units of (10^-2)m.
+    gribapi.grib_set(grib, "XR",
+                     cs.false_easting / _TRANSVERSE_MERCATOR_DEGREES_UNITS)
+    gribapi.grib_set(grib, "YR",
+                     cs.false_northing / _TRANSVERSE_MERCATOR_DEGREES_UNITS)
+
+    # GRIBAPI expects a signed int for scaleFactorAtReferencePoint
+    # but it should accept a float, so work around this.
+    # See https://software.ecmwf.int/issues/browse/SUP-1100
+    value = cs.scale_factor_at_central_meridian
+    key_type = gribapi.grib_get_native_type(grib,
+                                            "scaleFactorAtReferencePoint")
+    if key_type is not float:
+        value = fixup_float32_as_int32(value)
+    gribapi.grib_set(grib, "scaleFactorAtReferencePoint", value)
+
+
 def grid_definition_section(cube, grib):
     """
     Set keys within the grid definition section of the provided grib message,
@@ -375,22 +463,26 @@ def grid_definition_section(cube, grib):
     cs = x_coord.coord_system  # N.B. already checked same cs for x and y.
     regular_x_and_y = is_regular(x_coord) and is_regular(y_coord)
 
-    if isinstance(cs, iris.coord_systems.GeogCS):
+    if isinstance(cs, GeogCS):
         if not regular_x_and_y:
             raise iris.exceptions.TranslationError(
                 'Saving an irregular latlon grid to GRIB (PDT3.4) is not '
                 'yet supported.')
 
         grid_definition_template_0(cube, grib)
 
-    elif isinstance(cs, iris.coord_systems.RotatedGeogCS):
+    elif isinstance(cs, RotatedGeogCS):
         # Rotated coordinate system cases.
         # Choose between GDT 3.1 and 3.5 according to coordinate regularity.
         if regular_x_and_y:
             grid_definition_template_1(cube, grib)
         else:
             grid_definition_template_5(cube, grib)
 
+    elif isinstance(cs, TransverseMercator):
+        # Transverse Mercator coordinate system (template 3.12).
+        grid_definition_template_12(cube, grib)
+
     else:
         raise ValueError('Grib saving is not supported for coordinate system: '
                          '{:s}'.format(cs))

lib/iris/tests/unit/fileformats/grib/save_rules/__init__.py

@@ -39,11 +39,11 @@ def grib_set_trap(grib, name, value):
             # Record a key setting on the mock passed as the 'grib message id'.
             grib.keys[name] = value
 
+        self.mock_gribapi.grib_set = grib_set_trap
         self.mock_gribapi.grib_set_long = grib_set_trap
         self.mock_gribapi.grib_set_float = grib_set_trap
         self.mock_gribapi.grib_set_double = grib_set_trap
         self.mock_gribapi.grib_set_long_array = grib_set_trap
-        self.mock_gribapi.grib_set = grib_set_trap
         self.mock_gribapi.grib_set_array = grib_set_trap
 
         # Create a mock 'grib message id', with a 'keys' dict for settings.