Gdt40 save

pyproject.toml

@@ -123,6 +123,7 @@ disallow_untyped_defs = false  # 964 errors
 disallow_incomplete_defs = false  # 3 errors
 check_untyped_defs = false  # 100 errors
 no_implicit_reexport = false  # 134 errors
+disallow_untyped_decorators = false  # 1 error?
 
 exclude = [
     'noxfile\.py',

src/iris_grib/__init__.py

@@ -19,6 +19,7 @@
 from iris.exceptions import TranslationError
 
 from . import _save_rules
+from ._grib2_convert import TEMPLATE_RECORD
 from ._load_convert import convert as load_convert
 from .message import GribMessage
 
@@ -32,6 +33,7 @@
     # TODO: publish grib1 loading controls, when ready to announce
     # "GRIB1_LOADING_MODE",
     # "Grib1LoadingMode",
+    "TEMPLATE_RECORD",
     "load_cubes",
     "load_pairs_from_fields",
     "save_grib2",
@@ -349,15 +351,35 @@ def save_pairs_from_cube(cube):
         slice of the input and each``field`` is an eccodes message "id".
         N.B. the message "id"s are integer handles.
     """
+    dim_coords = True
     x_coords = cube.coords(axis="x", dim_coords=True)
     y_coords = cube.coords(axis="y", dim_coords=True)
+    if not x_coords and not y_coords:
+        # Try aux-coords if there are no dim-coords : needed for gaussian data
+        # NB in this case we won't rely on axis logic, but simply assume that the
+        #  required coordinates are always "latitude" and "longitude".
+        #  This excludes rotated gaussian grids, but we don't support those (yet).
+        dim_coords = False
+        x_coords = cube.coords("longitude")
+        y_coords = cube.coords("latitude")
+
     if len(x_coords) != 1 or len(y_coords) != 1:
         raise TranslationError("Did not find one (and only one) x or y coord")
 
+    (x_coord,) = x_coords
+    (y_coord,) = y_coords
+
     # Save each latlon slice2D in the cube
-    for slice2D in cube.slices([y_coords[0], x_coords[0]]):
-        grib_message = eccodes.codes_grib_new_from_samples("GRIB2")
-        _save_rules.run(slice2D, grib_message, cube)
+    slice_param = [y_coord, x_coord] if dim_coords else [x_coord]
+    for slice2D in cube.slices(slice_param):
+        if _save_rules.is_grid_definition_template_40(cube, x_coord, y_coord):
+            template_name = "reduced_gg_sfc_grib2"
+        else:
+            template_name = "GRIB2"
+        grib_message = eccodes.codes_grib_new_from_samples(template_name)
+        result = _save_rules.run(slice2D, grib_message, cube, x_coord, y_coord)
+        if result is not None:
+            grib_message = result
         yield (slice2D, grib_message)
 
 

src/iris_grib/_grib1_legacy/grib1_load_rules.py

@@ -8,12 +8,21 @@
 Used exclusively by the 'legacy' GRIB loader in `_grib1_legacy.grib_wrapper`.
 """
 
+import numpy as np
+
+from iris_grib._grib2_convert import ensure_surface_air_pressure_name
+
 from cf_units import CALENDAR_GREGORIAN, Unit
 
 from iris.aux_factory import HybridPressureFactory
 from iris.coords import AuxCoord, CellMethod, DimCoord
 from iris.exceptions import TranslationError
-from iris.fileformats.rules import ConversionMetadata, Factory, Reference
+from iris.fileformats.rules import (
+    ConversionMetadata,
+    Factory,
+    Reference,
+    ReferenceTarget,
+)
 
 
 def grib1_convert(grib):
@@ -168,6 +177,17 @@ def grib1_convert(grib):
         )
         units = "???"
 
+    # Special!!
+    if grib.table2Version == 128 and grib.indicatorOfParameter == 134:
+        standard_name = "air_pressure"
+        units = "Pa"
+        references.append(
+            ReferenceTarget("ref_surface_pressure", ensure_surface_air_pressure_name)
+        )
+        dim_coords_and_dims.append(
+            (DimCoord(np.arange(grib._x_points.shape[0]), long_name="gaussian_grid"), 0)
+        )
+
     ##################################################
     ##### decode TIME -> aux_coords (=always scalar)
 

src/iris_grib/_grib2_convert.py

@@ -11,8 +11,10 @@
 
 from collections import namedtuple
 from collections.abc import Iterable
+from contextlib import contextmanager
 from datetime import datetime, timedelta
 import math
+import threading
 import warnings
 
 import cartopy.crs as ccrs
@@ -1223,6 +1225,13 @@ def grid_definition_template_40_reduced(section, metadata, cs):
     metadata["aux_coords_and_dims"].append((y_coord, 0))
     metadata["aux_coords_and_dims"].append((x_coord, 0))
 
+    # In order to enable a reference surface to be attached to a resulting cube
+    #  with a factory, it is necessary for the cubes to have *also* have a DimCoord on
+    #  the unstructured dimension, so they can match their dimensions.
+    metadata["dim_coords_and_dims"].append(
+        (DimCoord(np.arange(x_points.shape[0]), long_name="gaussian_grid"), 0)
+    )
+
 
 def grid_definition_template_90(section, metadata):
     """
@@ -1379,6 +1388,47 @@ def grid_definition_template_140(section, metadata):
     metadata["dim_coords_and_dims"].append((x_coord, x_dim))
 
 
+class TemplateRecorder(threading.local):
+    """
+    Object to control the recording of grid template number on loading.
+
+    This provides a setting, controlled by `set` and `context` methods, which causes all
+    loaded cubes to record the grid template definition number as a cube attribute.
+
+    When present, the 'GRIB_GRID_TEMPLATE' attribute is used to control how certain
+    types of data are saved -- notably, for gaussian grids.
+
+    e.g. ``TEMPLATE_RECORD.set(True)`` or ``with TEMPLATE_RECORD.context(True): ...``.
+    """
+
+    def __init__(self):
+        self._record = False  # default to 'off' : not recording template numbers
+
+    def __repr__(self):
+        content = ", ".join(f"{key}={value}" for key, value in self.__dict__.items())
+        msg = f"TemplateRecorder({content})"
+        return msg
+
+    def set(self, record: bool = False):
+        self._record = record
+
+    @contextmanager
+    def context(self, record: bool = False):
+        old_value = self._record
+        try:
+            self._record = record
+            yield
+        finally:
+            self._record = old_value
+
+    def __bool__(self):
+        return bool(self._record)
+
+
+#: The unique object which controls recording of grid template numbers.
+TEMPLATE_RECORD = TemplateRecorder()
+
+
 def grid_definition_section(section, metadata):
     """
     Translate section 3 from the GRIB2 message.
@@ -1442,6 +1492,11 @@ def grid_definition_section(section, metadata):
         msg = "Grid definition template [{}] is not supported".format(template)
         raise TranslationError(msg)
 
+    # If enabled, always record the 'original' grib template number, as an attribute.
+    if TEMPLATE_RECORD:
+        # This can also potentially control saving, in some cases.
+        metadata["attributes"]["GRIB2_GRID_TEMPLATE"] = template
+
 
 ###############################################################################
 #
@@ -1527,7 +1582,7 @@ def translate_phenomenon(
     # Look for fields at surface level first.
     if (
         typeOfFirstFixedSurface == 1
-        and scaledValueOfFirstFixedSurface == 0
+        and scaledValueOfFirstFixedSurface in [0, None]
         and typeOfSecondFixedSurface == _TYPE_OF_FIXED_SURFACE_MISSING
     ):
         # Land surface products for model terrain height: