Revert pv level indexing change (#165)

lbdreyer · pp-mo · commit 2cf099a7b90d · 2019-12-02T10:46:46.000Z
diff --git a/iris_grib/_load_convert.py b/iris_grib/_load_convert.py
@@ -1508,14 +1508,17 @@ def hybrid_factories(section, metadata):
                                 {'long_name': 'sigma'},
                                 Reference('ref_surface_pressure')]
 
-            # Create the level pressure scalar coordinate.
+            # Create the level height/pressure scalar coordinate.
+            # scaledValue represents the level number, which is used to select
+            # the sigma and delta values as follows:
+            # sigma, delta = PV[i], PV[NV/2+i] : where i=1..level_number
             pv = section['pv']
-            offset = scaledValue - 1
+            offset = scaledValue
             coord = DimCoord(pv[offset], long_name=level_value_name,
                              units=level_value_units)
             metadata['aux_coords_and_dims'].append((coord, None))
             # Create the sigma scalar coordinate.
-            offset = NV // 2 + scaledValue - 1
+            offset = NV // 2 + scaledValue
             coord = AuxCoord(pv[offset], long_name='sigma')
             metadata['aux_coords_and_dims'].append((coord, None))
             # Create the associated factory reference.
diff --git a/iris_grib/_save_rules.py b/iris_grib/_save_rules.py
@@ -950,10 +950,9 @@ def set_fixed_surfaces(cube, grib, full3d_cube=None):
                    'maximum value = {}.')
             raise ValueError(msg.format(max_valid_nlevels, n_levels))
         # In sample data we have seen, there seems to be an extra missing data
-        # value after each set of n-levels coefficients.
-        # For now, we retain that encoding. This should not cause problems as
-        # values are indexed according to model-level-number,
-        # I.E. sigma, delta = PV[i], PV[NV/2+i]
+        # value *before* each set of n-levels coefficients.
+        # Note: values are indexed according to model_level_number,
+        # I.E. sigma, delta = PV[i], PV[NV/2+i] : where i=1..n_levels
         n_coeffs = n_levels + 1
         coeffs_array = np.zeros(n_coeffs * 2, dtype=np.float32)
         for n_lev, height, sigma in zip(model_levels,
@@ -962,11 +961,11 @@ def set_fixed_surfaces(cube, grib, full3d_cube=None):
             # Record all the level coefficients coming from the 'full' cube.
             # Note: if some model levels are missing, we must still have the
             # coeffs at the correct index according to the model_level_number
-            # value, i.e. at [level - 1] and [NV // 2 + level - 1].
+            # value, i.e. at [level] and [NV // 2 + level].
             # Thus, we can *not* paste the values in a block: each one needs to
             # go in the index corresponding to its 'model_level_number' value.
-            coeffs_array[n_lev - 1] = height
-            coeffs_array[n_coeffs + n_lev - 1] = sigma
+            coeffs_array[n_lev] = height
+            coeffs_array[n_coeffs + n_lev] = sigma
         pv_values = [float(el) for el in coeffs_array]
         # eccodes does not support writing numpy.int64, cast to python int
         gribapi.grib_set(grib, "NV", int(n_coeffs * 2))
diff --git a/iris_grib/tests/integration/load_convert/test_load_hybrid_coords.py b/iris_grib/tests/integration/load_convert/test_load_hybrid_coords.py
@@ -35,14 +35,13 @@ def test_hybrid_height_coords_values(self):
                               [1, 11, 21])
 
         # check sigma values correctly loaded from indices 1, 11, 21.
-        # NOTE: level[0] == 1, so sigma[0] == 1.0 :  This makes sense !
         self.assertArrayAllClose(cube.coord('sigma').points,
-                                 [1.0, 0.911, 0.694],
+                                 [0.998, 0.894, 0.667],
                                  atol=0.001)
 
         # check height values too.
         self.assertArrayAllClose(cube.coord('level_height').points,
-                                 [0., 800.,  2933.],
+                                 [20., 953.3,  3220.],
                                  atol=0.5)
 
 
