Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix bug where cf.Field.subspace doesn't always correctly handle global cyclic subspaces #830

Merged
merged 15 commits into from
Nov 14, 2024
4 changes: 3 additions & 1 deletion Changelog.rst
Original file line number Diff line number Diff line change
Expand Up @@ -33,13 +33,15 @@ version NEXTVERSION
* Fix bug where `cf.normalize_slice` doesn't correctly
handle certain cyclic slices
(https://github.com/NCAS-CMS/cf-python/issues/774)
* Fix bug where `cf.Field.subspace` doesn't always correctly handle
global or near-global cyclic subspaces
(https://github.com/NCAS-CMS/cf-python/issues/828)
* New dependency: ``h5netcdf>=1.3.0``
* New dependency: ``h5py>=3.10.0``
* New dependency: ``s3fs>=2024.2.0``
* Changed dependency: ``1.11.2.0<=cfdm<1.11.3.0``
* Changed dependency: ``cfunits>=3.3.7``


----

version 3.16.2
Expand Down
1 change: 0 additions & 1 deletion cf/cfimplementation.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,6 @@
TiePointIndex,
)
from .data import Data

from .data.array import (
BoundsFromNodesArray,
CellConnectivityArray,
Expand Down
1 change: 0 additions & 1 deletion cf/data/data.py
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,6 @@
from ..units import Units
from .collapse import Collapse
from .creation import generate_axis_identifiers, to_dask

from .dask_utils import (
_da_ma_allclose,
cf_asanyarray,
Expand Down
154 changes: 153 additions & 1 deletion cf/dimensioncoordinate.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,11 @@
_inplace_enabled,
_inplace_enabled_define_and_cleanup,
)
from .functions import _DEPRECATION_ERROR_ATTRIBUTE, _DEPRECATION_ERROR_KWARGS
from .functions import (
_DEPRECATION_ERROR_ATTRIBUTE,
_DEPRECATION_ERROR_KWARGS,
bounds_combination_mode,
)
from .timeduration import TimeDuration
from .units import Units

Expand Down Expand Up @@ -246,6 +250,154 @@ def increasing(self):
"""
return self.direction()

@_inplace_enabled(default=False)
def anchor(self, value, cell=False, parameters=None, inplace=False):
"""Anchor the coordinate values.

By default, the coordinate values are transformed so that the
first coordinate is the closest to *value* from above (below)
for increasing (decreasing) coordinates.

If the *cell* parameter is True, then the coordinate values
are transformed so that the first cell either contains
*value*; or is the closest to cell to *value* from above
(below) for increasing (decreasing) coordinates.

.. versionadded:: NEXTVERSION

.. seealso:: `period`, `roll`

:Parameters:

value: scalar array_like
Anchor the coordinate values for the selected cyclic
axis to the *value*. May be any numeric scalar object
that can be converted to a `Data` object (which
includes `numpy` and `Data` objects). If *value* has
units then they must be compatible with those of the
coordinates, otherwise it is assumed to have the same
units as the coordinates.

The coordinate values are transformed so the first
corodinate is the closest to *value* from above (for
increasing coordinates), or the closest to *value* from
above (for decreasing coordinates)

* Increasing coordinates with positive period, P,
are transformed so that *value* lies in the
half-open range (L-P, F], where F and L are the
transformed first and last coordinate values,
respectively.

..

* Decreasing coordinates with positive period, P,
are transformed so that *value* lies in the
half-open range (L+P, F], where F and L are the
transformed first and last coordinate values,
respectively.

*Parameter example:*
If the original coordinates are ``0, 5, ..., 355``
(evenly spaced) and the period is ``360`` then
``value=0`` implies transformed coordinates of ``0,
5, ..., 355``; ``value=-12`` implies transformed
coordinates of ``-10, -5, ..., 345``; ``value=380``
implies transformed coordinates of ``380, 385, ...,
735``.

*Parameter example:*
If the original coordinates are ``355, 350, ..., 0``
(evenly spaced) and the period is ``360`` then
``value=355`` implies transformed coordinates of
``355, 350, ..., 0``; ``value=0`` implies
transformed coordinates of ``0, -5, ..., -355``;
``value=392`` implies transformed coordinates of
``390, 385, ..., 35``.

cell: `bool`, optional
If True, then the coordinate values are transformed so
that the first cell either contains *value*, or is the
closest to cell to *value* from above (below) for
increasing (decreasing) coordinates.

If False (the default) then the coordinate values are
transformed so that the first coordinate is the closest
to *value* from above (below) for increasing
(decreasing) coordinates.

parameters: `dict`, optional
If a `dict` is provided then it will be updated
in-place with parameters which describe thethe
anchoring process.

{{inplace: `bool`, optional}}

:Returns:

`{{class}}` or `None`
The anchored dimension coordinates, or `None` if the
operation was in-place.

"""
d = _inplace_enabled_define_and_cleanup(self)

period = d.period()
if period is None:
raise ValueError(f"Cyclic {d!r} has no period")

value = d._Data.asdata(value)
if not value.Units:
value = value.override_units(d.Units)
elif not value.Units.equivalent(d.Units):
raise ValueError(
f"Anchor value has incompatible units: {value.Units!r}"
)

if cell:
c = d.upper_bounds.persist()
else:
d.persist(inplace=True)
c = d.get_data(_fill_value=False)

if d.increasing:
# Adjust value so it's in the range [c[0], c[0]+period)
n = ((c[0] - value) / period).ceil()
value1 = value + n * period
shift = c.size - np.argmax((c - value1 >= 0).array)
d.roll(0, shift, inplace=True)
if cell:
d0 = d[0].upper_bounds
else:
d0 = d.get_data(_fill_value=False)[0]

n = ((value - d0) / period).ceil()
else:
# Adjust value so it's in the range (c[0]-period, c[0]]
n = ((c[0] - value) / period).floor()
value1 = value + n * period
shift = c.size - np.argmax((value1 - c >= 0).array)
d.roll(0, shift, inplace=True)
if cell:
d0 = d[0].upper_bounds
else:
d0 = d.get_data(_fill_value=False)[0]

n = ((value - d0) / period).floor()

n.persist(inplace=True)
if n:
nperiod = n * period
with bounds_combination_mode("OR"):
d += nperiod
else:
nperiod = 0

if parameters is not None:
parameters.update({"shift": shift, "nperiod": nperiod})

return d

def direction(self):
"""Return True if the dimension coordinate values are
increasing, otherwise return False.
Expand Down
140 changes: 44 additions & 96 deletions cf/mixin/fielddomain.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@
bounds_combination_mode,
normalize_slice,
)
from ..query import Query
from ..query import Query, wi
from ..units import Units

logger = logging.getLogger(__name__)
Expand Down Expand Up @@ -440,50 +440,41 @@ def _indices(self, config, data_axes, ancillary_mask, kwargs):
if debug:
logger.debug(" 1-d CASE 2:") # pragma: no cover

size = item.size
if item.increasing:
anchor0 = value.value[0]
anchor1 = value.value[1]
anchor = value.value[0]
else:
anchor0 = value.value[1]
anchor1 = value.value[0]

a = self.anchor(axis, anchor0, dry_run=True)["roll"]
b = self.flip(axis).anchor(axis, anchor1, dry_run=True)[
"roll"
]

size = item.size
if abs(anchor1 - anchor0) >= item.period():
if value.operator == "wo":
set_start_stop = 0
else:
set_start_stop = -a

start = set_start_stop
stop = set_start_stop
elif a + b == size:
b = self.anchor(axis, anchor1, dry_run=True)["roll"]
if (b == a and value.operator == "wo") or not (
b == a or value.operator == "wo"
):
set_start_stop = -a
else:
set_start_stop = 0
anchor = value.value[1]

item = item.persist()
parameters = {}
item = item.anchor(anchor, parameters=parameters)
n = np.roll(np.arange(size), parameters["shift"], 0)
if value.operator == "wi":
sadielbartholomew marked this conversation as resolved.
Show resolved Hide resolved
n = n[item == value]
if not n.size:
raise ValueError(
f"No indices found from: {identity}={value!r}"
)

start = set_start_stop
stop = set_start_stop
start = n[0]
stop = n[-1] + 1
else:
if value.operator == "wo":
start = b - size
stop = -a + size
else:
start = -a
stop = b - size
# "wo" operator
n = n[item == wi(*value.value)]
if n.size == size:
raise ValueError(
f"No indices found from: {identity}={value!r}"
)

if start == stop == 0:
raise ValueError(
f"No indices found from: {identity}={value!r}"
)
if n.size:
start = n[-1] + 1
stop = start - n.size
else:
start = size - parameters["shift"]
stop = start + size
if stop > size:
stop -= size

index = slice(start, stop, 1)

Expand Down Expand Up @@ -1287,77 +1278,34 @@ def anchor(
self, "anchor", kwargs
) # pragma: no cover

da_key, axis = self.domain_axis(axis, item=True)
axis = self.domain_axis(axis, key=True)

if dry_run:
f = self
else:
f = _inplace_enabled_define_and_cleanup(self)

dim = f.dimension_coordinate(filter_by_axis=(da_key,), default=None)
dim = f.dimension_coordinate(filter_by_axis=(axis,), default=None)
if dim is None:
raise ValueError(
"Can't shift non-cyclic "
f"{f.constructs.domain_axis_identity(da_key)!r} axis"
f"{f.constructs.domain_axis_identity(axis)!r} axis"
)

period = dim.period()
if period is None:
raise ValueError(f"Cyclic {dim.identity()!r} axis has no period")

value = f._Data.asdata(value)
if not value.Units:
value = value.override_units(dim.Units)
elif not value.Units.equivalent(dim.Units):
raise ValueError(
f"Anchor value has incompatible units: {value.Units!r}"
)

axis_size = axis.get_size()

if axis_size <= 1:
# Don't need to roll a size one axis
if dry_run:
return {"axis": da_key, "roll": 0, "nperiod": 0}

return f

c = dim.get_data(_fill_value=False)

if dim.increasing:
# Adjust value so it's in the range [c[0], c[0]+period)
n = ((c[0] - value) / period).ceil()
value1 = value + n * period

shift = axis_size - np.argmax((c - value1 >= 0).array)
if not dry_run:
f.roll(da_key, shift, inplace=True)

# Re-get dim
dim = f.dimension_coordinate(filter_by_axis=(da_key,))
# TODO CHECK n for dry run or not
n = ((value - dim.data[0]) / period).ceil()
else:
# Adjust value so it's in the range (c[0]-period, c[0]]
n = ((c[0] - value) / period).floor()
value1 = value + n * period

shift = axis_size - np.argmax((value1 - c >= 0).array)

if not dry_run:
f.roll(da_key, shift, inplace=True)

# Re-get dim
dim = f.dimension_coordinate(filter_by_axis=(da_key,))
# TODO CHECK n for dry run or not
n = ((value - dim.data[0]) / period).floor()
parameters = {"axis": axis}
dim = dim.anchor(value, parameters=parameters)

if dry_run:
return {"axis": da_key, "roll": shift, "nperiod": n * period}
return parameters

f.roll(axis, parameters["shift"], inplace=True)

if n:
if parameters["nperiod"]:
# Get the rolled dimension coordinate and adjust the
# values by the non-zero integer multiple of 'period'
dim = f.dimension_coordinate(filter_by_axis=(axis,))
with bounds_combination_mode("OR"):
dim += n * period
dim += parameters["nperiod"]

return f

Expand Down
1 change: 0 additions & 1 deletion cf/regrid/regrid.py
Original file line number Diff line number Diff line change
Expand Up @@ -2465,7 +2465,6 @@ def create_esmpy_weights(
from netCDF4 import Dataset

from .. import __version__

from ..data.array.locks import netcdf_lock

if (
Expand Down
Loading