Skip to content

Commit 7c2acd9

Browse files
davidhasselldavidhassell
authored
Merge pull request #785 from davidhassell/field-combination
Field combination
2 parents efe80fc + b7c1723 commit 7c2acd9

File tree

7 files changed

+297
-101
lines changed

7 files changed

+297
-101
lines changed

‎Changelog.rst‎

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -3,11 +3,16 @@ version NEXTVERSION
33

44
**2024-??-??**
55

6+
* New method: `cf.Field.is_discrete_axis`
7+
(https://github.com/NCAS-CMS/cf-python/issues/784)
68
* Include the UM version as a field property when reading UM files
79
(https://github.com/NCAS-CMS/cf-python/issues/777)
8-
* Fix bug where `cf.example_fields` returned a `list`
9-
of Fields rather than a `Fieldlist`
10+
* Fix bug where `cf.example_fields` returned a `list` of Fields rather
11+
than a `Fieldlist`
1012
(https://github.com/NCAS-CMS/cf-python/issues/725)
13+
* Fix bug where combining UGRID fields erroneously creates an extra
14+
axis and broadcasts over it
15+
(https://github.com/NCAS-CMS/cf-python/issues/784)
1116
* Fix bug where `cf.normalize_slice` doesn't correctly
1217
handle certain cyclic slices
1318
(https://github.com/NCAS-CMS/cf-python/issues/774)

‎cf/cellmethod.py‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -486,8 +486,7 @@ def expand_intervals(self, inplace=False, i=False):
486486
@_deprecated_kwarg_check("i", version="3.0.0", removed_at="4.0.0")
487487
@_inplace_enabled(default=False)
488488
def change_axes(self, axis_map, inplace=False, i=False):
489-
"""Change the axes of the cell method according to a given
490-
mapping.
489+
"""Replace the axes of the cell method.
491490
492491
:Parameters:
493492

‎cf/field.py‎

Lines changed: 168 additions & 97 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
import logging
2-
from collections import namedtuple
2+
from dataclasses import dataclass
33
from functools import reduce
44
from operator import mul as operator_mul
55
from os import sep
@@ -190,11 +190,29 @@
190190
"__ge__",
191191
)
192192

193-
_xxx = namedtuple(
194-
"data_dimension", ["size", "axis", "key", "coord", "coord_type", "scalar"]
195-
)
196193

197-
# _empty_set = set()
194+
@dataclass()
195+
class _Axis_characterisation:
196+
"""Characterise a domain axis.
197+
198+
Used by `_binary_operation` to help with ascertaining if there is
199+
a common axis in two fields.
200+
201+
.. versionaddedd:: NEXTVERSION
202+
203+
"""
204+
205+
# The size of the axis, an integer.
206+
size: int = -1
207+
# The domain axis identifier. E.g. 'domainaxis0'
208+
axis: str = ""
209+
# The coordinate constructs that characterize the axis
210+
coords: tuple = ()
211+
# The identifiers of the coordinate
212+
# constructs. E.g. ('dimensioncoordinate1',)
213+
keys: tuple = ()
214+
# Whether or not the axis is spanned by the field's data array
215+
axis_in_data_axes: bool = True
198216

199217

200218
class Field(mixin.FieldDomain, mixin.PropertiesData, cfdm.Field):
@@ -985,80 +1003,127 @@ def _binary_operation(self, other, method):
9851003
data_axes = f.get_data_axes()
9861004
for axis in f.domain_axes(todict=True):
9871005
identity = None
988-
key = None
989-
coord = None
990-
coord_type = None
9911006

992-
key, coord = f.dimension_coordinate(
993-
item=True, default=(None, None), filter_by_axis=(axis,)
994-
)
995-
if coord is not None:
996-
# This axis of the domain has a dimension
997-
# coordinate
998-
identity = coord.identity(strict=True, default=None)
999-
if identity is None:
1000-
# Dimension coordinate has no identity, but it
1001-
# may have a recognised axis.
1002-
for ctype in ("T", "X", "Y", "Z"):
1003-
if getattr(coord, ctype, False):
1004-
identity = ctype
1005-
break
1006-
1007-
if identity is None and relaxed_identities:
1008-
identity = coord.identity(relaxed=True, default=None)
1009-
else:
1010-
key, coord = f.auxiliary_coordinate(
1011-
item=True,
1012-
default=(None, None),
1007+
if self.is_discrete_axis(axis):
1008+
# This is a discrete axis whose identity is
1009+
# inferred from all of its auxiliary coordinates
1010+
x = {}
1011+
for key, aux_coord in f.auxiliary_coordinates(
10131012
filter_by_axis=(axis,),
1014-
axis_mode="exact",
1013+
axis_mode="and",
1014+
todict=True,
1015+
).items():
1016+
identity = aux_coord.identity(
1017+
strict=True, default=None
1018+
)
1019+
if identity is None and relaxed_identities:
1020+
identity = aux_coord.identity(
1021+
relaxed=True, default=None
1022+
)
1023+
if identity is not None:
1024+
x[identity] = key
1025+
1026+
if x:
1027+
# Get the sorted identities (sorted so that
1028+
# they're comparable between fields) and their
1029+
# corresponding keys.
1030+
#
1031+
# E.g. {2:3, 4:6, 1:7} -> (1, 2, 4), (7, 3, 6)
1032+
identity, keys = tuple(zip(*sorted(x.items())))
1033+
coords = tuple(
1034+
f.auxiliary_coordinate(key) for key in keys
1035+
)
1036+
else:
1037+
identity = None
1038+
keys = ()
1039+
coords = ()
1040+
else:
1041+
key, dim_coord = f.dimension_coordinate(
1042+
item=True, default=(None, None), filter_by_axis=(axis,)
10151043
)
1016-
if coord is not None:
1017-
# This axis of the domain does not have a
1018-
# dimension coordinate but it does have
1019-
# exactly one 1-d auxiliary coordinate, so
1020-
# that will do.
1021-
identity = coord.identity(strict=True, default=None)
1044+
if dim_coord is not None:
1045+
# This non-discrete axis has a dimension
1046+
# coordinate
1047+
identity = dim_coord.identity(
1048+
strict=True, default=None
1049+
)
1050+
if identity is None:
1051+
# Dimension coordinate has no identity,
1052+
# but it may have a recognised axis.
1053+
for ctype in ("T", "X", "Y", "Z"):
1054+
if getattr(dim_coord, ctype, False):
1055+
identity = ctype
1056+
break
10221057

10231058
if identity is None and relaxed_identities:
1024-
identity = coord.identity(
1059+
identity = dim_coord.identity(
10251060
relaxed=True, default=None
10261061
)
10271062

1063+
keys = (key,)
1064+
coords = (dim_coord,)
1065+
else:
1066+
key, aux_coord = f.auxiliary_coordinate(
1067+
item=True,
1068+
default=(None, None),
1069+
filter_by_axis=(axis,),
1070+
axis_mode="exact",
1071+
)
1072+
if aux_coord is not None:
1073+
# This non-discrete axis does not have a
1074+
# dimension coordinate but it does have
1075+
# exactly one 1-d auxiliary coordinate, so
1076+
# that will do.
1077+
coords = (aux_coord,)
1078+
identity = aux_coord.identity(
1079+
strict=True, default=None
1080+
)
1081+
if identity is None and relaxed_identities:
1082+
identity = aux_coord.identity(
1083+
relaxed=True, default=None
1084+
)
1085+
10281086
if identity is None:
10291087
identity = i
1030-
else:
1031-
coord_type = coord.construct_type
10321088

1033-
out[identity] = _xxx(
1089+
out[identity] = _Axis_characterisation(
10341090
size=f.domain_axis(axis).get_size(),
10351091
axis=axis,
1036-
key=key,
1037-
coord=coord,
1038-
coord_type=coord_type,
1039-
scalar=(axis not in data_axes),
1092+
keys=keys,
1093+
coords=coords,
1094+
axis_in_data_axes=axis in data_axes,
10401095
)
10411096

10421097
for identity, y in tuple(out1.items()):
1043-
asdas = True
1044-
if y.scalar and identity in out0 and isinstance(identity, str):
1098+
missing_axis_inserted = False
1099+
if (
1100+
not y.axis_in_data_axes
1101+
and identity in out0
1102+
and isinstance(identity, str)
1103+
):
10451104
a = out0[identity]
10461105
if a.size > 1:
1106+
# Put missing axis into data axes
10471107
field1.insert_dimension(y.axis, position=0, inplace=True)
1048-
asdas = False
1108+
missing_axis_inserted = True
10491109

1050-
if y.scalar and asdas:
1110+
if not missing_axis_inserted and not y.axis_in_data_axes:
10511111
del out1[identity]
10521112

10531113
for identity, a in tuple(out0.items()):
1054-
asdas = True
1055-
if a.scalar and identity in out1 and isinstance(identity, str):
1114+
missing_axis_inserted = False
1115+
if (
1116+
not a.axis_in_data_axes
1117+
and identity in out1
1118+
and isinstance(identity, str)
1119+
):
10561120
y = out1[identity]
10571121
if y.size > 1:
1122+
# Put missing axis into data axes
10581123
field0.insert_dimension(a.axis, position=0, inplace=True)
1059-
asdas = False
1124+
missing_axis_inserted = True
10601125

1061-
if a.scalar and asdas:
1126+
if not missing_axis_inserted and not a.axis_in_data_axes:
10621127
del out0[identity]
10631128

10641129
squeeze1 = []
@@ -1069,15 +1134,14 @@ def _binary_operation(self, other, method):
10691134
axes_added_from_field1 = []
10701135

10711136
# Dictionary of size > 1 axes from field1 which will replace
1072-
# matching size 1 axes in field0. E.g. {'domainaxis1':
1073-
# data_dimension(
1074-
# size=8,
1075-
# axis='domainaxis1',
1076-
# key='dimensioncoordinate1',
1077-
# coord=<CF DimensionCoordinate: longitude(8) degrees_east>,
1078-
# coord_type='dimension_coordinate',
1079-
# scalar=False
1080-
# )
1137+
# matching size 1 axes in field0.
1138+
#
1139+
# E.g. {'domainaxis1': _Axis_characterisation(
1140+
# size=8,
1141+
# axis='domainaxis1',
1142+
# keys=('dimensioncoordinate1',),
1143+
# coords=(CF DimensionCoordinate: longitude(8) degrees_east>,),
1144+
# axis_in_data_axes=True)
10811145
# }
10821146
axes_to_replace_from_field1 = {}
10831147

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

1181-
# Ensure matching axis directions
1182-
if y.coord.direction() != a.coord.direction():
1183-
other.flip(y.axis, inplace=True)
1245+
for a_key, a_coord, y_key, y_coord in zip(
1246+
a.keys, a.coords, y.keys, y.coords
1247+
):
1248+
# Ensure matching axis directions
1249+
if y_coord.direction() != a_coord.direction():
1250+
other.flip(y.axis, inplace=True)
11841251

1185-
# Check for matching coordinate values
1186-
if not y.coord._equivalent_data(a.coord):
1187-
raise ValueError(
1188-
f"Can't combine size {y.size} {identity!r} axes with "
1189-
f"non-matching coordinate values"
1190-
)
1252+
# Check for matching coordinate values
1253+
if not y_coord._equivalent_data(a_coord):
1254+
raise ValueError(
1255+
f"Can't combine size {y.size} {identity!r} axes with "
1256+
f"non-matching coordinate values"
1257+
)
11911258

1192-
# Check coord refs
1193-
refs1 = field1.get_coordinate_reference(construct=y.key, key=True)
1194-
refs0 = field0.get_coordinate_reference(construct=a.key, key=True)
1259+
# Check coord refs
1260+
refs1 = field1.get_coordinate_reference(
1261+
construct=y_key, key=True
1262+
)
1263+
refs0 = field0.get_coordinate_reference(
1264+
construct=a_key, key=True
1265+
)
11951266

1196-
n_refs = len(refs1)
1197-
n0_refs = len(refs0)
1267+
n_refs = len(refs1)
1268+
n0_refs = len(refs0)
11981269

1199-
if n_refs != n0_refs:
1200-
raise ValueError(
1201-
f"Can't combine {self.__class__.__name__!r} with "
1202-
f"{other.__class__.__name__!r} because the coordinate "
1203-
f"references have different lengths: {n_refs} and "
1204-
f"{n0_refs}."
1205-
)
1270+
if n_refs != n0_refs:
1271+
raise ValueError(
1272+
f"Can't combine {self.__class__.__name__!r} with "
1273+
f"{other.__class__.__name__!r} because the coordinate "
1274+
f"references have different lengths: {n_refs} and "
1275+
f"{n0_refs}."
1276+
)
12061277

1207-
n_equivalent_refs = 0
1208-
for ref1 in refs1:
1209-
for ref0 in refs0[:]:
1210-
if field1._equivalent_coordinate_references(
1211-
field0, key0=ref1, key1=ref0, axis_map=axis_map
1212-
):
1213-
n_equivalent_refs += 1
1214-
refs0.remove(ref0)
1215-
break
1278+
n_equivalent_refs = 0
1279+
for ref1 in refs1:
1280+
for ref0 in refs0[:]:
1281+
if field1._equivalent_coordinate_references(
1282+
field0, key0=ref1, key1=ref0, axis_map=axis_map
1283+
):
1284+
n_equivalent_refs += 1
1285+
refs0.remove(ref0)
1286+
break
12161287

1217-
if n_equivalent_refs != n_refs:
1218-
raise ValueError(
1219-
f"Can't combine {self.__class__.__name__!r} with "
1220-
f"{other.__class__.__name__!r} because the fields have "
1221-
"incompatible coordinate references."
1222-
)
1288+
if n_equivalent_refs != n_refs:
1289+
raise ValueError(
1290+
f"Can't combine {self.__class__.__name__!r} with "
1291+
f"{other.__class__.__name__!r} because the fields "
1292+
"have incompatible coordinate references."
1293+
)
12231294

12241295
# Change the domain axis sizes in field0 so that they match
12251296
# the broadcasted result data

0 commit comments

Comments
 (0)