Skip to content

Commit

Permalink
Merge pull request #785 from davidhassell/field-combination
Browse files Browse the repository at this point in the history 
Field combination
  • Loading branch information
@davidhassell
davidhassell authored Jun 18, 2024
2 parents efe80fc + b7c1723 commit 7c2acd9
Show file tree
Hide file tree
Showing 7 changed files with 297 additions and 101 deletions.
9 changes: 7 additions & 2 deletions Changelog.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,11 +3,16 @@ 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`
* 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)
Expand Down
3 changes: 1 addition & 2 deletions cf/cellmethod.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -486,8 +486,7 @@ def expand_intervals(self, inplace=False, i=False):
@_deprecated_kwarg_check("i", version="3.0.0", removed_at="4.0.0")
@_inplace_enabled(default=False)
def change_axes(self, axis_map, inplace=False, i=False):
"""Change the axes of the cell method according to a given
mapping.
"""Replace the axes of the cell method.
:Parameters:
Expand Down
265 changes: 168 additions & 97 deletions cf/field.py
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
import logging
from collections import namedtuple
from dataclasses import dataclass
from functools import reduce
from operator import mul as operator_mul
from os import sep
Expand Down Expand Up @@ -190,11 +190,29 @@
"__ge__",
)

_xxx = namedtuple(
"data_dimension", ["size", "axis", "key", "coord", "coord_type", "scalar"]
)

# _empty_set = set()
@dataclass()
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, an integer.
size: int = -1
# The domain axis identifier. E.g. 'domainaxis0'
axis: str = ""
# The coordinate constructs that characterize the axis
coords: 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


class Field(mixin.FieldDomain, mixin.PropertiesData, cfdm.Field):
Expand Down Expand Up @@ -985,80 +1003,127 @@ 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, coord = f.dimension_coordinate(
item=True, default=(None, None), filter_by_axis=(axis,)
)
if coord is not None:
# This axis of the domain has a dimension
# coordinate
identity = 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(coord, ctype, False):
identity = ctype
break

if identity is None and relaxed_identities:
identity = coord.identity(relaxed=True, default=None)
else:
key, coord = f.auxiliary_coordinate(
item=True,
default=(None, None),
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="exact",
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:
# Get the sorted identities (sorted so that
# they're comparable between fields) and their
# corresponding keys.
#
# 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:
key, dim_coord = f.dimension_coordinate(
item=True, default=(None, None), filter_by_axis=(axis,)
)
if coord is not None:
# This axis of the domain does not have a
# dimension coordinate but it does have
# exactly one 1-d auxiliary coordinate, so
# that will do.
identity = coord.identity(strict=True, default=None)
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 = coord.identity(
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),
filter_by_axis=(axis,),
axis_mode="exact",
)
if aux_coord is not None:
# 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
)
if identity is None and relaxed_identities:
identity = aux_coord.identity(
relaxed=True, default=None
)

if identity is None:
identity = i
else:
coord_type = coord.construct_type

out[identity] = _xxx(
out[identity] = _Axis_characterisation(
size=f.domain_axis(axis).get_size(),
axis=axis,
key=key,
coord=coord,
coord_type=coord_type,
scalar=(axis not in data_axes),
keys=keys,
coords=coords,
axis_in_data_axes=axis in data_axes,
)

for identity, y in tuple(out1.items()):
asdas = True
if y.scalar and identity in out0 and isinstance(identity, str):
missing_axis_inserted = False
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
field1.insert_dimension(y.axis, position=0, inplace=True)
asdas = False
missing_axis_inserted = True

if y.scalar and asdas:
if not missing_axis_inserted and not y.axis_in_data_axes:
del out1[identity]

for identity, a in tuple(out0.items()):
asdas = True
if a.scalar and identity in out1 and isinstance(identity, str):
missing_axis_inserted = False
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)
asdas = False
missing_axis_inserted = True

if a.scalar and asdas:
if not missing_axis_inserted and not a.axis_in_data_axes:
del out0[identity]

squeeze1 = []
Expand All @@ -1069,15 +1134,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=<CF DimensionCoordinate: longitude(8) degrees_east>,
# 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 = {}

Expand Down Expand Up @@ -1178,48 +1242,55 @@ def _binary_operation(self, other, method):
f"{a.size} {identity!r} axis"
)

# Ensure matching axis directions
if y.coord.direction() != a.coord.direction():
other.flip(y.axis, inplace=True)
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 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)
# Check coord refs
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)
n_refs = len(refs1)
n0_refs = len(refs0)

if n_refs != n0_refs:
raise ValueError(
f"Can't combine {self.__class__.__name__!r} with "
f"{other.__class__.__name__!r} because the coordinate "
f"references have different lengths: {n_refs} and "
f"{n0_refs}."
)
if n_refs != n0_refs:
raise ValueError(
f"Can't combine {self.__class__.__name__!r} with "
f"{other.__class__.__name__!r} because the coordinate "
f"references have different lengths: {n_refs} and "
f"{n0_refs}."
)

n_equivalent_refs = 0
for ref1 in refs1:
for ref0 in refs0[:]:
if field1._equivalent_coordinate_references(
field0, key0=ref1, key1=ref0, axis_map=axis_map
):
n_equivalent_refs += 1
refs0.remove(ref0)
break
n_equivalent_refs = 0
for ref1 in refs1:
for ref0 in refs0[:]:
if field1._equivalent_coordinate_references(
field0, key0=ref1, key1=ref0, axis_map=axis_map
):
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."
)
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."
)

# Change the domain axis sizes in field0 so that they match
# the broadcasted result data
Expand Down
Loading

0 comments on commit 7c2acd9

Please sign in to comment.