Repository: NCAS-CMS/cf-python

Issue:
Allow `relaxed_units` in `cf.aggregate` to cope with invalid (in addition to missing) units
Currently, the  `relaxed_units` keyword to `cf.aggregate` copes with missing units (i.e. `cf.units()`), but not invalid but otherwise equal units:
```python
>>> # 3.14.0 behaviour
>>> f = cf.example_field(0)
>>> bad_units = cf.Units("bad-units")
>>  f.override_units(bad_units, inplace=True)
>>> g = f[:2]
>>> h = f[2:]
>>> i = cf.aggregate([g, h], relaxed_units=True)
[<CF Field: specific_humidity(latitude(2), longitude(8)) bad-units>,
 <CF Field: specific_humidity(latitude(3), longitude(8)) bad-units>]
```
It would be good to extend the functionality to include this case:
```python
>>> # Desired new behaviour
>>> i = cf.aggregate([g, h], relaxed_units=True)
[<CF Field: specific_humidity(latitude(5), longitude(8)) bad-units>]
```  

Discussion:

Patch:
diff --git a/Changelog.rst b/Changelog.rst
index f6345b30b..aae2752d8 100644
--- a/Changelog.rst
+++ b/Changelog.rst
@@ -3,6 +3,8 @@ version 3.14.1
 
 **2023-??-??**
 
