diff --git a/Changelog.rst b/Changelog.rst index 14723cf7a4..99d4c88c4b 100644 --- a/Changelog.rst +++ b/Changelog.rst @@ -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 diff --git a/cf/cfimplementation.py b/cf/cfimplementation.py index db8601627c..118e73ce7a 100644 --- a/cf/cfimplementation.py +++ b/cf/cfimplementation.py @@ -26,7 +26,6 @@ TiePointIndex, ) from .data import Data - from .data.array import ( BoundsFromNodesArray, CellConnectivityArray, diff --git a/cf/data/data.py b/cf/data/data.py index 61abf86bc7..d965a09473 100644 --- a/cf/data/data.py +++ b/cf/data/data.py @@ -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, diff --git a/cf/dimensioncoordinate.py b/cf/dimensioncoordinate.py index 347a0d7dd5..3124035e66 100644 --- a/cf/dimensioncoordinate.py +++ b/cf/dimensioncoordinate.py @@ -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 @@ -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. diff --git a/cf/mixin/fielddomain.py b/cf/mixin/fielddomain.py index 1e9e1f67e1..833f86ad3a 100644 --- a/cf/mixin/fielddomain.py +++ b/cf/mixin/fielddomain.py @@ -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__) @@ -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": + 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) @@ -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 diff --git a/cf/regrid/regrid.py b/cf/regrid/regrid.py index a06d1cc960..eb7ee6656a 100644 --- a/cf/regrid/regrid.py +++ b/cf/regrid/regrid.py @@ -2465,7 +2465,6 @@ def create_esmpy_weights( from netCDF4 import Dataset from .. import __version__ - from ..data.array.locks import netcdf_lock if ( diff --git a/cf/test/test_DimensionCoordinate.py b/cf/test/test_DimensionCoordinate.py index 9244a202b1..b39181b245 100644 --- a/cf/test/test_DimensionCoordinate.py +++ b/cf/test/test_DimensionCoordinate.py @@ -737,6 +737,70 @@ def test_DimensionCoordinate_cell_characteristics(self): with self.assertRaises(ValueError): d.get_cell_characteristics() + def test_DimensionCoordinate_anchor(self): + """Test the DimensionCoordinate.anchor""" + f = cf.example_field(0) + d = f.dimension_coordinate("X") + + self.assertEqual(d.period(), 360) + + e = d.anchor(-1) + self.assertIsInstance(e, cf.DimensionCoordinate) + self.assertTrue(e.equals(d)) + + # Increasing + e = d.anchor(15) + self.assertTrue(e[0].equals(d[0])) + self.assertEqual(e[0].array, d[0].array) + + e = d.anchor(30) + self.assertEqual(e[0].array, d[1].array) + + e = d.anchor(361) + self.assertEqual(e[0].array, d[0].array + 360) + + e = d.anchor(721) + self.assertEqual(e[0].array, d[0].array + 720) + + e = d.anchor(15, cell=True) + self.assertEqual(e[0].array, d[0].array) + + e = d.anchor(30, cell=True) + self.assertEqual(e[0].array, d[0].array) + + e = d.anchor(361, cell=True) + self.assertEqual(e[0].array, d[0].array + 360) + + e = d.anchor(721, cell=True) + self.assertEqual(e[0].array, d[0].array + 720) + + # Decreasing + d = d[::-1] + + e = d.anchor(721) + self.assertEqual(e[0].array, d[0].array + 360) + + e = d.anchor(361) + self.assertEqual(e[0].array, d[0].array) + + e = d.anchor(30) + self.assertEqual(e[0].array, d[-1].array) + + e = d.anchor(15) + self.assertEqual(e[0].array, d[0].array - 360) + + e = d.anchor(721, cell=True) + self.assertEqual(e[0].array, d[0].array + 360) + + e = d.anchor(361, cell=True) + self.assertEqual(e[0].array, d[0].array) + + e = d.anchor(30, cell=True) + self.assertEqual(e[0].array, d[0].array - 360) + + e = d.anchor(15, cell=True) + self.assertEqual(e[0].array, d[0].array - 360) + if __name__ == "__main__": print("Run date:", datetime.datetime.now()) diff --git a/cf/test/test_Field.py b/cf/test/test_Field.py index 92392b29f5..b1df9ef1f1 100644 --- a/cf/test/test_Field.py +++ b/cf/test/test_Field.py @@ -1214,6 +1214,30 @@ def test_Field_indices(self): (x == [-80, -40, 0, 40, 80, 120, 160, 200, 240.0]).all() ) + indices = f.indices(grid_longitude=cf.wi(-90, 270)) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 9)) + x = g.dimension_coordinate("X").array + self.assertTrue( + (x == [-80, -40, 0, 40, 80, 120, 160, 200, 240.0]).all() + ) + + indices = f.indices(grid_longitude=cf.wi(-180, 180)) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 9)) + x = g.dimension_coordinate("X").array + self.assertTrue( + (x == [-160, -120, -80, -40, 0, 40, 80, 120, 160]).all() + ) + + indices = f.indices(grid_longitude=cf.wi(-170, 170)) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 9)) + x = g.dimension_coordinate("X").array + self.assertTrue( + (x == [-160, -120, -80, -40, 0, 40, 80, 120, 160]).all() + ) + with self.assertRaises(ValueError): # No X coordinate values lie inside the range [90, 100] f.indices(grid_longitude=cf.wi(90, 100)) @@ -1337,6 +1361,34 @@ def test_Field_indices(self): x = g.dimension_coordinate("X").array self.assertTrue((x == [-200, -160, -120, -80, -40, 0, 40]).all()) + indices = f.indices(grid_longitude=cf.wo(1, 5)) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 9)) + x = g.dimension_coordinate("X").array + self.assertTrue( + (x == [-320, -280, -240, -200, -160, -120, -80, -40, 0]).all() + ) + + indices = f.indices(grid_longitude=cf.wo(41, 45)) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 9)) + x = g.dimension_coordinate("X").array + self.assertTrue( + (x == [-280, -240, -200, -160, -120, -80, -40, 0, 40]).all() + ) + + indices = f.indices(grid_longitude=cf.wo(-5, -1)) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 9)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == [0, 40, 80, 120, 160, 200, 240, 280, 320]).all()) + + indices = f.indices(grid_longitude=cf.wo(-45, -41)) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 9)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == [-40, 0, 40, 80, 120, 160, 200, 240, 280]).all()) + with self.assertRaises(ValueError): # No X coordinate values lie outside the range [-90, 370] f.indices(grid_longitude=cf.wo(-90, 370)) @@ -1372,6 +1424,33 @@ def test_Field_indices(self): x = g.dimension_coordinate("X").array self.assertTrue((x == [80, 120, 160, 200, 240, 280, 320]).all()) + indices = f.indices(grid_longitude=cf.wo(-45, 45)) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 6)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == [80, 120, 160, 200, 240, 280]).all()) + + indices = f.indices(grid_longitude=cf.wo(35, 85)) + g = f[indices] + x = g.dimension_coordinate("X").array + self.assertEqual(g.shape, (1, 10, 7)) + self.assertTrue((x == [-240, -200, -160, -120, -80, -40, 0]).all()) + + indices = f.indices(grid_longitude=cf.wo(35, 85)) + g = f[indices] + x = g.dimension_coordinate("X").array + self.assertEqual(g.shape, (1, 10, 7)) + self.assertTrue((x == [-240, -200, -160, -120, -80, -40, 0]).all()) + + with self.assertRaises(ValueError): + f.indices(grid_longitude=cf.wo(0, 360)) + + with self.assertRaises(ValueError): + f.indices(grid_longitude=cf.wo(-90, 270)) + + with self.assertRaises(ValueError): + f.indices(grid_longitude=cf.wo(-180, 180)) + # 2-d lon = f.construct("longitude").array lon = np.transpose(lon)