Merge pull request SciTools#3 from kaedonkers/save_rules_prod_def_15

Save rules prod def 15

iris_grib/_load_convert.py

@@ -143,7 +143,7 @@
     5: 'Spectral interpolation using the 4 source grid grid-point values '
        'surrounding the nominal grid-point',
     6: 'Neighbour-budget interpolation using the 4 source grid grid-point '
-       'values surrounding the nominal gridpoint'
+       'values surrounding the nominal grid-point'
 }
 
 # Class containing details of a probability analysis.
@@ -2178,25 +2178,31 @@ def product_definition_template_15(section, metadata, frt_coord):
     # Check unique keys for this template.
     spatial_processing_code = section['spatialProcessing']
 
-    if spatial_processing_code != 0:
-        # For now, we only support the simplest case, representing a statistic
-        # over the whole notional area of a cell.
+    # Only a limited number of spatial processing codes are supported
+    if spatial_processing_code not in _SPATIAL_PROCESSING_TYPES.keys():
         msg = ('Product definition section 4 contains an unsupported '
                'spatial processing type [{}]'.format(spatial_processing_code))
         raise TranslationError(msg)
 
-    # NOTE: PDT 4.15 alse defines a 'numberOfPointsUsed' key, but we think this
-    # is irrelevant to the currently supported spatial-processing types.
+    # NOTE: PDT 4.15 also defines a 'numberOfPointsUsed' key, but we have
+    # not developed a way to interpret this in a CF-compliant way and code
+    # table 4.15 includes this information in each description of types 0-6.
 
-    # Process parts in common with pdt 4.0.
+    # Process parts in common with PDT 4.0.
     product_definition_template_0(section, metadata, frt_coord)
 
-    # Decode the statistic method name.
-    cell_method_name = statistical_method_name(section)
+    # Add spatial processing type as an attribute.
+    metadata['attributes']['spatial_processing_type'] = spatial_processing_code
 
-    # Record an 'area' cell-method using this statistic.
-    metadata['cell_methods'] = [CellMethod(coords=('area',),
-                                           method=cell_method_name)]
+    # Only statistics on non-interpolated data can currently be represented
+    # with cell-methods.
+    if spatial_processing_code == 0:
+        # Decode the statistical method name.
+        cell_method_name = statistical_method_name(section)
+
+        # Record an 'area' cell-method using this statistic.
+        metadata['cell_methods'] = [CellMethod(coords=('area',),
+                                               method=cell_method_name)]
 
 
 def satellite_common(section, metadata):
@@ -2317,7 +2323,7 @@ def product_definition_template_40(section, metadata, frt_coord):
     # Reference GRIB2 Code Table 4.230.
     constituent_type = section['constituentType']
 
-    # Add the constituent type as  an attribute.
+    # Add the constituent type as an attribute.
     metadata['attributes']['WMO_constituent_type'] = constituent_type
 
 

iris_grib/_save_rules.py

@@ -41,17 +41,14 @@
 
 from ._iris_mercator_support import confirm_extended_mercator_supported
 from . import grib_phenom_translation as gptx
-from ._load_convert import (_STATISTIC_TYPE_NAMES, _TIME_RANGE_UNITS,
-                            _SPATIAL_PROCESSING_TYPES)
+from ._load_convert import (_STATISTIC_TYPE_NAMES, _TIME_RANGE_UNITS)
 from iris.util import is_regular, regular_step
 
 
 # Invert code tables from :mod:`iris_grib._load_convert`.
 _STATISTIC_TYPE_NAMES = {val: key for key, val in
                          _STATISTIC_TYPE_NAMES.items()}
 _TIME_RANGE_UNITS = {val: key for key, val in _TIME_RANGE_UNITS.items()}
-_SPATIAL_PROCESSING_TYPES = {val: key for key, val in
-                             _SPATIAL_PROCESSING_TYPES.items()}
 
 
 def fixup_float32_as_int32(value):
@@ -1210,7 +1207,7 @@ def _product_definition_template_8_10_and_11(cube, grib, full3d_cube=None):
     Set keys within the provided grib message based on common aspects of
     Product Definition Templates 4.8 and 4.11.
 
-    Templates 4.8  and 4.11 are used to represent aggregations over a time
+    Templates 4.8 and 4.11 are used to represent aggregations over a time
     interval.
 
     """
@@ -1285,11 +1282,20 @@ def product_definition_template_15(cube, grib, full3d_cube=None):
     arrive at the given data value from the source data.
 
     """
+    # Type of spatial processing (see code table 4.15)
+    spatial_processing_code = cube.attributes['spatial_processing_type']
+
+    # Only spatial processing with no interpolation is supported at save time.
+    if spatial_processing_code != 0:
+        msg = ('Cannot save Product Definition Type 4.15 with spatial '
+               'processing type {}'.format(spatial_processing_code))
+        raise ValueError(msg)
+
     gribapi.grib_set(grib, "productDefinitionTemplateNumber", 15)
     product_definition_template_common(cube, grib, full3d_cube)
 
-    # Check that there is one and only one cell method related to the
-    # spatial processing.
+    # Check that there is one and only one cell method (statistical
+    # process) related to the spatial processing.
     if cube.cell_methods:
         spatial_cell_methods = [
             cell_method for cell_method in cube.cell_methods if 'area' in
@@ -1306,11 +1312,17 @@ def product_definition_template_15(cube, grib, full3d_cube=None):
                              "the spatial processing related cell method. "
                              "Expected ('area',), got {!r}".format(
                                  cell_method.coord_names))
+    else:
+        raise ValueError('Cannot save Product Definition Template 4.15 '
+                         'without cell method information.')
+
+    # For spatial processing code 0 there are zero points used in the
+    # interpolation.
+    number_of_points = 0
 
-        # Type of spatial processing (see code table 4.15)
-        spatial_processing = _SPATIAL_PROCESSING_TYPES.get(cell_method.method,
-                                                           255)
-        gribapi.grib_set(grib, "spatialProcessing", spatial_processing)
+    gribapi.grib_set(grib, "statisticalProcess", cell_method)
+    gribapi.grib_set(grib, "spatialProcessing", spatial_processing_code)
+    gribapi.grib_set(grib, "numberOfPointsUsed", number_of_points)
 
 
 def product_definition_template_40(cube, grib, full3d_cube=None):