Repository: NCAS-CMS/cfdm

Issue:
Incorrect `formula_terms` in output netCDF files in some cases
When a domain ancillary construct with bounds is the same netCDF variable as the parent dimension coordinate constructs, the domain ancillary construct is erroneously written to an output file in addition to the dimension coordinate construct.

```python
>>> import cfdm
>>> cfdm.environment(paths=False)
Platform: Linux-5.15.0-60-generic-x86_64-with-glibc2.35
HDF5 library: 1.12.1
netcdf library: 4.8.1
Python: 3.10.8
netCDF4: 1.6.0
numpy: 1.23.5
cfdm.core: 1.10.0.2
cftime: 1.6.2
netcdf_flattener: 1.2.0
cfdm: 1.10.0.2
```

Discussion:

Patch:
diff --git a/Changelog.rst b/Changelog.rst
index b18f7dd25..721dfade1 100644
--- a/Changelog.rst
+++ b/Changelog.rst
@@ -7,7 +7,9 @@ Version 1.10.0.2
   (https://github.com/NCAS-CMS/cfdm/issues/228)
 * Relocate the call to `NetCDFRead._customize_read_vars` to earlier in
   `NetCDFRead.read` (https://github.com/NCAS-CMS/cfdm/issues/233)
- 
+* Fixed bug that caused incorrect formula terms in output netCDF files
+  in some cases (https://github.com/NCAS-CMS/cfdm/issues/242)
+
 ----
 
 Version 1.10.0.1
diff --git a/cfdm/examplefield.py b/cfdm/examplefield.py
index 7d9acfb68..6a15fbd88 100644
--- a/cfdm/examplefield.py
+++ b/cfdm/examplefield.py
@@ -466,13 +466,19 @@ def example_field(n, _implementation=_implementation):
 
         # domain_ancillary
         c = DomainAncillary()
-        c.set_properties({"units": "m"})
+        c.set_properties(
+            {
+                "units": "m",
+                "computed_standard_name": "altitude",
+                "standard_name": "atmosphere_hybrid_height_coordinate",
+            }
+        )
         c.nc_set_variable("a")
-        data = Data([10.0], units="m", dtype="f8")
+        data = Data([20.0], units="m", dtype="f8")
         c.set_data(data)
         b = Bounds()
         b.nc_set_variable("a_bounds")
-        data = Data([[5.0, 15.0]], units="m", dtype="f8")
+        data = Data([[0, 36.6666717529297]], units="m", dtype="f8")
         b.set_data(data)
         c.set_bounds(b)
         f.set_construct(
@@ -482,11 +488,11 @@ def example_field(n, _implementation=_implementation):
         # domain_ancillary
         c = DomainAncillary()
         c.nc_set_variable("b")
-        data = Data([20.0], dtype="f8")
+        data = Data([0.997741281986237], dtype="f8")
         c.set_data(data)
         b = Bounds()
         b.nc_set_variable("b_bounds")
-        data = Data([[14.0, 26.0]], dtype="f8")
+        data = Data([[1, 0.995860934257507]], dtype="f8")
         b.set_data(data)
         c.set_bounds(b)
         f.set_construct(
@@ -946,16 +952,17 @@ def example_field(n, _implementation=_implementation):
         c = DimensionCoordinate()
         c.set_properties(
             {
+                "units": "m",
                 "computed_standard_name": "altitude",
                 "standard_name": "atmosphere_hybrid_height_coordinate",
             }
         )
         c.nc_set_variable("atmosphere_hybrid_height_coordinate")
-        data = Data([1.5], dtype="f8")
+        data = Data([20.0], dtype="f8")
         c.set_data(data)
         b = Bounds()
         b.nc_set_variable("atmosphere_hybrid_height_coordinate_bounds")
-        data = Data([[1.0, 2.0]], dtype="f8")
+        data = Data([[0, 36.6666717529297]], dtype="f8")
         b.set_data(data)
         c.set_bounds(b)
         f.set_construct(
diff --git a/cfdm/read_write/netcdf/netcdfread.py b/cfdm/read_write/netcdf/netcdfread.py
index 201fe8b7e..29aa0ef19 100644
--- a/cfdm/read_write/netcdf/netcdfread.py
+++ b/cfdm/read_write/netcdf/netcdfread.py
@@ -4569,6 +4569,9 @@ class NetCDFRead(IORead):
         elif nodes:
             attribute = "nodes"
 
+        # Make sure that the bounds attribute is removed
+        properties.pop(attribute, None)
+
         if dimension:
             properties.pop("compress", None)
             c = self.implementation.initialise_DimensionCoordinate()
@@ -4778,7 +4781,6 @@ class NetCDFRead(IORead):
 
         # Store the netCDF variable name
         self.implementation.nc_set_variable(cell_measure, ncvar)
-
         if ncvar in g["external_variables"]:
             # The cell measure variable is in an unknown external file
             self.implementation.nc_set_external(construct=cell_measure)


Test patch:
diff --git a/cfdm/test/test_Constructs.py b/cfdm/test/test_Constructs.py
index 6cae2d2a4..350e8d987 100644
--- a/cfdm/test/test_Constructs.py
+++ b/cfdm/test/test_Constructs.py
@@ -390,7 +390,7 @@ class ConstructsTest(unittest.TestCase):
                 e = d.inverse_filter()
                 self.assertEqual(len(e), len(c) - len(d))
 
-        self.assertEqual(len(c.filter_by_property(standard_name=None)), 8)
+        self.assertEqual(len(c.filter_by_property(standard_name=None)), 9)
 
         with self.assertRaises(ValueError):
             c.filter_by_property("too many", "modes")
diff --git a/cfdm/test/test_read_write.py b/cfdm/test/test_read_write.py
index 1a2009ff4..0bc53bae2 100644
--- a/cfdm/test/test_read_write.py
+++ b/cfdm/test/test_read_write.py
@@ -7,6 +7,7 @@ import subprocess
 import tempfile
 import unittest
 
+import netCDF4
 import numpy as np
 
 faulthandler.enable()  # to debug seg faults and timeouts
@@ -938,6 +939,51 @@ class read_writeTest(unittest.TestCase):
         self.assertFalse(g.array.count())
         self.assertTrue(g.construct("grid_latitude").array.count())
 
+    def test_read_write_domain_ancillary(self):
+        """Test when domain ancillary equals dimension coordinate."""
+        f = cfdm.example_field(1)
+
+        # Check the domain ancillary does indeed equal the dimension
+        # coordinate
+        self.assertTrue(
+            f.domain_ancillary("atmosphere_hybrid_height_coordinate").equals(
+                f.dimension_coordinate("atmosphere_hybrid_height_coordinate"),
+                ignore_type=True,
+            )
+        )
+
+        cfdm.write(f, tmpfile)
+        g = cfdm.read(tmpfile)
+        self.assertEqual(len(g), 1)
+        g = g[0]
+        self.assertTrue(f.equals(g))
+
+        nc = netCDF4.Dataset(tmpfile, "r")
+
+        z = nc.variables["atmosphere_hybrid_height_coordinate"]
+        # In the following test we are checking that it is not 'a: a
+        # b: b orog: surface_altitude'
+        self.assertEqual(
+            z.getncattr("formula_terms"),
+            "a: atmosphere_hybrid_height_coordinate b: b orog: surface_altitude",
+        )
+
+        zb = nc.variables["atmosphere_hybrid_height_coordinate_bounds"]
+        # In the following test we are checking that it is not 'a:
+        # a_bounds b: b orog: surface_altitude'
+        self.assertEqual(
+            zb.getncattr("formula_terms"),
+            "a: atmosphere_hybrid_height_coordinate_bounds b: b_bounds orog: surface_altitude",
+        )
+
+        nc.close()
+
+        # Test the full round trip
+        cfdm.write(g, tmpfile0)
+        h = cfdm.read(tmpfile0)
+        self.assertEqual(len(h), 1)
+        self.assertTrue(f.equals(h[0]))
+
 
 if __name__ == "__main__":
     print("Run date:", datetime.datetime.now())