diff --git a/Changelog.rst b/Changelog.rst index 506baebb63..3a4b567d30 100644 --- a/Changelog.rst +++ b/Changelog.rst @@ -3,11 +3,18 @@ version NEXTVERSION **2024-??-??** +* New method: `cf.Field.is_discrete_axis` + (https://github.com/NCAS-CMS/cf-python/issues/784) * Include the UM version as a field property when reading UM files (https://github.com/NCAS-CMS/cf-python/issues/777) -* Fix bug where `cf.example_fields` returned a `list` - of Fields rather than a `Fieldlist` +* Include the UM version as a field property when reading UM files + (https://github.com/NCAS-CMS/cf-python/issues/777) +* Fix bug where `cf.example_fields` returned a `list` of Fields rather + than a `Fieldlist` (https://github.com/NCAS-CMS/cf-python/issues/725) +* Fix bug where combining UGRID fields erroneously creates an extra + axis and broadcasts over it + (https://github.com/NCAS-CMS/cf-python/issues/784) * Fix bug where `cf.normalize_slice` doesn't correctly handle certain cyclic slices (https://github.com/NCAS-CMS/cf-python/issues/774) diff --git a/cf/field.py b/cf/field.py index 0c38be525f..84b70f8308 100644 --- a/cf/field.py +++ b/cf/field.py @@ -191,26 +191,29 @@ "__ge__", ) -#_xxx = namedtuple( -# "data_dimension", ["size", "axis", "keys", "coords", "coord_type", "scalar#"] -#) @dataclass() -class _data_dimension: - """TODO +class _Axis_characterisation: + """Characterise a domain axis. + + Used by `_binary_operation` to help with ascertaining if there is + a common axis in two fields. .. versionaddedd:: NEXTVERSION """ + + # The size of the axis, a positive integer. size: int = -1 + # The domain axis identifier. E.g. 'domainaxis0' axis: str = "" - keys: tuple = () + # The coordinate constructs that characterize the axis coords: tuple = () - coord_type: tuple = () + # The identifiers of the coordinate + # constructs. E.g. ('dimensioncoordinate1',) + keys: tuple = () + # Whether or not the axis is spanned by the field's data array axis_in_data_axes: bool = True - - -# _empty_set = set() class Field(mixin.FieldDomain, mixin.PropertiesData, cfdm.Field): @@ -1001,56 +1004,62 @@ def _binary_operation(self, other, method): data_axes = f.get_data_axes() for axis in f.domain_axes(todict=True): identity = None - key = None - coord = None - coord_type = None - key, dim_coord = f.dimension_coordinate( - item=True, default=(None, None), filter_by_axis=(axis,) - ) - if dim_coord is not None: - # This axis of the domain has a dimension - # coordinate - identity = dim_coord.identity(strict=True, default=None) - if identity is None: - # Dimension coordinate has no identity, but it - # may have a recognised axis. - for ctype in ("T", "X", "Y", "Z"): - if getattr(dim_coord, ctype, False): - identity = ctype - break - - if identity is None and relaxed_identities: - identity = dim_coord.identity(relaxed=True, default=None) - - keys = (key,) - coords = (dim_coord,) + if self.is_discrete_axis(axis): + # This is a discrete axis whose identity is + # inferred from all of its auxiliary coordinates + x = {} + for key, aux_coord in f.auxiliary_coordinates( + filter_by_axis=(axis,), + axis_mode="and", + todict=True, + ).items(): + identity = aux_coord.identity( + strict=True, default=None + ) + if identity is None and relaxed_identities: + identity = aux_coord.identity( + relaxed=True, default=None + ) + if identity is not None: + x[identity] = key + + if x: + # E.g. {2:3, 4:6, 1:7} -> (1, 2, 4), (7, 3, 6) + identity, keys = tuple(zip(*sorted(x.items()))) + coords = tuple( + f.auxiliary_coordinate(key) for key in keys + ) + else: + identity = None + keys = () + coords = () else: - discrete_axis = f.has_property('featureType') or f.domain_topologies(todict=True) - if discrete_axis: - x = {} - for key, aux_coord in f.auxiliary_coordinates( - filter_by_axis=(axis,), - axis_mode="exact", todict=True - ).items(): - identity = aux_coord.identity(strict=True, default=None) - if identity is None and relaxed_identities: - identity = aux_coord.identity( - relaxed=True, default=None - ) - if identity is not None: - x[identity] = key - - if x: - # E.g. {2:3, 4:6, 1:7} -> (1, 2, 4), (7, 3, 6) - identity, keys = tuple(zip(*sorted(x.items()))) - coords = tuple(f.auxiliary_coordinate(key) - for key in keys) - else: - identity = None - keys = None - coords = None - else: + key, dim_coord = f.dimension_coordinate( + item=True, default=(None, None), filter_by_axis=(axis,) + ) + if dim_coord is not None: + # This non-discrete axis has a dimension + # coordinate + identity = dim_coord.identity( + strict=True, default=None + ) + if identity is None: + # Dimension coordinate has no identity, + # but it may have a recognised axis. + for ctype in ("T", "X", "Y", "Z"): + if getattr(dim_coord, ctype, False): + identity = ctype + break + + if identity is None and relaxed_identities: + identity = dim_coord.identity( + relaxed=True, default=None + ) + + keys = (key,) + coords = (dim_coord,) + else: key, aux_coord = f.auxiliary_coordinate( item=True, default=(None, None), @@ -1058,35 +1067,37 @@ def _binary_operation(self, other, method): axis_mode="exact", ) if aux_coord is not None: - # This non-discrete axis of the domain - # does not have a dimension coordinate but - # it does have exactly one 1-d auxiliary - # coordinate, so that will do. + # This non-discrete axis does not have a + # dimension coordinate but it does have + # exactly one 1-d auxiliary coordinate, so + # that will do. coords = (aux_coord,) - identity = aux_coord.identity(strict=True, default=None) + identity = aux_coord.identity( + strict=True, default=None + ) if identity is None and relaxed_identities: identity = aux_coord.identity( relaxed=True, default=None ) - if identity is None: identity = i - else: - coord_type = coords[0].construct_type - out[identity] = _data_dimension( + out[identity] = _Axis_characterisation( size=f.domain_axis(axis).get_size(), axis=axis, keys=keys, coords=coords, - coord_type=coord_type, axis_in_data_axes=axis in data_axes, ) for identity, y in tuple(out1.items()): missing_axis_inserted = False - if not y.axis_in_data_axes and identity in out0 and isinstance(identity, str): + if ( + not y.axis_in_data_axes + and identity in out0 + and isinstance(identity, str) + ): a = out0[identity] if a.size > 1: # Put missing axis into data axes @@ -1098,14 +1109,18 @@ def _binary_operation(self, other, method): for identity, a in tuple(out0.items()): missing_axis_inserted = False - if not a.axis_in_data_axes and identity in out1 and isinstance(identity, str): + if ( + not a.axis_in_data_axes + and identity in out1 + and isinstance(identity, str) + ): y = out1[identity] if y.size > 1: # Put missing axis into data axes field0.insert_dimension(a.axis, position=0, inplace=True) missing_axis_inserted = True - if not missing_axis_inserted and not a.axis_in_data_axes: + if not missing_axis_inserted and not a.axis_in_data_axes: del out0[identity] squeeze1 = [] @@ -1116,15 +1131,14 @@ def _binary_operation(self, other, method): axes_added_from_field1 = [] # Dictionary of size > 1 axes from field1 which will replace - # matching size 1 axes in field0. E.g. {'domainaxis1': - # data_dimension( - # size=8, - # axis='domainaxis1', - # key='dimensioncoordinate1', - # coord=, - # coord_type='dimension_coordinate', - # scalar=False - # ) + # matching size 1 axes in field0. + # + # E.g. {'domainaxis1': _Axis_characterisation( + # size=8, + # axis='domainaxis1', + # keys=('dimensioncoordinate1',), + # coords=(CF DimensionCoordinate: longitude(8) degrees_east>,), + # axis_in_data_axes=True) # } axes_to_replace_from_field1 = {} @@ -1225,25 +1239,31 @@ def _binary_operation(self, other, method): f"{a.size} {identity!r} axis" ) - for a_key, a_coord, y_key, y_coord in zip(a.keys, a.coords, y.keys, y.coords): + for a_key, a_coord, y_key, y_coord in zip( + a.keys, a.coords, y.keys, y.coords + ): # Ensure matching axis directions if y_coord.direction() != a_coord.direction(): other.flip(y.axis, inplace=True) - + # Check for matching coordinate values if not y_coord._equivalent_data(a_coord): raise ValueError( f"Can't combine size {y.size} {identity!r} axes with " f"non-matching coordinate values" ) - + # Check coord refs - refs1 = field1.get_coordinate_reference(construct=y_key, key=True) - refs0 = field0.get_coordinate_reference(construct=a_key, key=True) - + refs1 = field1.get_coordinate_reference( + construct=y_key, key=True + ) + refs0 = field0.get_coordinate_reference( + construct=a_key, key=True + ) + n_refs = len(refs1) n0_refs = len(refs0) - + if n_refs != n0_refs: raise ValueError( f"Can't combine {self.__class__.__name__!r} with " @@ -1251,7 +1271,7 @@ def _binary_operation(self, other, method): f"references have different lengths: {n_refs} and " f"{n0_refs}." ) - + n_equivalent_refs = 0 for ref1 in refs1: for ref0 in refs0[:]: @@ -1261,12 +1281,12 @@ def _binary_operation(self, other, method): n_equivalent_refs += 1 refs0.remove(ref0) break - + if n_equivalent_refs != n_refs: raise ValueError( f"Can't combine {self.__class__.__name__!r} with " - f"{other.__class__.__name__!r} because the fields have " - "incompatible coordinate references." + f"{other.__class__.__name__!r} because the fields " + "have incompatible coordinate references." ) # Change the domain axis sizes in field0 so that they match @@ -1473,12 +1493,7 @@ def _conform_cell_methods(self): axis_map = {} for cm in self.cell_methods(todict=True).values(): - - axes = cm.get_axis_identities(None) - if axes is None: - axes = cm.get_axes(None) - - for axis in axes: + for axis in cm.get_axes(()): if axis in axis_map: continue diff --git a/cf/mixin/fielddomain.py b/cf/mixin/fielddomain.py index 665658bb66..45208eea20 100644 --- a/cf/mixin/fielddomain.py +++ b/cf/mixin/fielddomain.py @@ -2347,6 +2347,98 @@ def iscyclic(self, *identity, **filter_kwargs): return axis in self.cyclic() + def is_discrete_axis(self, *identity, **filter_kwargs): + """Return True if the given axis is discrete. + + In general, a discrete axis is any axis that does not + correspond to a continuous physical quantity, but only the + following types of discrete axis are identified here: + + * The feature instance axis of a discreate sampling geometry + (DSG) domain. + + * An axis spanned by the domain topology construct of an + unstructured grid. + + * The axis with geometry cells. + + .. versionaddedd:: NEXTVERSION + + .. seealso:: `domain_axis` + + :Parameters: + + identity: `tuple`, optional + Select domain axis constructs that have an identity, + defined by their `!identities` methods, that matches + any of the given values. + + Additionally, the values are matched against construct + identifiers, with or without the ``'key%'`` prefix. + + Additionally, if for a given ``value``, + ``f.coordinates(value, filter_by_naxes=(1,))`` returns + 1-d coordinate constructs that all span the same + domain axis construct then that domain axis construct + is selected. See `coordinates` for details. + + Additionally, if there is a `Field` data array and a + value matches the integer position of an array + dimension, then the corresponding domain axis + construct is selected. + + If no values are provided then all domain axis + constructs are selected. + + {{value match}} + + {{displayed identity}} + + {{filter_kwargs: optional}} + + **Examples** + + >>> f = cf.example_{{class_lower}}(8) + >>> f.is_discrete_axis('X') + True + >>> f.is_discrete_axis('T') + False + + """ + # Get the axis key + axis = self.domain_axis(*identity, key=True, **filter_kwargs) + + # DSG + if self.has_property("featureType") and not self.dimension_coordinate( + filter_by_axis=(axis,), default=False + ): + ctypes = ("X", "Y", "T") + n = 0 + for aux in self.auxiliary_coordinates( + filter_by_axis=(axis,), axis_mode="and", todict=True + ).values(): + if aux.ctype in ctypes: + n += 1 + + if n == len(ctypes): + return True + + # UGRID + if self.domain_topologies( + filter_by_axis=(axis,), axis_mode="exact", todict=True + ): + return True + + # Geometries + for aux in self.auxiliary_coordinates( + filter_by_axis=(axis,), axis_mode="exact", todict=True + ).values(): + if aux.get_geometry(None): + return True + + # Still here? Then the axis is not discrete. + return False + def match_by_rank(self, *ranks): """Whether or not the number of domain axis constructs satisfies conditions. diff --git a/cf/test/test_Field.py b/cf/test/test_Field.py index 62048773e2..13366c52f2 100644 --- a/cf/test/test_Field.py +++ b/cf/test/test_Field.py @@ -749,6 +749,11 @@ def test_Field__add__(self): with self.assertRaises(TypeError): f + ("a string",) + # Addition with a UGRID discrete axis + f = cf.example_field(8) + g = f + f + self.assertEqual(g.shape, f.shape) + def test_Field__mul__(self): f = self.f.copy().squeeze() @@ -2883,6 +2888,28 @@ def test_Field_cyclic_iscyclic(self): self.assertEqual(f2.cyclic(), set(("domainaxis2",))) self.assertTrue(f2.iscyclic("X")) + def test_Field_is_discrete_axis(self): + """Test the `is_discrete_axis` Field method.""" + # No discrete axes + f = cf.example_field(1) + for axis in f.domain_axes(todict=True): + self.assertFalse(f.is_discrete_axis(axis)) + + # UGRID + f = cf.example_field(8) + self.assertTrue(f.is_discrete_axis("X")) + self.assertFalse(f.is_discrete_axis("T")) + + # DSG + f = cf.example_field(4) + self.assertTrue(f.is_discrete_axis("cf_role=timeseries_id")) + self.assertFalse(f.is_discrete_axis("ncdim%timeseries")) + + # Geometries + f = cf.example_field(6) + self.assertTrue(f.is_discrete_axis("cf_role=timeseries_id")) + self.assertFalse(f.is_discrete_axis("time")) + if __name__ == "__main__": print("Run date:", datetime.datetime.now()) diff --git a/docs/source/class/cf.Domain.rst b/docs/source/class/cf.Domain.rst index 3117abd118..e3d076733f 100644 --- a/docs/source/class/cf.Domain.rst +++ b/docs/source/class/cf.Domain.rst @@ -207,6 +207,7 @@ Domain axes ~cf.Domain.direction ~cf.Domain.directions ~cf.Domain.iscyclic + ~cf.Domain.is_discrete_axis Subspacing ---------- diff --git a/docs/source/class/cf.Field.rst b/docs/source/class/cf.Field.rst index c0a0275daa..8bbe0f5743 100644 --- a/docs/source/class/cf.Field.rst +++ b/docs/source/class/cf.Field.rst @@ -518,6 +518,7 @@ Domain axes ~cf.Field.direction ~cf.Field.directions ~cf.Field.iscyclic + ~cf.Field.is_discrete_axis ~cf.Field.isperiodic ~cf.Field.item_axes ~cf.Field.items_axes