Merge launch_ancils feature branch

docs/iris/src/whatsnew/contributions_3.0.0/bugfix_2019-Nov-21_cell_measure_attributes.txt

@@ -0,0 +1,2 @@
+* Fixed a bug where the attributes of cell measures in netcdf-CF files were discarded on
+  loading.  They now appear on the CellMeasure in the loaded cube.

docs/iris/src/whatsnew/contributions_3.0.0/newfeature_2019-Nov-21_netcdf_ancillary_data.txt

@@ -0,0 +1,2 @@
+* CF Ancillary Variables are now loaded from and saved to netcdf-CF files.
+

lib/iris/fileformats/netcdf.py

@@ -460,7 +460,10 @@ def __setstate__(self, state):
 
 def _assert_case_specific_facts(engine, cf, cf_group):
     # Initialise pyke engine "provides" hooks.
-    engine.provides["coordinates"] = []
+    # These are used to patch non-processed element attributes after rules activation.
+    engine.cube_parts["coordinates"] = []
+    engine.cube_parts["cell_measures"] = []
+    engine.cube_parts["ancillary_variables"] = []
 
     # Assert facts for CF coordinates.
     for cf_name in cf_group.coordinates.keys():
@@ -480,6 +483,12 @@ def _assert_case_specific_facts(engine, cf, cf_group):
             _PYKE_FACT_BASE, "cell_measure", (cf_name,)
         )
 