@@ -63,9 +62,9 @@ def test_hybrid_pressure_coords_values(self):
         self.assertArrayEqual(cube.coord('model_level_number').points,
                               [1, 51, 91])
         self.assertArrayAllClose(cube.coord('sigma').points,
-                                 [0., 0.036, 0.998], atol=0.001)
+                                 [0., 0.045, 1.], atol=0.001)
         self.assertArrayAllClose(cube.coord('level_pressure').points,
-                                 [0., 18191.0, 0.00316], rtol=0.001)
+                                 [2.00004, 18716.9688, 0.], rtol=0.0001)
         self.assertArrayAllClose(
             cube.coord('surface_air_pressure')[:2, :3].points,
             [[103493.8, 103493.8, 103493.8],
@@ -74,9 +73,9 @@ def test_hybrid_pressure_coords_values(self):
         # Also check a few values from the derived coord.
         self.assertArrayAllClose(
             cube.coord('air_pressure')[:, :3, 0].points,
-            [[0., 0., 0.],
-             [21940.3, 21936.9, 21932.8],
-             [103248.5, 103156.0, 103041.0]], atol=0.1)
+            [[2., 2., 2.],
+             [23389.3, 23385.1, 23379.9],
+             [103493.8, 103401.0, 103285.8]], atol=0.1)
 
 
 if __name__ == '__main__':
diff --git a/iris_grib/tests/integration/save_rules/test_save_hybrid_coords.py b/iris_grib/tests/integration/save_rules/test_save_hybrid_coords.py
@@ -57,8 +57,7 @@ def test_save(self):
             # Note: gaps here are because we took model levels = (1, 3, 5).
             self.assertArrayAllClose(
                 msgs[0].sections[4]['pv'],
-                [5., 0, 45., 0, 111.667, 0,
-                 0.999, 0, 0.995, 0, 0.987, 0],
+                [0, 5., 0, 45., 0, 111.667, 0, 0.999, 0, 0.995, 0, 0.987],
                 atol=0.0015)
 
             # Check message #2-of-3 has the correctly encoded hybrid height.
@@ -114,8 +113,8 @@ def test_save(self):
             # Note: HUGE gaps here because we took model levels = (1, 51, 91).
             self.assertEqual(msgs[0].sections[4]['NV'], 184)
             pv_expected = np.zeros(184, dtype=np.float64)
-            pv_expected[[0, 50, 90]] = [0., 18191.03, 0.003]
-            pv_expected[[92, 142, 182]] = [0., 0.036, 0.998]
+            pv_expected[[1, 51, 91]] = [0., 18191.03, 0.003]
+            pv_expected[[93, 143, 183]] = [0., 0.036, 0.998]
             self.assertArrayAllClose(
                 msgs[0].sections[4]['pv'], pv_expected, atol=0.001)
 
diff --git a/iris_grib/tests/results/unit/load_cubes/load_cubes/reduced_raw.cml b/iris_grib/tests/results/unit/load_cubes/load_cubes/reduced_raw.cml
@@ -18,7 +18,7 @@
         </auxCoord>
       </coord>
       <coord>
-        <dimCoord id="c9091dc6" long_name="level_pressure" points="[0.0]" shape="(1,)" units="Unit('Pa')" value_type="float64"/>
+        <dimCoord id="c9091dc6" long_name="level_pressure" points="[2.00004005432]" shape="(1,)" units="Unit('Pa')" value_type="float64"/>
       </coord>
       <coord datadims="[0]">
         <auxCoord id="f913a8b3" points="[0.0, 18.0, 36.0, ..., 306.0, 324.0, 342.0]" shape="(13280,)" standard_name="longitude" units="Unit('degrees')" value_type="float64">