+* New keyword parameter to `cf.Data.concatenate`: ``relaxed_units``
+  (https://github.com/NCAS-CMS/cf-python/issues/602)
 * New method: `cf.Field.del_properties`
   (https://github.com/NCAS-CMS/cf-python/issues/598)
 * New keyword parameter to `cf.unique_constructs`:
diff --git a/cf/aggregate.py b/cf/aggregate.py
index 94abada77..0e9dd1a1f 100644
--- a/cf/aggregate.py
+++ b/cf/aggregate.py
@@ -177,9 +177,10 @@ class _Meta:
             relaxed_units: `bool`, optional
                 If True then assume that field and metadata constructs
                 with the same identity but missing units actually have
-                equivalent (but unspecified) units, so that aggregation
-                may occur. By default such field constructs are not
-                aggregatable.
+                equivalent (but unspecified) units, so that
+                aggregation may occur. Also assumes that invalid but
+                otherwise equal units are equal. By default such field
+                constructs are not aggregatable.
 
             allow_no_identity: `bool`, optional
                 If True then assume that field and metadata constructs
@@ -237,6 +238,7 @@ class _Meta:
         )
 
         self.relaxed_identities = relaxed_identities
+        self.relaxed_units = relaxed_units
         self.strict_identities = strict_identities
         self.field_identity = field_identity
         self.ncvar_identities = ncvar_identities
@@ -854,21 +856,32 @@ class _Meta:
         _canonical_units = self._canonical_units
 
         if identity in _canonical_units:
-            if var_units:
+            if var_units.isvalid:
+                if var_units:
+                    for u in _canonical_units[identity]:
+                        if var_units.equivalent(u):
+                            return u
+
+                    # Still here?
+                    _canonical_units[identity].append(var_units)
+                elif relaxed_units or variable.dtype.kind in ("S", "U"):
+                    return _no_units
+
+            elif relaxed_units:
                 for u in _canonical_units[identity]:
-                    if var_units.equivalent(u):
+                    if u.isvalid:
+                        continue
+
+                    if var_units.__dict__ == u.__dict__:
                         return u
 
                 # Still here?
                 _canonical_units[identity].append(var_units)
-
-            elif relaxed_units or variable.dtype.kind in ("S", "U"):
-                var_units = _no_units
         else:
-            if var_units:
+            if var_units or (relaxed_units and not var_units.isvalid):
                 _canonical_units[identity] = [var_units]
             elif relaxed_units or variable.dtype.kind in ("S", "U"):
-                var_units = _no_units
+                return _no_units
 
         # Still here?
         return var_units
@@ -1196,7 +1209,11 @@ class _Meta:
         # * whether or not there is a data array
         Type = f.construct_type
         Identity = self.identity
-        Units = self.units.formatted(definition=True)
+        if self.units.isvalid:
+            Units = self.units.formatted(definition=True)
+        else:
+            Units = self.units.units
+
         Cell_methods = self.cell_methods
         Data = self.has_data
         #        signature_append = signature.append
@@ -1501,7 +1518,8 @@ def aggregate(
             If True then assume that field and metadata constructs
             with the same identity but missing units actually have
             equivalent (but unspecified) units, so that aggregation
-            may occur. By default such field constructs are not
+            may occur. Also assumes that invalid but otherwise equal
+            units are equal. By default such field constructs are not
             aggregatable.
 
         allow_no_identity: `bool`, optional
@@ -1883,7 +1901,7 @@ def aggregate(
 
             continue
 
-        if not meta[0].units.isvalid:
+        if not relaxed_units and not meta[0].units.isvalid:
             if is_log_level_info(logger):
                 x = ", ".join(set(repr(m.units) for m in meta))
                 logger.info(
@@ -2046,6 +2064,7 @@ def aggregate(
                         atol=atol,
                         verbose=verbose,
                         concatenate=concatenate,
+                        relaxed_units=relaxed_units,
                         copy=(copy or not exclude),
                     )
 
@@ -2945,7 +2964,14 @@ def _ok_coordinate_arrays(meta, axis, overlap, contiguous, verbose=None):
 
 @_manage_log_level_via_verbosity
 def _aggregate_2_fields(
-    m0, m1, rtol=None, atol=None, verbose=None, concatenate=True, copy=True
+    m0,
+    m1,
+    rtol=None,
+    atol=None,
+    verbose=None,
+    concatenate=True,
+    relaxed_units=False,
+    copy=True,
 ):
     """Aggregate two fields, returning the _Meta object of the
     aggregated field.
@@ -3132,6 +3158,7 @@ def _aggregate_2_fields(
             data = Data.concatenate(
                 (construct0.get_data(), construct1.get_data()),
                 axis,
+                relaxed_units=relaxed_units,
             )
             construct0.set_data(data, copy=False)
             if construct0.has_bounds():
@@ -3141,6 +3168,7 @@ def _aggregate_2_fields(
                         construct1.bounds.get_data(_fill_value=False),
                     ),
                     axis,
+                    relaxed_units=relaxed_units,
                 )
                 construct0.bounds.set_data(data, copy=False)
         else:
@@ -3151,6 +3179,7 @@ def _aggregate_2_fields(
                     construct0.get_data(_fill_value=False),
                 ),
                 axis,
+                relaxed_units=relaxed_units,
             )
             construct0.set_data(data)
             if construct0.has_bounds():
@@ -3160,6 +3189,7 @@ def _aggregate_2_fields(
                         construct0.bounds.get_data(_fill_value=False),
                     ),
                     axis,
+                    relaxed_units=relaxed_units,
                 )
                 construct0.bounds.set_data(data)
 
@@ -3210,12 +3240,16 @@ def _aggregate_2_fields(
         if direction0:
             # The fields are increasing along the aggregating axis
             data = Data.concatenate(
-                (parent0.get_data(), parent1.get_data()), axis
+                (parent0.get_data(), parent1.get_data()),
+                axis,
+                relaxed_units=relaxed_units,
             )
         else:
             # The fields are decreasing along the aggregating axis
             data = Data.concatenate(
-                (parent1.get_data(), parent0.get_data()), axis
+                (parent1.get_data(), parent0.get_data()),
+                axis,
+                relaxed_units=relaxed_units,
             )
 
         # Update the size of the aggregating axis in parent0
diff --git a/cf/data/data.py b/cf/data/data.py
index 2c5442d79..3b730176d 100644
--- a/cf/data/data.py
+++ b/cf/data/data.py
@@ -3532,7 +3532,7 @@ class Data(DataClassDeprecationsMixin, Container, cfdm.Data):
         return d
 
     @classmethod
-    def concatenate(cls, data, axis=0, cull_graph=True):
+    def concatenate(cls, data, axis=0, cull_graph=True, relaxed_units=False):
         """Join a sequence of data arrays together.
 
         .. seealso:: `cull_graph`
@@ -3556,6 +3556,27 @@ class Data(DataClassDeprecationsMixin, Container, cfdm.Data):
 
             {{cull_graph: `bool`, optional}}
 
+            relaxed_units: `bool`, optional
+                If True then allow the concatenation of data arrays
+                with invalid but otherwise equal units. By default, if
+                any data array has invalid units then the
+                concatenation will fail.
+
+                A `Units` object is considered to be invalid if its
+                `!isvalid` attribute is `False`:
+
+                >>> d = cf.Data(9, 'metre')
+                >>> d.Units.isvalid
+                True
+                >>> d = cf.Data(9)
+                >>> d.Units.isvalid
+                True
+                >>> d = cf.Data(9, 'bad-units')
+                >>> d.Units.isvalid
+                False
+
+                .. versionadded:: 3.14.1
+
         :Returns:
 
             `Data`
@@ -3618,12 +3639,21 @@ class Data(DataClassDeprecationsMixin, Container, cfdm.Data):
                 data1.insert_dimension(inplace=True)
 
             # Check and conform, if necessary, the units of all inputs
-            if not units0.equivalent(data1.Units):
+            units1 = data1.Units
+            if (
+                relaxed_units
+                and not units0.isvalid
+                and not units1.isvalid
+                and units0.__dict__ == units1.__dict__
+            ):
+                # Allow identical invalid units to be equal
+                pass
+            elif not units0.equivalent(units1):
                 raise ValueError(
                     "Can't concatenate: All the input arrays must have "
                     "equivalent units"
                 )
-            elif not units0.equals(data1.Units):  # conform for consistency
+            elif not units0.equals(units1):
                 if not copied:
                     data1 = data1.copy()
 
@@ -10321,7 +10351,7 @@ class Data(DataClassDeprecationsMixin, Container, cfdm.Data):
         **Performance**
 
         An unnecessary task is one which does not contribute to the
-        computed result. Such tasks will be automatically removed
+        computed result. Such tasks are always automatically removed
         (culled) at compute time, but removing them beforehand might
         improve performance by reducing the amount of work done in
         later steps.


Test patch:
diff --git a/cf/test/test_Data.py b/cf/test/test_Data.py
index 305399857..cdc72acd7 100644
--- a/cf/test/test_Data.py
+++ b/cf/test/test_Data.py
@@ -1173,6 +1173,35 @@ class DataTest(unittest.TestCase):
         self.assertEqual(f.shape, f_np.shape)
         self.assertTrue((f.array == f_np).all())
 
+        # Check concatenation with one invalid units
+        d.override_units(cf.Units("foo"), inplace=1)
+        with self.assertRaises(ValueError):
+            f = cf.Data.concatenate([d, e], axis=1, relaxed_units=True)
+
+        with self.assertRaises(ValueError):
+            f = cf.Data.concatenate([d, e], axis=1)
+
+        # Check concatenation with both invalid units
+        d.override_units(cf.Units("foo"), inplace=1)
+        e.override_units(cf.Units("foo"), inplace=1)
+        f = cf.Data.concatenate([d, e], axis=1, relaxed_units=True)
+        with self.assertRaises(ValueError):
+            f = cf.Data.concatenate([d, e], axis=1)
+
+        e.override_units(cf.Units("foobar"), inplace=1)
+        with self.assertRaises(ValueError):
+            f = cf.Data.concatenate([d, e], axis=1, relaxed_units=True)
+
+        with self.assertRaises(ValueError):
+            f = cf.Data.concatenate([d, e], axis=1)
+
+        e.override_units(cf.Units("metre"), inplace=1)
+        with self.assertRaises(ValueError):
+            f = cf.Data.concatenate([d, e], axis=1, relaxed_units=True)
+
+        with self.assertRaises(ValueError):
+            f = cf.Data.concatenate([d, e], axis=1)
+
     def test_Data__contains__(self):
         """Test containment checking against Data."""
         d = cf.Data([[0, 1, 2], [3, 4, 5]], units="m", chunks=2)
diff --git a/cf/test/test_aggregate.py b/cf/test/test_aggregate.py
index 51db26a6b..e80a98878 100644
--- a/cf/test/test_aggregate.py
+++ b/cf/test/test_aggregate.py
@@ -322,6 +322,18 @@ class aggregateTest(unittest.TestCase):
         g = cf.aggregate([a, b, c, d], exist_all=True)
         self.assertEqual(len(g), 2)
 
+    def test_aggregate_relaxed_units(self):
+        f = cf.example_field(0)
+        bad_units = cf.Units("bad-units")
+        f.override_units(bad_units, inplace=True)
+        g = f[:2]
+        h = f[2:]
+        i = cf.aggregate([g, h], relaxed_units=True)
+        self.assertEqual(len(i), 1)
+        i = i[0]
+        self.assertEqual(i.Units.__dict__, bad_units.__dict__)
+        self.assertTrue((i.array == f.array).all())
+
 
 if __name__ == "__main__":
     print("Run date:", datetime.datetime.now())