+    # Assert facts for CF ancillary variables.
+    for cf_name in cf_group.ancillary_variables.keys():
+        engine.add_case_specific_fact(
+            _PYKE_FACT_BASE, "ancillary_variable", (cf_name,)
+        )
+
     # Assert facts for CF grid_mappings.
     for cf_name in cf_group.grid_mappings.keys():
         engine.add_case_specific_fact(
@@ -587,7 +596,7 @@ def _load_cube(engine, cf, cf_var, filename):
     # Initialise pyke engine rule processing hooks.
     engine.cf_var = cf_var
     engine.cube = cube
-    engine.provides = {}
+    engine.cube_parts = {}
     engine.requires = {}
     engine.rule_triggered = set()
     engine.filename = filename
@@ -598,31 +607,38 @@ def _load_cube(engine, cf, cf_var, filename):
     # Run pyke inference engine with forward chaining rules.
     engine.activate(_PYKE_RULE_BASE)
 
-    # Populate coordinate attributes with the untouched attributes from the
-    # associated CF-netCDF variable.
-    coordinates = engine.provides.get("coordinates", [])
-
+    # Having run the rules, now populate the attributes of all the cf elements with the
+    # "unused" attributes from the associated CF-netCDF variable.
+    # That is, all those that aren't CF reserved terms.
     def attribute_predicate(item):
         return item[0] not in _CF_ATTRS
 
-    for coord, cf_var_name in coordinates:
-        tmpvar = filter(
-            attribute_predicate, cf.cf_group[cf_var_name].cf_attrs_unused()
-        )
+    def add_unused_attributes(iris_object, cf_var):
+        tmpvar = filter(attribute_predicate, cf_var.cf_attrs_unused())
         for attr_name, attr_value in tmpvar:
-            _set_attributes(coord.attributes, attr_name, attr_value)
+            _set_attributes(iris_object.attributes, attr_name, attr_value)
+
+    def fix_attributes_all_elements(role_name):
+        elements_and_names = engine.cube_parts.get(role_name, [])
+
+        for iris_object, cf_var_name in elements_and_names:
+            add_unused_attributes(iris_object, cf.cf_group[cf_var_name])
+
+    # Populate the attributes of all coordinates, cell-measures and ancillary-vars.
+    fix_attributes_all_elements("coordinates")
+    fix_attributes_all_elements("ancillary_variables")
+    fix_attributes_all_elements("cell_measures")
 
-    tmpvar = filter(attribute_predicate, cf_var.cf_attrs_unused())
-    # Attach untouched attributes of the associated CF-netCDF data variable to
-    # the cube.
-    for attr_name, attr_value in tmpvar:
-        _set_attributes(cube.attributes, attr_name, attr_value)
+    # Also populate attributes of the top-level cube itself.
+    add_unused_attributes(cube, cf_var)
 
+    # Work out reference names for all the coords.
     names = {
         coord.var_name: coord.standard_name or coord.var_name or "unknown"
         for coord in cube.coords()
     }
 
+    # Add all the cube cell methods.
     cube.cell_methods = [
         iris.coords.CellMethod(
             method=method.method,
@@ -662,7 +678,7 @@ def coord_from_term(term):
             # Convert term names to coordinates (via netCDF variable names).
             name = engine.requires["formula_terms"].get(term, None)
             if name is not None:
-                for coord, cf_var_name in engine.provides["coordinates"]:
+                for coord, cf_var_name in engine.cube_parts["coordinates"]:
                     if cf_var_name == name:
                         return coord
                 warnings.warn(

lib/iris/tests/results/netcdf/TestNetCDFSave__ancillaries/flag.cdl

@@ -0,0 +1,46 @@
+dimensions:
+	grid_latitude = 9 ;
+	grid_longitude = 11 ;
+	time = 7 ;
+variables:
+	int64 air_potential_temperature(time, grid_latitude, grid_longitude) ;
+		air_potential_temperature:standard_name = "air_potential_temperature" ;
+		air_potential_temperature:units = "K" ;
+		air_potential_temperature:grid_mapping = "rotated_latitude_longitude" ;
+		air_potential_temperature:coordinates = "air_pressure forecast_period" ;
+		air_potential_temperature:ancillary_variables = "quality_flag" ;
+	int rotated_latitude_longitude ;
+		rotated_latitude_longitude:grid_mapping_name = "rotated_latitude_longitude" ;
+		rotated_latitude_longitude:longitude_of_prime_meridian = 0. ;
+		rotated_latitude_longitude:earth_radius = 6371229. ;
+		rotated_latitude_longitude:grid_north_pole_latitude = 37.5 ;
+		rotated_latitude_longitude:grid_north_pole_longitude = 177.5 ;
+		rotated_latitude_longitude:north_pole_grid_longitude = 0. ;
+	double time(time) ;
+		time:axis = "T" ;
+		time:units = "hours since 1970-01-01 00:00:00" ;
+		time:standard_name = "time" ;
+		time:calendar = "gregorian" ;
+	double grid_latitude(grid_latitude) ;
+		grid_latitude:axis = "Y" ;
+		grid_latitude:units = "degrees" ;
+		grid_latitude:standard_name = "grid_latitude" ;
+	double grid_longitude(grid_longitude) ;
+		grid_longitude:axis = "X" ;
+		grid_longitude:units = "degrees" ;
+		grid_longitude:standard_name = "grid_longitude" ;
+	double air_pressure ;
+		air_pressure:units = "Pa" ;
+		air_pressure:standard_name = "air_pressure" ;
+	double forecast_period(time) ;
+		forecast_period:units = "hours" ;
+		forecast_period:standard_name = "forecast_period" ;
+	byte quality_flag(time, grid_latitude, grid_longitude) ;
+		quality_flag:long_name = "quality_flag" ;
+		quality_flag:flag_meanings = "PASS FAIL MISSING" ;
+		quality_flag:flag_values = 1b, 2b, 9b ;
+
+// global attributes:
+		:source = "Iris test case" ;
+		:Conventions = "CF-1.7" ;
+}

lib/iris/tests/test_netcdf.py

@@ -32,13 +32,21 @@
 from iris.fileformats.netcdf import load_cubes as nc_load_cubes
 import iris.std_names
 import iris.util
+from iris.coords import AncillaryVariable, CellMeasure
 import iris.coord_systems as icoord_systems
 import iris.tests.stock as stock
 from iris._lazy_data import is_lazy_data
 
 
 @tests.skip_data
 class TestNetCDFLoad(tests.IrisTest):
+    def setUp(self):
+        self.tmpdir = None
+
+    def tearDown(self):
+        if self.tmpdir is not None:
+            shutil.rmtree(self.tmpdir)
+
     def test_monotonic(self):
         cubes = iris.load(
             tests.get_data_path(
@@ -243,6 +251,153 @@ def test_cell_methods(self):
 
         self.assertCML(cubes, ("netcdf", "netcdf_cell_methods.cml"))
 
+    def test_ancillary_variables(self):
+        # Note: using a CDL string as a test data reference, rather than a binary file.
+        ref_cdl = """
+            netcdf cm_attr {
+            dimensions:
+                axv = 3 ;
+            variables:
+                int64 qqv(axv) ;
+                    qqv:long_name = "qq" ;
+                    qqv:units = "1" ;
+                    qqv:ancillary_variables = "my_av" ;
+                int64 axv(axv) ;
+                    axv:units = "1" ;
+                    axv:long_name = "x" ;
+                double my_av(axv) ;
+                    my_av:units = "1" ;
+                    my_av:long_name = "refs" ;
+                    my_av:custom = "extra-attribute";
+            data:
+                axv = 1, 2, 3;
+                my_av = 11., 12., 13.;
+            }
+            """
+        self.tmpdir = tempfile.mkdtemp()
+        cdl_path = os.path.join(self.tmpdir, "tst.cdl")
+        nc_path = os.path.join(self.tmpdir, "tst.nc")
+        # Write CDL string into a temporary CDL file.
+        with open(cdl_path, "w") as f_out:
+            f_out.write(ref_cdl)
+        # Use ncgen to convert this into an actual (temporary) netCDF file.
+        command = "ncgen -o {} {}".format(nc_path, cdl_path)
+        check_call(command, shell=True)
+        # Load with iris.fileformats.netcdf.load_cubes, and check expected content.
+        cubes = list(nc_load_cubes(nc_path))
+        self.assertEqual(len(cubes), 1)
+        avs = cubes[0].ancillary_variables()
+        self.assertEqual(len(avs), 1)
+        expected = AncillaryVariable(
+            np.ma.array([11.0, 12.0, 13.0]),
+            long_name="refs",
+            var_name="my_av",
+            units="1",
+            attributes={"custom": "extra-attribute"},
+        )
+        self.assertEqual(avs[0], expected)
+
+    def test_status_flags(self):
+        # Note: using a CDL string as a test data reference, rather than a binary file.
+        ref_cdl = """
+            netcdf cm_attr {
+            dimensions:
+                axv = 3 ;
+            variables:
+                int64 qqv(axv) ;
+                    qqv:long_name = "qq" ;
+                    qqv:units = "1" ;
+                    qqv:ancillary_variables = "my_av" ;
+                int64 axv(axv) ;
+                    axv:units = "1" ;
+                    axv:long_name = "x" ;
+                byte my_av(axv) ;
+                    my_av:long_name = "qq status_flag" ;
+                    my_av:flag_values = 1b, 2b ;
+                    my_av:flag_meanings = "a b" ;
+            data:
+                axv = 11, 21, 31;
+                my_av = 1b, 1b, 2b;
+            }
+            """
+        self.tmpdir = tempfile.mkdtemp()
+        cdl_path = os.path.join(self.tmpdir, "tst.cdl")
+        nc_path = os.path.join(self.tmpdir, "tst.nc")
+        # Write CDL string into a temporary CDL file.
+        with open(cdl_path, "w") as f_out:
+            f_out.write(ref_cdl)
+        # Use ncgen to convert this into an actual (temporary) netCDF file.
+        command = "ncgen -o {} {}".format(nc_path, cdl_path)
+        check_call(command, shell=True)
+        # Load with iris.fileformats.netcdf.load_cubes, and check expected content.
+        cubes = list(nc_load_cubes(nc_path))
+        self.assertEqual(len(cubes), 1)
+        avs = cubes[0].ancillary_variables()
+        self.assertEqual(len(avs), 1)
+        expected = AncillaryVariable(
+            np.ma.array([1, 1, 2], dtype=np.int8),
+            long_name="qq status_flag",
+            var_name="my_av",
+            units="no_unit",
+            attributes={
+                "flag_values": np.array([1, 2], dtype=np.int8),
+                "flag_meanings": "a b",
+            },
+        )
+        self.assertEqual(avs[0], expected)
+
+    def test_cell_measures(self):
+        # Note: using a CDL string as a test data reference, rather than a binary file.
+        ref_cdl = """
+            netcdf cm_attr {
+            dimensions:
+                axv = 3 ;
+                ayv = 2 ;
+            variables:
+                int64 qqv(ayv, axv) ;
+                    qqv:long_name = "qq" ;
+                    qqv:units = "1" ;
+                    qqv:cell_measures = "area: my_areas" ;
+                int64 ayv(ayv) ;
+                    ayv:units = "1" ;
+                    ayv:long_name = "y" ;
+                int64 axv(axv) ;
+                    axv:units = "1" ;
+                    axv:long_name = "x" ;
+                double my_areas(ayv, axv) ;
+                    my_areas:units = "m2" ;
+                    my_areas:long_name = "standardised cell areas" ;
+                    my_areas:custom = "extra-attribute";
+            data:
+                axv = 11, 12, 13;
+                ayv = 21, 22;
+                my_areas = 110., 120., 130., 221., 231., 241.;
+            }
+            """
+        self.tmpdir = tempfile.mkdtemp()
+        cdl_path = os.path.join(self.tmpdir, "tst.cdl")
+        nc_path = os.path.join(self.tmpdir, "tst.nc")
+        # Write CDL string into a temporary CDL file.
+        with open(cdl_path, "w") as f_out:
+            f_out.write(ref_cdl)
+        # Use ncgen to convert this into an actual (temporary) netCDF file.
+        command = "ncgen -o {} {}".format(nc_path, cdl_path)
+        check_call(command, shell=True)
+        # Load with iris.fileformats.netcdf.load_cubes, and check expected content.
+        cubes = list(nc_load_cubes(nc_path))
+        self.assertEqual(len(cubes), 1)
+        cms = cubes[0].cell_measures()
+        self.assertEqual(len(cms), 1)
+        expected = CellMeasure(
+            np.ma.array([[110.0, 120.0, 130.0], [221.0, 231.0, 241.0]]),
+            measure="area",
+            var_name="my_areas",
+            long_name="standardised cell areas",
+            units="m2",
+            attributes={"custom": "extra-attribute"},
+        )
+        self.assertEqual(cms[0], expected)
+
     def test_deferred_loading(self):
         # Test exercising CF-netCDF deferred loading and deferred slicing.
         # shape (31, 161, 320)
@@ -305,14 +460,21 @@ def test_default_units(self):
             variables:
                 int64 qqv(ayv, axv) ;
                     qqv:long_name = "qq" ;
+                    qqv:ancillary_variables = "my_av" ;
+                    qqv:cell_measures = "area: my_areas" ;
                 int64 ayv(ayv) ;
                     ayv:long_name = "y" ;
                 int64 axv(axv) ;
                     axv:units = "1" ;
                     axv:long_name = "x" ;
+                double my_av(axv) ;
+                    my_av:long_name = "refs" ;
+                double my_areas(ayv, axv) ;
+                    my_areas:long_name = "areas" ;
             data:
                 axv = 11, 12, 13;
                 ayv = 21, 22;
+                my_areas = 110., 120., 130., 221., 231., 241.;
             }
             """
         self.tmpdir = tempfile.mkdtemp()
@@ -330,6 +492,12 @@ def test_default_units(self):
         self.assertEqual(cubes[0].units, as_unit("unknown"))
         self.assertEqual(cubes[0].coord("y").units, as_unit("unknown"))
         self.assertEqual(cubes[0].coord("x").units, as_unit(1))
+        self.assertEqual(
+            cubes[0].ancillary_variable("refs").units, as_unit("unknown")
+        )
+        self.assertEqual(
+            cubes[0].cell_measure("areas").units, as_unit("unknown")
+        )
 
     def test_units(self):
         # Test exercising graceful cube and coordinate units loading.
@@ -1108,6 +1276,21 @@ def test_aliases(self):
             iris.save([testcube_1, testcube_2], filename)
             self.assertCDL(filename)
 
+    def test_flag(self):
+        testcube = stock.realistic_3d()
+        flag = iris.coords.AncillaryVariable(
+            np.ones(testcube.shape, dtype=np.int8),
+            long_name="quality_flag",
+            attributes={
+                "flag_meanings": "PASS FAIL MISSING",
+                "flag_values": np.array([1, 2, 9], dtype=np.int8),
+            },
+        )
+        testcube.add_ancillary_variable(flag, (0, 1, 2))
+        with self.temp_filename(suffix=".nc") as filename:
+            iris.save(testcube, filename)
+            self.assertCDL(filename)
+
 
 class TestNetCDF3SaveInteger(tests.IrisTest):
     def setUp(self):