From dae03c9605a87ae6f201950581b625322b749664 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Tue, 9 Apr 2024 14:58:26 +0100 Subject: [PATCH 01/30] dev --- cf/field.py | 5 +-- cf/mixin/fielddomain.py | 77 +++++++++++++++++++++++++++++++++++++---- 2 files changed, 73 insertions(+), 9 deletions(-) diff --git a/cf/field.py b/cf/field.py index 165d85d427..56bdc3ec31 100644 --- a/cf/field.py +++ b/cf/field.py @@ -8877,8 +8877,9 @@ def indices(self, *mode, **kwargs): returned subspace. Note that if a multi-dimensional metadata construct is being used to define the indices - then some missing data may still be - inserted at unselected locations. + then some missing data may still need + to be inserted at unselected + locations. ``'envelope'`` The returned subspace is the smallest that contains all of the selected diff --git a/cf/mixin/fielddomain.py b/cf/mixin/fielddomain.py index 85cf503e55..3c0bc40ffb 100644 --- a/cf/mixin/fielddomain.py +++ b/cf/mixin/fielddomain.py @@ -250,7 +250,8 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): debug = is_log_level_debug(logger) compress = mode == "compress" - envelope = mode == "envelope" + halo_p1 = mode == "halo+1" + envelope = mode == "envelope" or halo_p1 full = mode == "full" domain_axes = self.domain_axes(todict=True) @@ -318,6 +319,7 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): ) mask = {} + cyclic = {} for canonical_axes, axes_key_construct_value_id in parsed.items(): axes, keys, constructs, points, identities = tuple( @@ -452,6 +454,8 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): ind = (d[index].array,) index = slice(None) + cyclic[axis] = True + elif item is not None: # 1-d CASE 3: All other 1-d cases if debug: @@ -529,7 +533,7 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): item_match &= m item_match = item_match.compute() - if np.ma.isMA: + if np.ma.isMA(item_match): ind = np.ma.where(item_match) else: ind = np.where(item_match) @@ -559,7 +563,7 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): # already selected for which the given value is in # fact outside of the cell. This could happen if the # cells are not rectangular (e.g. for curvilinear - # latitudes and longitudes array). + # latitudes and longitudes arrays). if n_items == constructs[0].ndim == len(bounds) == 2: point2 = [] for v, construct in zip(points, transposed_constructs): @@ -636,9 +640,11 @@ def _point_not_in_cell(nodes_x, nodes_y, point): if indices[axis] == slice(None): if compress: - # Create a compressed index for this axis + # Create a compressed index for this axis size = stop - start + 1 index = sorted(set(ind[i])) + if compress_p1: + elif envelope: # Create an envelope index for this axis stop += 1 @@ -652,7 +658,7 @@ def _point_not_in_cell(nodes_x, nodes_y, point): index = slice(None) else: raise ValueError( - "Must have mode full, envelope or compress" + "Must have mode full, envelope, or compress" ) # pragma: no cover indices[axis] = index @@ -662,19 +668,76 @@ def _point_not_in_cell(nodes_x, nodes_y, point): ind[i] -= start create_mask = ( - ancillary_mask + not halo_p1 + and ancillary_mask and data_axes and ind.shape[1] < masked_subspace_size ) else: create_mask = False + if compress_p1: + for axis, i in indices.items(): + if i == slice(None): + continue + + size = domain_axes[axis].get_size() + if isinstance(i, slice): + if cyclic.get(axis): + pass + else: + i = normalize_index(i, (size,))[0] + step = i.step + + i = np.arange(size)[i].tolist() + if step > 0: + stop = i[-1] + if stop <= size - n: + i.append(stop + 1) + + start = i[0] + if start > 0: + i.insert(0, start -1) + else: + + + i = normalize_index(i, (size,))[0].tolist() + start, stop, step = i.indices(size) + if step > 0: + if start >= step: + start -= step + + if stop + step <= size: + stop += step + else: + # step < 0 + if start - step < size: + start -= step + + if stop >= -step : + stop += step + else: + stop = None + + i = slice(start, stop , step) + elif np.iterable(i): + i = normalize_index(i, (size,))[0].tolist() + mx = i.max() + if mx < size: + i.append(mx + 1) + + mn = i.min() + if mn > 0: + i.insert(0, mn -1 ) + + indices[axis] = i + # Create an ancillary mask for these axes if debug: logger.debug( f" create_mask = {create_mask}" ) # pragma: no cover - + if create_mask: mask[canonical_axes] = _create_ancillary_mask_component( mask_component_shape, ind, compress From 5f774a68b4e4608b42f5c04929c9752c65f875e0 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Tue, 9 Apr 2024 23:54:32 +0100 Subject: [PATCH 02/30] dev --- cf/domain.py | 18 +++- cf/field.py | 154 ++++++++++++++++++--------------- cf/mixin/fielddomain.py | 185 ++++++++++++++++++++++++++-------------- 3 files changed, 219 insertions(+), 138 deletions(-) diff --git a/cf/domain.py b/cf/domain.py index 262dd56b7f..e35316418a 100644 --- a/cf/domain.py +++ b/cf/domain.py @@ -857,9 +857,21 @@ def indices(self, *mode, **kwargs): : time(1) = [2019-01-01 00:00:00] """ - if len(mode) > 1: + n_mode = len(mode) + if not n_mode: + halo = 0 + elif n_mode == 1: + try: + halo = int(mode[0]) + except ValueError: + halo = 0 + else: + mode = () + elif n_mode == 2: + halo = mode[1] + else: raise ValueError( - "Can't provide more than one positional argument. " + "Can't provide more than two positional argument. " f"Got: {', '.join(repr(x) for x in mode)}" ) @@ -872,7 +884,7 @@ def indices(self, *mode, **kwargs): # Get the indices for every domain axis in the domain, without # any auxiliary masks. - domain_indices = self._indices(mode, None, False, kwargs) + domain_indices = self._indices(mode, halo, None, False, kwargs) return domain_indices["indices"] diff --git a/cf/field.py b/cf/field.py index 56bdc3ec31..b358e5afad 100644 --- a/cf/field.py +++ b/cf/field.py @@ -9026,9 +9026,21 @@ def indices(self, *mode, **kwargs): removed_at="4.0.0", ) # pragma: no cover - if len(mode) > 1: + n_mode = len(mode) + if not n_mode: + halo = 0 + elif n_mode == 1: + try: + halo = int(mode[0]) + except ValueError: + halo = 0 + else: + mode = () + elif n_mode == 2: + halo = mode[1] + else: raise ValueError( - "Can't provide more than one positional argument. " + "Can't provide more than two positional arguments. " f"Got: {', '.join(repr(x) for x in mode)}" ) @@ -9039,13 +9051,15 @@ def indices(self, *mode, **kwargs): elif "full" in mode: mode = "full" else: - raise ValueError(f"Invalid value for 'mode' argument: {mode[0]!r}") + raise ValueError( + f"Invalid value for 'mode' positional argument: {mode[0]!r}" + ) data_axes = self.get_data_axes() # Get the indices for every domain axis in the domain, # including any ancillary masks - domain_indices = self._indices(mode, data_axes, True, kwargs) + domain_indices = self._indices(mode, halo, data_axes, True, kwargs) # Initialise the output indices with any ancillary masks. # Ensure that each ancillary mask is broadcastable to the @@ -9097,92 +9111,92 @@ def set_data( ): """Set the field construct data. - .. versionadded:: 3.0.0 + .. versionadded:: 3.0.0 - .. seealso:: `data`, `del_data`, `get_data`, `has_data`, - `set_construct` + .. seealso:: `data`, `del_data`, `get_data`, `has_data`, + `set_construct` - :Parameters: + :Parameters: - data: `Data` - The data to be inserted. + data: `Data` + The data to be inserted. - {{data_like}} + {{data_like}} - axes: (sequence of) `str` or `int`, optional - Set the domain axes constructs that are spanned by the - data. If unset, and the *set_axes* parameter is True, then - an attempt will be made to assign existing domain axis - constructs to the data. + axes: (sequence of) `str` or `int`, optional + Set the domain axes constructs that are spanned by the + data. If unset, and the *set_axes* parameter is True, then + an attempt will be made to assign existing domain axis + constructs to the data. - The contents of the *axes* parameter is mapped to domain - axis constructs by translating each element into a domain - axis construct key via the `domain_axis` method. + The contents of the *axes* parameter is mapped to domain + axis constructs by translating each element into a domain + axis construct key via the `domain_axis` method. - *Parameter example:* - ``axes='domainaxis1'`` + *Parameter example:* + ``axes='domainaxis1'`` - *Parameter example:* - ``axes='X'`` + *Parameter example:* + ``axes='X'`` - *Parameter example:* - ``axes=['latitude']`` + *Parameter example:* + ``axes=['latitude']`` - *Parameter example:* - ``axes=['X', 'longitude']`` + *Parameter example:* + ``axes=['X', 'longitude']`` - *Parameter example:* - ``axes=[1, 0]`` + *Parameter example:* + ``axes=[1, 0]`` - set_axes: `bool`, optional - If False then do not set the domain axes constructs that - are spanned by the data, even if the *axes* parameter has - been set. By default the axes are set either according to - the *axes* parameter, or if any domain axis constructs - exist then an attempt will be made to assign existing - domain axis constructs to the data. + set_axes: `bool`, optional + If False then do not set the domain axes constructs that + are spanned by the data, even if the *axes* parameter has + been set. By default the axes are set either according to + the *axes* parameter, or if any domain axis constructs + exist then an attempt will be made to assign existing + domain axis constructs to the data. - If the *axes* parameter is `None` and no domain axis - constructs exist then no attempt is made to assign domain - axes constructs to the data, regardless of the value of - *set_axes*. + If the *axes* parameter is `None` and no domain axis + constructs exist then no attempt is made to assign domain + axes constructs to the data, regardless of the value of + s *set_axes*. - copy: `bool`, optional - If True then set a copy of the data. By default the data - are copied. + copy: `bool`, optional + If True then set a copy of the data. By default the data + are copied. - {{inplace: `bool`, optional (default True)}} + {{inplace: `bool`, optional (default True)}} - .. versionadded:: 3.7.0 + .. versionadded:: 3.7.0 - :Returns: + :Returns: - `None` or `Field` - If the operation was in-place then `None` is returned, - otherwise return a new `Field` instance containing the new - data. + `None` or `Field` + If the operation was in-place then `None` is returned, + otherwise return a new `Field` instance containing the new + data. - **Examples** + **Examples** - >>> f = cf.Field() - >>> f.set_data([1, 2, 3]) - >>> f.has_data() - True - >>> f.get_data() - - >>> f.data - - >>> f.del_data() - - >>> g = f.set_data([4, 5, 6], inplace=False) - >>> g.data - - >>> f.has_data() - False - >>> print(f.get_data(None)) - None - >>> print(f.del_data(None)) - None + >>> f = cf.Field() + >>> f.set_data([1, 2, 3]) + >>> f.has_data() + True + >>> f.get_data() + + >>> f.data + + >>> f.del_data() + + >>> g = f.set_data([4, 5, 6], inplace=False) + >>> g.data + + >>> f.has_data() + False + >>> print(f.get_data(None)) + None + >>> print(f.del_data(None)) + None """ data = self._Data(data, copy=False) diff --git a/cf/mixin/fielddomain.py b/cf/mixin/fielddomain.py index 3c0bc40ffb..c0ea6260de 100644 --- a/cf/mixin/fielddomain.py +++ b/cf/mixin/fielddomain.py @@ -4,6 +4,7 @@ import dask.array as da import numpy as np from cfdm import is_log_level_debug, is_log_level_info +from dask.array.slicing import normalize_index from ..data import Data from ..decorators import ( @@ -198,7 +199,7 @@ def _equivalent_coordinate_references( return True - def _indices(self, mode, data_axes, ancillary_mask, kwargs): + def _indices(self, mode, halo, data_axes, ancillary_mask, kwargs): """Create indices that define a subspace of the field or domain construct. @@ -216,6 +217,13 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): The mode of operation. See the *mode* parameter of `indices` for details. + halo: `int` + The size of the halo to be added to the subspace, for + each axis specified by *kwargs*. The creation of + ancillary masks is automatically disabled when a + non-zero halo is provided, even when *ancillary_mask* + is True. + data_axes: sequence of `str`, or `None` The domain axis identifiers of the data axes, or `None` if there is no data array. @@ -250,10 +258,30 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): debug = is_log_level_debug(logger) compress = mode == "compress" - halo_p1 = mode == "halo+1" - envelope = mode == "envelope" or halo_p1 + envelope = mode == "envelope" full = mode == "full" + if halo: + try: + halo = int(halo) + except ValueError: + raise ValueError( + "'halo' positional argument must be convertible to an" + f"integer. Got {halo!r}" + ) + + if halo <0 : + raise ValueError( + "'halo' positional argument must be a non-negative " + f"integer. Got {halo!r}" + ) + + # Ancillary masks are automatically disabled when a + # non-zero halo is provided + ancillary_mask = False +# if full: +# raise ValueError("Can't halo full TODO") + domain_axes = self.domain_axes(todict=True) # Initialise the index for each axis @@ -319,7 +347,7 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): ) mask = {} - cyclic = {} + cyclic_index = {} for canonical_axes, axes_key_construct_value_id in parsed.items(): axes, keys, constructs, points, identities = tuple( @@ -447,6 +475,7 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): ) index = slice(start, stop, 1) + cyclic_index[axis] = index if full: d = self._Data(da.arange(size)) @@ -454,8 +483,6 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): ind = (d[index].array,) index = slice(None) - cyclic[axis] = True - elif item is not None: # 1-d CASE 3: All other 1-d cases if debug: @@ -640,11 +667,9 @@ def _point_not_in_cell(nodes_x, nodes_y, point): if indices[axis] == slice(None): if compress: - # Create a compressed index for this axis + # Create a compressed index for this axis size = stop - start + 1 index = sorted(set(ind[i])) - if compress_p1: - elif envelope: # Create an envelope index for this axis stop += 1 @@ -668,76 +693,106 @@ def _point_not_in_cell(nodes_x, nodes_y, point): ind[i] -= start create_mask = ( - not halo_p1 - and ancillary_mask + ancillary_mask and data_axes and ind.shape[1] < masked_subspace_size ) else: create_mask = False - if compress_p1: - for axis, i in indices.items(): - if i == slice(None): - continue - + # -------------------------------------------------------- + # TODO + # -------------------------------------------------------- + if halo: + # Note: For non-zero halos, 'ancillary_mask' is always + # False, and therefore 'create_mask' is also + # always False. This means that we don't need t + # worry about 'ind' becoming inconsistent with + # 'indices', because the former is only + # subsequently used if create_mask is True. + for axis in item_axes: size = domain_axes[axis].get_size() - if isinstance(i, slice): - if cyclic.get(axis): - pass + + index = cyclic_index.get(axis) + if index is not None: + # Cyclic index, that is a slice with step 1 or -1. + start = index.start + stop = index.stop + step = index.step + + if start < 0: + # Increasing cyclic slice + start = max(start - halo, stop - size) + stop = min(stop + halo, size + start) else: - i = normalize_index(i, (size,))[0] - step = i.step - - i = np.arange(size)[i].tolist() - if step > 0: - stop = i[-1] - if stop <= size - n: - i.append(stop + 1) - - start = i[0] - if start > 0: - i.insert(0, start -1) - else: - + # Decreasing cyclic slice + start = max(start + halo, size + stop) + stop = min(stop - halo, size - start) + + index = slice(start, stop, step) + else: + # Non-cyclic index, that is either a slice or + # a list/1-d array of int/bool. + index = normalize_index(indices[axis], (size,))[0] + if isinstance(index, slice): + index = np.arange(size)[index] + + index = index.tolist() + + # 'index' is now a list of positive integers + mn = min(index) + mx = max(index) + if first <= last: + index[:0] = range(max(0, first - halo), first) + index.extend(range(last + 1, max(last + 1 + halo, size-1)) + else: + index[:0] = range(min(first + halo, size-1), first, -1) + index.extend(range(last - 1, min(0, last - 1 - halo), -1)) + + if first <= last: + left = range(max(0, first - halo), first) + right = range(last + 1, max(last + 1 + halo, size-1)) + else: + left = range(min(first + halo, size-1), first, -1) + right = range(last - 1, min(0, last - 1 - halo), -1)) + + index[:0] = left + index.extend(right) + first = index[0] + last = index[-1] + if first <= last: + index[:0] = range(max(0, first - halo), first) + index.extend(range(last + 1, max(last + 1 + halo, size-1)) + else: + index[:0] = range(min(first + halo, size-1), first, -1) + index.extend(range(last - 1, min(0, last - 1 - halo), -1)) + + if first <= last: + left = range(max(0, first - halo), first) + right = range(last + 1, max(last + 1 + halo, size-1)) + else: + left = range(min(first + halo, size-1), first, -1) + right = range(last - 1, min(0, last - 1 - halo), -1)) + + index[:0] = left + index.extend(right) - i = normalize_index(i, (size,))[0].tolist() - start, stop, step = i.indices(size) - if step > 0: - if start >= step: - start -= step - - if stop + step <= size: - stop += step - else: - # step < 0 - if start - step < size: - start -= step - - if stop >= -step : - stop += step - else: - stop = None - - i = slice(start, stop , step) - elif np.iterable(i): - i = normalize_index(i, (size,))[0].tolist() - mx = i.max() - if mx < size: - i.append(mx + 1) - - mn = i.min() - if mn > 0: - i.insert(0, mn -1 ) - - indices[axis] = i - +# if not ( +# 0 <= index[0] < size and 0 <= index[-1] < size +# ): +# raise IndexError( +# f"Halo of size {halo} exceeds at least one " +# "array edge" +# ) + + indices[axis] = index + # Create an ancillary mask for these axes if debug: logger.debug( f" create_mask = {create_mask}" ) # pragma: no cover - + if create_mask: mask[canonical_axes] = _create_ancillary_mask_component( mask_component_shape, ind, compress From f61ccbd662ce933cc63a262b9b8981be4a66e9b3 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Wed, 10 Apr 2024 10:34:59 +0100 Subject: [PATCH 03/30] dev --- cf/functions.py | 191 +++++++++++++++++++++++++------------- cf/mixin/fielddomain.py | 121 ++++++++++++------------ cf/test/test_functions.py | 34 +++++++ 3 files changed, 223 insertions(+), 123 deletions(-) diff --git a/cf/functions.py b/cf/functions.py index 70eb6a81df..c06c402b56 100644 --- a/cf/functions.py +++ b/cf/functions.py @@ -1966,68 +1966,47 @@ def parse_indices(shape, indices, cyclic=False, keepdims=True): for i, (index, size) in enumerate(zip(parsed_indices, shape)): if cyclic and isinstance(index, slice): # Check for a cyclic slice - start = index.start - stop = index.stop - step = index.step - if start is None or stop is None: - step = 0 - elif step is None: - step = 1 - - if step > 0: - if 0 < start < size and 0 <= stop <= start: - # 6:0:1 => -4:0:1 - # 6:1:1 => -4:1:1 - # 6:3:1 => -4:3:1 - # 6:6:1 => -4:6:1 - start = size - start - elif -size <= start < 0 and -size <= stop <= start: - # -4:-10:1 => -4:1:1 - # -4:-9:1 => -4:1:1 - # -4:-7:1 => -4:3:1 - # -4:-4:1 => -4:6:1 - # -10:-10:1 => -10:0:1 - stop += size - elif step < 0: - if -size <= start < 0 and start <= stop < 0: - # -4:-1:-1 => 6:-1:-1 - # -4:-2:-1 => 6:-2:-1 - # -4:-4:-1 => 6:-4:-1 - # -10:-2:-1 => 0:-2:-1 - # -10:-10:-1 => 0:-10:-1 - start += size - elif 0 <= start < size and start < stop < size: - # 0:6:-1 => 0:-4:-1 - # 3:6:-1 => 3:-4:-1 - # 3:9:-1 => 3:-1:-1 - stop -= size - - if step > 0 and -size <= start < 0 and 0 <= stop <= size + start: - # x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] - # x[ -1:0:1] => [9] - # x[ -1:1:1] => [9, 0] - # x[ -1:3:1] => [9, 0, 1, 2] - # x[ -1:9:1] => [9, 0, 1, 2, 3, 4, 5, 6, 7, 8] - # x[ -4:0:1] => [6, 7, 8, 9] - # x[ -4:1:1] => [6, 7, 8, 9, 0] - # x[ -4:3:1] => [6, 7, 8, 9, 0, 1, 2] - # x[ -4:6:1] => [6, 7, 8, 9, 0, 1, 2, 3, 4, 5] - # x[ -9:0:1] => [1, 2, 3, 4, 5, 6, 7, 8, 9] - # x[ -9:1:1] => [1, 2, 3, 4, 5, 6, 7, 8, 9, 0] - # x[-10:0:1] => [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] - index = slice(0, stop - start, step) - roll[i] = -start - - elif step < 0 and 0 <= start < size and start - size <= stop < 0: - # x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] - # x[0: -4:-1] => [0, 9, 8, 7] - # x[6: -1:-1] => [6, 5, 4, 3, 2, 1, 0] - # x[6: -2:-1] => [6, 5, 4, 3, 2, 1, 0, 9] - # x[6: -4:-1] => [6, 5, 4, 3, 2, 1, 0, 9, 8, 7] - # x[0: -2:-1] => [0, 9] - # x[0:-10:-1] => [0, 9, 8, 7, 6, 5, 4, 3, 2, 1] - index = slice(start - stop - 1, None, step) - roll[i] = -1 - stop + try: + index = normalize_cyclic_slice(index, size) + except IndexError: + pass + else: + # Got a cyclic slice + start = index.start + stop = index.stop + step = index.step + if ( + step > 0 + and -size <= start < 0 + and 0 <= stop <= size + start + ): + # x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] + # x[ -1:0:1] => [9] + # x[ -1:1:1] => [9, 0] + # x[ -1:3:1] => [9, 0, 1, 2] + # x[ -1:9:1] => [9, 0, 1, 2, 3, 4, 5, 6, 7, 8] + # x[ -4:0:1] => [6, 7, 8, 9] + # x[ -4:1:1] => [6, 7, 8, 9, 0] + # x[ -4:3:1] => [6, 7, 8, 9, 0, 1, 2] + # x[ -4:6:1] => [6, 7, 8, 9, 0, 1, 2, 3, 4, 5] + # x[ -9:0:1] => [1, 2, 3, 4, 5, 6, 7, 8, 9] + # x[ -9:1:1] => [1, 2, 3, 4, 5, 6, 7, 8, 9, 0] + # x[-10:0:1] => [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] + index = slice(0, stop - start, step) + roll[i] = -start + + elif ( + step < 0 and 0 <= start < size and start - size <= stop < 0 + ): + # x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] + # x[0: -4:-1] => [0, 9, 8, 7] + # x[6: -1:-1] => [6, 5, 4, 3, 2, 1, 0] + # x[6: -2:-1] => [6, 5, 4, 3, 2, 1, 0, 9] + # x[6: -4:-1] => [6, 5, 4, 3, 2, 1, 0, 9, 8, 7] + # x[0: -2:-1] => [0, 9] + # x[0:-10:-1] => [0, 9, 8, 7, 6, 5, 4, 3, 2, 1] + index = slice(start - stop - 1, None, step) + roll[i] = -1 - stop elif keepdims and isinstance(index, Integral): # Convert an integral index to a slice @@ -2050,6 +2029,94 @@ def parse_indices(shape, indices, cyclic=False, keepdims=True): return parsed_indices, roll +def normalize_cyclic_slice(index, size): + """Normalise a cyclic slice. + + If the *index* is found to be a cyclic slice then a normalised + copy if it is returned. If a non-cyclic slice *index* is found + then an `IndexError` is raised. + + .. versionadded:: NEXTRELEASE + + :Parameters: + + index: + The index to be normalised. + + size: `int` + The size of the axis to which the *index* applies. + + :Returns: + + `slice` + The normalised cyclic slice. + + **Examples** + + >>> cf.normalize_cyclic_slice(slice(-2, 3), 8) + slice(-2, 3, 1) + >>> cf.normalize_cyclic_slice(slice(6, 3), 8) + slice(-2, 3, 1) + + >>> cf.normalize_cyclic_slice(slice(2, -3, -1), 8) + slice(2, -3, -1) + >>> cf.normalize_cyclic_slice(slice(2, 5, -1), 8) + slice(2, -3, -1) + + >>> cf.normalize_cyclic_slice(slice(1, 6), 8) + IndexError: slice(1, 6, None) is not a cyclic slice + >>> cf.normalize_cyclic_slice([1, 2, 3, 4], 8) + IndexError: [1, 2, 3, 4] is not a cyclic slice + + """ + if not isinstance(index, slice): + step = 0 + else: + start = index.start + stop = index.stop + step = index.step + if start is None or stop is None: + step = 0 + elif step is None: + step = 1 + + if step > 0: + if 0 < start < size and 0 <= stop <= start: + # 6:0:1 => -4:0:1 + # 6:1:1 => -4:1:1 + # 6:3:1 => -4:3:1 + # 6:6:1 => -4:6:1 + start -= size + elif -size <= start < 0 and -size <= stop <= start: + # -4:-10:1 => -4:1:1 + # -4:-9:1 => -4:1:1 + # -4:-7:1 => -4:3:1 + # -4:-4:1 => -4:6:1 + # -10:-10:1 => -10:0:1 + stop += size + elif step < 0: + if -size <= start < 0 and start <= stop < 0: + # -4:-1:-1 => 6:-1:-1 + # -4:-2:-1 => 6:-2:-1 + # -4:-4:-1 => 6:-4:-1 + # -10:-2:-1 => 0:-2:-1 + # -10:-10:-1 => 0:-10:-1 + start += size + elif 0 <= start < size and start < stop < size: + # 0:6:-1 => 0:-4:-1 + # 3:6:-1 => 3:-4:-1 + # 3:9:-1 => 3:-1:-1 + stop -= size + + if not ( + (step > 0 and start < 0 and stop > 0) + or (step < 0 and start > 0 and stop < 0) + ): + raise IndexError(f"{index!r} is not a cyclic slice") + + return slice(start, stop, step) + + def get_subspace(array, indices): """Return a subspace defined by the given indices of an array. diff --git a/cf/mixin/fielddomain.py b/cf/mixin/fielddomain.py index c0ea6260de..4cdc9d38d5 100644 --- a/cf/mixin/fielddomain.py +++ b/cf/mixin/fielddomain.py @@ -5,6 +5,7 @@ import numpy as np from cfdm import is_log_level_debug, is_log_level_info from dask.array.slicing import normalize_index +from dask.base import is_dask_collection from ..data import Data from ..decorators import ( @@ -13,7 +14,11 @@ _inplace_enabled_define_and_cleanup, _manage_log_level_via_verbosity, ) -from ..functions import _DEPRECATION_ERROR_KWARGS, bounds_combination_mode +from ..functions import ( + _DEPRECATION_ERROR_KWARGS, + bounds_combination_mode, + normalize_cyclic_slice, +) from ..query import Query from ..units import Units @@ -270,7 +275,7 @@ def _indices(self, mode, halo, data_axes, ancillary_mask, kwargs): f"integer. Got {halo!r}" ) - if halo <0 : + if halo < 0: raise ValueError( "'halo' positional argument must be a non-negative " f"integer. Got {halo!r}" @@ -279,8 +284,6 @@ def _indices(self, mode, halo, data_axes, ancillary_mask, kwargs): # Ancillary masks are automatically disabled when a # non-zero halo is provided ancillary_mask = False -# if full: -# raise ValueError("Can't halo full TODO") domain_axes = self.domain_axes(todict=True) @@ -474,6 +477,7 @@ def _indices(self, mode, halo, data_axes, ancillary_mask, kwargs): f"No indices found from: {identity}={value!r}" ) + print(start, stop) index = slice(start, stop, 1) cyclic_index[axis] = index @@ -710,80 +714,75 @@ def _point_not_in_cell(nodes_x, nodes_y, point): # worry about 'ind' becoming inconsistent with # 'indices', because the former is only # subsequently used if create_mask is True. + print("item_axes=", item_axes) + print(indices) for axis in item_axes: + index = indices[axis] size = domain_axes[axis].get_size() - index = cyclic_index.get(axis) - if index is not None: - # Cyclic index, that is a slice with step 1 or -1. - start = index.start - stop = index.stop - step = index.step - - if start < 0: - # Increasing cyclic slice - start = max(start - halo, stop - size) - stop = min(stop + halo, size + start) - else: - # Decreasing cyclic slice - start = max(start + halo, size + stop) - stop = min(stop - halo, size - start) - - index = slice(start, stop, step) - else: - # Non-cyclic index, that is either a slice or - # a list/1-d array of int/bool. - index = normalize_index(indices[axis], (size,))[0] + try: + index = normalize_cyclic_slice(index, size) + except IndexError: + # Non-cyclic index that is either a slice or a + # list/1-d array of int/bool. + index = normalize_index((index,), (size,))[0] if isinstance(index, slice): + step = index.step + increasing = step is None or step > 0 index = np.arange(size)[index] + else: + if is_dask_collection(index): + # Convert dask to numpy, and bool to int. + index = np.asanyarray(index) + index = normalize_index((index,), (size,))[0] + + increasing = index[0] <= index[-1] + # Convert 'index' to a list of integers (which + # will all be non-negative) index = index.tolist() - # 'index' is now a list of positive integers + # Extend the list at each end mn = min(index) mx = max(index) - if first <= last: - index[:0] = range(max(0, first - halo), first) - index.extend(range(last + 1, max(last + 1 + halo, size-1)) - else: - index[:0] = range(min(first + halo, size-1), first, -1) - index.extend(range(last - 1, min(0, last - 1 - halo), -1)) - - if first <= last: - left = range(max(0, first - halo), first) - right = range(last + 1, max(last + 1 + halo, size-1)) + if increasing: + # Put smaller indices in the left hand + # halo, and larger ones in the right hand + # halo + left = range(max(*(0, mn - halo)), mn) + right = range(mx + 1, min(*(mx + 1 + halo, size))) else: - left = range(min(first + halo, size-1), first, -1) - right = range(last - 1, min(0, last - 1 - halo), -1)) + # Put larger indices in the left hand + # halo, and smaller ones in the right hand + # halo + left = range(min(*(size - 1, mx + halo)), mx, -1) + right = range( + mn - 1, max(*(mn - 1 - halo, -1)), -1 + ) index[:0] = left index.extend(right) - first = index[0] - last = index[-1] - if first <= last: - index[:0] = range(max(0, first - halo), first) - index.extend(range(last + 1, max(last + 1 + halo, size-1)) - else: - index[:0] = range(min(first + halo, size-1), first, -1) - index.extend(range(last - 1, min(0, last - 1 - halo), -1)) + else: + # Cyclic slice index + start = index.start + stop = index.stop + step = index.step + if abs(step) > 1: + raise ValueError( + "Can only add halos to cyclic slices with " + f"step 1 or -1. Got {index!r}" + ) - if first <= last: - left = range(max(0, first - halo), first) - right = range(last + 1, max(last + 1 + halo, size-1)) + if step > 0: + # Increasing cyclic slice + start = max(start - halo, stop - size) + stop = min(stop + halo, size + start) else: - left = range(min(first + halo, size-1), first, -1) - right = range(last - 1, min(0, last - 1 - halo), -1)) + # Decreasing cyclic slice + start = min(start + halo, size + stop) + stop = max(stop - halo, start - size) - index[:0] = left - index.extend(right) - -# if not ( -# 0 <= index[0] < size and 0 <= index[-1] < size -# ): -# raise IndexError( -# f"Halo of size {halo} exceeds at least one " -# "array edge" -# ) + index = slice(start, stop, step) indices[axis] = index diff --git a/cf/test/test_functions.py b/cf/test/test_functions.py index a5e583aaaa..74596d0b93 100644 --- a/cf/test/test_functions.py +++ b/cf/test/test_functions.py @@ -354,6 +354,40 @@ def test_size(self): def test_CFA(self): self.assertEqual(cf.CFA(), cf.__cfa_version__) + def test_normalize_cyclic_slice(self): + # Cyclic slices + self.assertEqual( + cf.normalize_cyclic_slice(slice(-2, 3), 8), slice(-2, 3, 1) + ) + self.assertEqual( + cf.normalize_cyclic_slice(slice(6, 3), 8), slice(-2, 3, 1) + ) + self.assertEqual( + cf.normalize_cyclic_slice(slice(6, 3, 2), 8), slice(-2, 3, 2) + ) + + self.assertEqual( + cf.normalize_cyclic_slice(slice(2, -3, -1), 8), slice(2, -3, -1) + ) + self.assertEqual( + cf.normalize_cyclic_slice(slice(2, 5, -1), 8), slice(2, -3, -1) + ) + self.assertEqual( + cf.normalize_cyclic_slice(slice(2, 5, -2), 8), slice(2, -3, -2) + ) + + # Non-cylic slices + for index in ( + slice(1, 6), + slice(6, 1, -1), + slice(None, 4, None), + slice(1, 6, 0), + [1, 2], + 5, + ): + with self.assertRaises(IndexError): + cf.normalize_cyclic_slice(index, 8) + if __name__ == "__main__": print("Run date:", datetime.datetime.now()) From 01482b1766ef03d1ed4287573b6e0106b2961706 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Wed, 10 Apr 2024 16:03:11 +0100 Subject: [PATCH 04/30] dev --- cf/docstring/docstring.py | 124 +++++++++++++++++++ cf/domain.py | 136 +++++++++++++++------ cf/field.py | 88 +++++--------- cf/mixin/fielddomain.py | 82 ++++++------- cf/test/test_Field.py | 85 +++++++++++++ docs/attribute/cf.Domain.indices.html | 165 -------------------------- docs/source/class/cf.Domain.rst | 8 +- docs/source/class/cf.Field.rst | 14 ++- 8 files changed, 395 insertions(+), 307 deletions(-) delete mode 100644 docs/attribute/cf.Domain.indices.html diff --git a/cf/docstring/docstring.py b/cf/docstring/docstring.py index 94d2c03935..277e9e00f4 100644 --- a/cf/docstring/docstring.py +++ b/cf/docstring/docstring.py @@ -613,6 +613,76 @@ "{{to_size: `int`, optional}}": """to_size: `int`, optional Pad the axis after so that the new axis has the given size.""", + # indices mode options + "{{indices mode options}}": """TODO + Specify the mode of operation (``mode``) and a halo to + be added to the subspaced axes (``halo``) with + positional arguments in format ``mode``, ``halo``, or + ``mode, halo``, or with no positional arguments at + all. + + A mode is given as a `str`, and a halo as a + non-negative `int` (or any object that can be + converted to one). + + ============== ====================================== + *mode* Description + ============== ====================================== + ```` If no positional arguments are + provided then assume the + ``'compress'`` basic mode of operation + with no halo added to the subspaced + axes. Note that this is equivalent to + ``'compress', 0``. + + ``mode`` Define the basic mode of operation + with no halo added to the subspaced + axes. One of: + + * ``'compress'``: Unselected locations + are removed to create the returned + subspace. Note that if a + multi-dimensional metadata construct + is being used to define the indices + then some missing data may still + need to be inserted at unselected + locations. + + * ``'envelope'``: The returned + subspace is the smallest subsapce + that contains all of the selected + indices. Missing data is inserted at + unselected locations within the + envelope. + + * ``'full'`` The returned subspace has + the same domain as the original + field construct. Missing data is + inserted at unselected locations. + + * ``'test'``: May be used on its own + or in addition to one of the other + positional arguments. Do not create + a subspace, but return `True` or + `False` depending on whether or not + it is possible to create specified + the subspace. + + Note that ``mode`` is equivalent to + ``mode, 0``. + + ``mode, halo`` Define a basic mode of operation (one + of ``'compress'``, ``'envelope'``, or + ``'full'``) and define a halo to be + added to the subspaced axes. Note that + ``mode, 0`` is equivalent to ``mode``. + + ``halo`` Assume the ``'compress'`` basic mode + of operation and define a halo to be + added to the subspaced axes. Note + that ``halo`` is equivalent to + ``'compress', halo``. + ============== ======================================""", # ---------------------------------------------------------------- # Method description substitutions (4 levels of indentation) # ---------------------------------------------------------------- @@ -663,4 +733,58 @@ The removed CFA-netCDF file name substitution. If the substitution was not defined then an empty dictionary is returned.""", + # indices mode options + "{{indices mode options}}": """ + A mode of operation is given as a `str`, and a halo as + a non-negative `int` (or any object that can be + converted to one): + + ============== ====================================== + *mode* Description + ============== ====================================== + ```` If no positional arguments are + provided then assume the + ``'compress'`` basic mode of operation + with no halo added to the subspaced + axes. Note that this is equivalent to + ``'compress', 0``. + + ``mode`` Define the basic mode of operation + with no halo added to the subspaced + axes. One of: + + * ``'compress'``: Unselected locations + are removed to create the subspace. + Note that if a multi-dimensional + metadata construct is being used to + define the indices then some missing + data may still need to be inserted + at unselected locations. + + * ``'envelope'``: The subspace is the + smallest subspace that contains all + of the selected locations. Missing + data is inserted at unselected + locations within the envelope. + + * ``'full'`` The subspace has the same + domain as the original + construct. Missing data is inserted + at unselected locations. + + Note that ``mode`` is equivalent to + ``mode, 0``. + + ``mode, halo`` Define a basic mode of operation (one + of ``'compress'``, ``'envelope'``, or + ``'full'``) and define a halo to be + added to the subspaced axes. Note that + ``mode, 0`` is equivalent to ``mode``. + + ``halo`` Assume the ``'compress'`` basic mode + of operation and define a halo to be + added to the subspaced axes. Note + that ``halo`` is equivalent to + ``'compress', halo``. + ============== ======================================""", } diff --git a/cf/domain.py b/cf/domain.py index e35316418a..67d77737c0 100644 --- a/cf/domain.py +++ b/cf/domain.py @@ -785,24 +785,57 @@ def indices(self, *mode, **kwargs): :Parameters: - mode: `str`, *optional* - There are two modes of operation, each of which provides - indices for a different type of subspace: + mode: optional + Specify the mode of operation (``mode``) and a halo to + be added to the subspaced axes (``halo``) with + positional arguments in format ``mode``, ``halo``, or + ``mode, halo``, or with no positional arguments at + all. + + A mode of operation is given as a `str`, and a halo as + a non-negative `int` (or any object that can be + converted to one): ============== ====================================== *mode* Description ============== ====================================== - ``'compress'`` This is the default mode. Unselected - locations are removed to create the - returned subspace. Note that if a - multi-dimensional metadata construct - is being used to define the indices - then some missing data may still be - inserted at unselected locations. - - ``'envelope'`` The returned subspace is the smallest - that contains all of the selected - indices. + ```` If no positional arguments are + provided then assume the + ``'compress'`` basic mode of operation + with no halo added to the subspaced + axes. Note that this is equivalent to + ``'compress', 0``. + + ``mode`` Define the basic mode of operation + with no halo added to the subspaced + axes. One of: + + * ``'compress'``: Unselected locations + are removed to create the returned + subspace. Note that if a + multi-dimensional metadata construct + is being used to define the indices + then some unrequested locations may + also be selected. + + * ``'envelope'``: The subspace is the + smallest subspace that contains all + of the selected locations. + + Note that ``mode`` is equivalent to + ``mode, 0``. + + ``mode, halo`` Define a basic mode of operation (one + of ``'compress'``, ``'envelope'``, or + ``'full'``) and define a halo to be + added to the subspaced axes. Note that + ``mode, 0`` is equivalent to ``mode``. + + ``halo`` Assume the ``'compress'`` basic mode + of operation and define a halo to be + added to the subspaced axes. Note + that ``halo`` is equivalent to + ``'compress', halo``. ============== ====================================== kwargs: *optional* @@ -1174,33 +1207,58 @@ def subspace(self, *mode, **kwargs): :Parameters: - mode: `str`, *optional* - There are two modes of operation, each of which provides - indices for a different type of subspace: + mode: optional + Specify the mode of operation (``mode``) and a halo to + be added to the subspaced axes (``halo``) with + positional arguments in format ``mode``, ``halo``, or + ``mode, halo``, or with no positional arguments at + all. + + A mode of operation is given as a `str`, and a halo as + a non-negative `int` (or any object that can be + converted to one): - ============== ========================================== + ============== ====================================== *mode* Description - ============== ========================================== - ``'compress'`` Return indices that identify only the - requested locations. - - This is the default mode. - - Note that if a multi-dimensional metadata - construct is being used to define the - indices then some unrequested locations - may also be selected. - - ``'envelope'`` The returned subspace is the smallest that - contains all of the requested locations. - - ``'test'`` May be used on its own or in addition to - one of the other positional arguments. Do - not create a subspace, but return `True` - or `False` depending on whether or not it - is possible to create the specified - subspace. - ============== ========================================== + ============== ====================================== + ```` If no positional arguments are + provided then assume the + ``'compress'`` basic mode of operation + with no halo added to the subspaced + axes. Note that this is equivalent to + ``'compress', 0``. + + ``mode`` Define the basic mode of operation + with no halo added to the subspaced + axes. One of: + + * ``'compress'``: Unselected locations + are removed to create the returned + subspace. Note that if a + multi-dimensional metadata construct + is being used to define the indices + then some unrequested locations may + also be selected. + + * ``'envelope'``: The subspace is the + smallest subspace that contains all + of the selected locations. + + Note that ``mode`` is equivalent to + ``mode, 0``. + + ``mode, halo`` Define a basic mode of operation (one + of ``'compress'``, ``'envelope'``, or + ``'full'``) and define a halo to be + added to the subspaced axes. Note that + ``mode, 0`` is equivalent to ``mode``. + + ``halo`` Assume the ``'compress'`` basic mode + of operation and define a halo to be + added to the subspaced axes. Note + that ``halo`` is equivalent to + ``'compress', halo``. + ============== ====================================== kwargs: *optional* A keyword name is an identity of a metadata construct, and diff --git a/cf/field.py b/cf/field.py index b358e5afad..5e7349f9a6 100644 --- a/cf/field.py +++ b/cf/field.py @@ -8860,38 +8860,28 @@ def indices(self, *mode, **kwargs): actually created, these masks are all automatically applied to the result. + **Halos** + + A halo on either "side" of an axis will be automatically + reduced if including the complete halo would extend the + subspace beyond the axis's extent. + + If a non-zero halo has been defined then no ancillary masks + will be created. + .. seealso:: `subspace`, `where`, `__getitem__`, `__setitem__`, `cf.Domain.indices` :Parameters: - mode: `str`, *optional* - There are three modes of operation, each of which - provides indices for a different type of subspace: - - ============== ====================================== - *mode* Description - ============== ====================================== - ``'compress'`` This is the default mode. Unselected - locations are removed to create the - returned subspace. Note that if a - multi-dimensional metadata construct - is being used to define the indices - then some missing data may still need - to be inserted at unselected - locations. - - ``'envelope'`` The returned subspace is the smallest - that contains all of the selected - indices. Missing data is inserted at - unselected locations within the - envelope. - - ``'full'`` The returned subspace has the same - domain as the original field - construct. Missing data is inserted at - unselected locations. - ============== ====================================== + mode: optional + Specify the mode of operation (``mode``) and a halo to + be added to the subspaced axes (``halo``) with + positional arguments in format ``mode``, ``halo``, or + ``mode, halo``, or with no positional arguments at + all. + + {{indices mode options}} kwargs: *optional* A keyword name is an identity of a metadata construct, @@ -13307,42 +13297,24 @@ def subspace(self): ``3:-2:-1``) is assumed to "wrap" around, rather then producing a null result. - .. seealso:: `indices`, `squeeze`, `where`, `__getitem__` :Parameters: - positional arguments: *optional* - There are three modes of operation, each of which provides - a different type of subspace: + mode: optional + Specify the mode of operation (``mode``) and a halo to + be added to the subspaced axes (``halo``) with + positional arguments in format ``mode``, ``halo``, or + ``mode, halo``, or with no positional arguments at + all. - ============== ========================================== - *argument* Description - ============== ========================================== - ``'compress'`` This is the default mode. Unselected - locations are removed to create the - returned subspace. Note that if a - multi-dimensional metadata construct is - being used to define the indices then some - missing data may still be inserted at - unselected locations. - - ``'envelope'`` The returned subspace is the smallest that - contains all of the selected - indices. Missing data is inserted at - unselected locations within the envelope. - - ``'full'`` The returned subspace has the same domain - as the original field construct. Missing - data is inserted at unselected locations. - - ``'test'`` May be used on its own or in addition to - one of the other positional arguments. Do - not create a subspace, but return `True` - or `False` depending on whether or not it - is possible to create specified the - subspace. - ============== ========================================== + In addition, an extra postional argument of ``'test'`` + is allowed. If provided, then the subspace is not + returned , rather `True` or `False` is returned + depending on whether or not it is possible for the + subspace to be created. + + {{indices mode options}} keyword parameters: *optional* A keyword name is an identity of a metadata construct, and diff --git a/cf/mixin/fielddomain.py b/cf/mixin/fielddomain.py index 4cdc9d38d5..88cb60a325 100644 --- a/cf/mixin/fielddomain.py +++ b/cf/mixin/fielddomain.py @@ -223,11 +223,8 @@ def _indices(self, mode, halo, data_axes, ancillary_mask, kwargs): `indices` for details. halo: `int` - The size of the halo to be added to the subspace, for - each axis specified by *kwargs*. The creation of - ancillary masks is automatically disabled when a - non-zero halo is provided, even when *ancillary_mask* - is True. + The halo to be added to the subpsaced axes. See the + *mode* parameter of `indices` for details. data_axes: sequence of `str`, or `None` The domain axis identifiers of the data axes, or @@ -266,25 +263,21 @@ def _indices(self, mode, halo, data_axes, ancillary_mask, kwargs): envelope = mode == "envelope" full = mode == "full" + # Parse halo if halo: try: halo = int(halo) except ValueError: - raise ValueError( - "'halo' positional argument must be convertible to an" - f"integer. Got {halo!r}" - ) + ok = False + else: + ok = halo >= 0 - if halo < 0: + if not ok: raise ValueError( - "'halo' positional argument must be a non-negative " - f"integer. Got {halo!r}" + "halo positional argument must be convertible to a " + f"non-negative integer. Got {halo!r}" ) - # Ancillary masks are automatically disabled when a - # non-zero halo is provided - ancillary_mask = False - domain_axes = self.domain_axes(todict=True) # Initialise the index for each axis @@ -350,7 +343,6 @@ def _indices(self, mode, halo, data_axes, ancillary_mask, kwargs): ) mask = {} - cyclic_index = {} for canonical_axes, axes_key_construct_value_id in parsed.items(): axes, keys, constructs, points, identities = tuple( @@ -477,9 +469,7 @@ def _indices(self, mode, halo, data_axes, ancillary_mask, kwargs): f"No indices found from: {identity}={value!r}" ) - print(start, stop) index = slice(start, stop, 1) - cyclic_index[axis] = index if full: d = self._Data(da.arange(size)) @@ -708,14 +698,6 @@ def _point_not_in_cell(nodes_x, nodes_y, point): # TODO # -------------------------------------------------------- if halo: - # Note: For non-zero halos, 'ancillary_mask' is always - # False, and therefore 'create_mask' is also - # always False. This means that we don't need t - # worry about 'ind' becoming inconsistent with - # 'indices', because the former is only - # subsequently used if create_mask is True. - print("item_axes=", item_axes) - print(indices) for axis in item_axes: index = indices[axis] size = domain_axes[axis].get_size() @@ -725,6 +707,16 @@ def _point_not_in_cell(nodes_x, nodes_y, point): except IndexError: # Non-cyclic index that is either a slice or a # list/1-d array of int/bool. + # + # E.g. for halo=1 and size=5: + # slice(1, 3) -> [0, 1, 2, 3] + # slice(1, 4, 2) -> [0, 1, 3, 4] + # slice(2, 0, -1) -> [3, 2, 1, 0] + # [1, 2] -> [0, 1, 2, 3] + # [1, 3] -> [0, 1, 3, 4] + # [2, 1] -> [3, 2, 1, 0] + # [3, 1] -> [4, 3, 1, 0] + # [False, True, False, True, False] -> [0, 1, 3, 4] index = normalize_index((index,), (size,))[0] if isinstance(index, slice): step = index.step @@ -738,23 +730,23 @@ def _point_not_in_cell(nodes_x, nodes_y, point): increasing = index[0] <= index[-1] - # Convert 'index' to a list of integers (which - # will all be non-negative) + # Convert 'index' to a list (which will + # exclusively contain non-negative integers) index = index.tolist() # Extend the list at each end mn = min(index) mx = max(index) if increasing: - # Put smaller indices in the left hand - # halo, and larger ones in the right hand - # halo + # "Increasing" sequence: Put smaller + # indices in the left hand halo, and + # larger ones in the right hand halo left = range(max(*(0, mn - halo)), mn) right = range(mx + 1, min(*(mx + 1 + halo, size))) else: - # Put larger indices in the left hand - # halo, and smaller ones in the right hand - # halo + # "Decreasing" sequence: Put larger + # indices in the left hand halo, and + # smaller ones in the right hand halo left = range(min(*(size - 1, mx + halo)), mx, -1) right = range( mn - 1, max(*(mn - 1 - halo, -1)), -1 @@ -764,13 +756,17 @@ def _point_not_in_cell(nodes_x, nodes_y, point): index.extend(right) else: # Cyclic slice index + # + # E.g. for halo=1 and size=5: + # slice(-1, 2) -> slice(-2, 3) + # slice(1, -2, -1) -> slice(2, -3, -1) start = index.start stop = index.stop step = index.step if abs(step) > 1: raise ValueError( - "Can only add halos to cyclic slices with " - f"step 1 or -1. Got {index!r}" + "A cyclic slice can only have halos if it " + f"has step 1 or -1. Got {index!r}" ) if step > 0: @@ -793,9 +789,15 @@ def _point_not_in_cell(nodes_x, nodes_y, point): ) # pragma: no cover if create_mask: - mask[canonical_axes] = _create_ancillary_mask_component( - mask_component_shape, ind, compress - ) + if halo: + logger.warn( + "Not applying ancillary masks to a subspace with " + "non-zero halos" + ) + else: + mask[canonical_axes] = _create_ancillary_mask_component( + mask_component_shape, ind, compress + ) indices = {"indices": indices, "mask": mask} diff --git a/cf/test/test_Field.py b/cf/test/test_Field.py index 13b36bf426..8cf38c0966 100644 --- a/cf/test/test_Field.py +++ b/cf/test/test_Field.py @@ -1561,6 +1561,91 @@ def test_Field_indices(self): self.assertEqual(g.construct("aux_x").array, 160) self.assertEqual(g.construct("aux_y").array, 3) + # Halos + x0 = f.dimension_coordinate("X").array + for index in (cf.wi(70, 200), slice(2, 6), [2, 3, 4, 5]): + indices = f.indices(0, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 4)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == [80, 120, 160, 200]).all()) + + indices = f.indices(1, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 6)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == [40, 80, 120, 160, 200, 240]).all()) + + indices = f.indices(999, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 9)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == x0).all()) + + # Halos: non-monotonic sequence + index = [2, 3, 4, 1] + indices = f.indices(0, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 4)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == [80, 120, 160, 40]).all()) + + indices = f.indices(1, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 6)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == [200, 80, 120, 160, 40, 0]).all()) + + indices = f.indices(2, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 7)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == [240, 200, 80, 120, 160, 40, 0]).all()) + + # Halos: cyclic slice increasing + index = cf.wi(-170, 40) + indices = f.indices(0, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 6)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == [-160, -120, -80, -40, 0, 40]).all()) + + indices = f.indices(1, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 8)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == [-200, -160, -120, -80, -40, 0, 40, 80]).all()) + + indices = f.indices(2, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 9)) + x = g.dimension_coordinate("X").array + self.assertTrue( + (x == [-240, -200, -160, -120, -80, -40, 0, 40, 80]).all() + ) + + # Halos: cyclic slice decreasing + index = slice(1, -5, -1) + indices = f.indices(0, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 6)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == [40, 0, -40, -80, -120, -160]).all()) + + indices = f.indices(1, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 8)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == [80, 40, 0, -40, -80, -120, -160, -200]).all()) + + indices = f.indices(2, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 9)) + x = g.dimension_coordinate("X").array + self.assertTrue( + (x == [120, 80, 40, 0, -40, -80, -120, -160, -200]).all() + ) + # Subspace has size 0 axis resulting from dask array index indices = f.indices(grid_latitude=cf.contains(-23.2)) with self.assertRaises(IndexError): diff --git a/docs/attribute/cf.Domain.indices.html b/docs/attribute/cf.Domain.indices.html deleted file mode 100644 index b0f069b511..0000000000 --- a/docs/attribute/cf.Domain.indices.html +++ /dev/null @@ -1,165 +0,0 @@ - - - - - - - - cf.Domain.indices — Documentation - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - -
-
- - -
- -
-

cf.Domain.indices

-
-
-Domain.indices = <function Domain.indices>[source]
-
- -
- - -
- -
-
-
-
- - - - - - - \ No newline at end of file diff --git a/docs/source/class/cf.Domain.rst b/docs/source/class/cf.Domain.rst index 526700e480..4d09b2b047 100644 --- a/docs/source/class/cf.Domain.rst +++ b/docs/source/class/cf.Domain.rst @@ -211,13 +211,19 @@ Domain axes Subspacing ---------- +.. autosummary:: + :nosignatures: + :toctree: ../method/ + :template: method.rst + + ~cf.Domain.indices + .. autosummary:: :nosignatures: :toctree: ../attribute/ :template: attribute.rst ~cf.Domain.subspace - ~cf.Domain.indices NetCDF ------ diff --git a/docs/source/class/cf.Field.rst b/docs/source/class/cf.Field.rst index f56f6eff61..c0a0275daa 100644 --- a/docs/source/class/cf.Field.rst +++ b/docs/source/class/cf.Field.rst @@ -169,7 +169,8 @@ Data ~cf.Field.ndim ~cf.Field.shape ~cf.Field.size - ~cf.Field.varray + ~cf.Field.to_dask_array + ~cf.Field.varray .. rubric:: *Rearranging elements* @@ -524,15 +525,20 @@ Domain axes Subspacing ---------- +.. autosummary:: + :nosignatures: + :toctree: ../method/ + :template: method.rst + + ~cf.Field.__getitem__ + ~cf.Field.indices + .. autosummary:: :nosignatures: :toctree: ../attribute/ :template: attribute.rst - ~cf.Field.to_dask_array - ~cf.Field.__getitem__ ~cf.Field.subspace - ~cf.Field.indices Mathematical operations ----------------------- From 3f5c3ec23cfeb33430d1f3a321ca853f16cc28ef Mon Sep 17 00:00:00 2001 From: David Hassell Date: Wed, 10 Apr 2024 20:44:53 +0100 Subject: [PATCH 05/30] dev --- cf/docstring/docstring.py | 173 +++++++++++++++----------------------- cf/domain.py | 116 +++---------------------- cf/field.py | 42 ++++----- 3 files changed, 104 insertions(+), 227 deletions(-) diff --git a/cf/docstring/docstring.py b/cf/docstring/docstring.py index 277e9e00f4..35716143d7 100644 --- a/cf/docstring/docstring.py +++ b/cf/docstring/docstring.py @@ -80,6 +80,19 @@ Whether `esmpy` logging is enabled or not is determined by `cf.regrid_logging`. If it is enabled then logging takes place after every call. By default logging is disabled.""", + # indices halos + "{{indices halos}}": """ + If a halo is defined via a positional argument, then the + returned index for each subspaced axis will be extended to + include that many extra elements at each "side" of the + axis. For instance, ``f.indices(X=slice(10, 20))`` will give + identical results to each of ``f.indices(0, X=slice(10, + 20))``, ``f.indices(1, X=slice(11, 19))``, ``f.indices(2, + X=slice(12, 18))``, etc. + + The number of extra elements will be automatically reduced if + including full amount defined by the halo would extend the + subspace beyond the axis's extent.""", # ---------------------------------------------------------------- # Method description substitutions (3 levels of indentation) # ---------------------------------------------------------------- @@ -613,76 +626,6 @@ "{{to_size: `int`, optional}}": """to_size: `int`, optional Pad the axis after so that the new axis has the given size.""", - # indices mode options - "{{indices mode options}}": """TODO - Specify the mode of operation (``mode``) and a halo to - be added to the subspaced axes (``halo``) with - positional arguments in format ``mode``, ``halo``, or - ``mode, halo``, or with no positional arguments at - all. - - A mode is given as a `str`, and a halo as a - non-negative `int` (or any object that can be - converted to one). - - ============== ====================================== - *mode* Description - ============== ====================================== - ```` If no positional arguments are - provided then assume the - ``'compress'`` basic mode of operation - with no halo added to the subspaced - axes. Note that this is equivalent to - ``'compress', 0``. - - ``mode`` Define the basic mode of operation - with no halo added to the subspaced - axes. One of: - - * ``'compress'``: Unselected locations - are removed to create the returned - subspace. Note that if a - multi-dimensional metadata construct - is being used to define the indices - then some missing data may still - need to be inserted at unselected - locations. - - * ``'envelope'``: The returned - subspace is the smallest subsapce - that contains all of the selected - indices. Missing data is inserted at - unselected locations within the - envelope. - - * ``'full'`` The returned subspace has - the same domain as the original - field construct. Missing data is - inserted at unselected locations. - - * ``'test'``: May be used on its own - or in addition to one of the other - positional arguments. Do not create - a subspace, but return `True` or - `False` depending on whether or not - it is possible to create specified - the subspace. - - Note that ``mode`` is equivalent to - ``mode, 0``. - - ``mode, halo`` Define a basic mode of operation (one - of ``'compress'``, ``'envelope'``, or - ``'full'``) and define a halo to be - added to the subspaced axes. Note that - ``mode, 0`` is equivalent to ``mode``. - - ``halo`` Assume the ``'compress'`` basic mode - of operation and define a halo to be - added to the subspaced axes. Note - that ``halo`` is equivalent to - ``'compress', halo``. - ============== ======================================""", # ---------------------------------------------------------------- # Method description substitutions (4 levels of indentation) # ---------------------------------------------------------------- @@ -734,7 +677,12 @@ substitution was not defined then an empty dictionary is returned.""", # indices mode options - "{{indices mode options}}": """ + "{{indices mode options}}": """Specify the mode of operation (``mode``) and a halo to + be added to the subspaced axes (``halo``) with + positional arguments in format ``mode``, ``halo``, or + ``mode, halo``, or with no positional arguments at + all. + A mode of operation is given as a `str`, and a halo as a non-negative `int` (or any object that can be converted to one): @@ -747,44 +695,61 @@ ``'compress'`` basic mode of operation with no halo added to the subspaced axes. Note that this is equivalent to - ``'compress', 0``. + either of ``'compress', 0`` and ``0``. ``mode`` Define the basic mode of operation with no halo added to the subspaced - axes. One of: - - * ``'compress'``: Unselected locations - are removed to create the subspace. - Note that if a multi-dimensional - metadata construct is being used to - define the indices then some missing - data may still need to be inserted - at unselected locations. - - * ``'envelope'``: The subspace is the - smallest subspace that contains all - of the selected locations. Missing - data is inserted at unselected - locations within the envelope. - - * ``'full'`` The subspace has the same - domain as the original - construct. Missing data is inserted - at unselected locations. - - Note that ``mode`` is equivalent to - ``mode, 0``. - - ``mode, halo`` Define a basic mode of operation (one - of ``'compress'``, ``'envelope'``, or - ``'full'``) and define a halo to be - added to the subspaced axes. Note that - ``mode, 0`` is equivalent to ``mode``. - - ``halo`` Assume the ``'compress'`` basic mode - of operation and define a halo to be + axes. Note that ``mode`` is equivalent + to ``mode, 0``. + + ``mode, halo`` Define a basic mode of operation, as + well as a halo to be added to the + subspaced axes. Note that ``mode, 0`` + is equivalent to ``mode``. + + ``halo`` Assume the default basic mode of + operation and define a halo to be added to the subspaced axes. Note that ``halo`` is equivalent to ``'compress', halo``. ============== ======================================""", + # indices valid modes Field + "{{indices valid modes Field}}": """ + Valid modes are: + + * ``'compress'`` This the default. + Unselected locations are removed to create the + subspace. Missing data may be inserted at + unselected locations, but only where it is + required to ensure a hyperrectangular result, and + not if a non-zero halo has been defined. + + * ``'envelope'`` + The subspace is the smallest subspace that + contains all of the selected locations. Missing + data is inserted at unselected locations within + the envelope, unless a non-zero halo has been + defined. + + * ``'full'`` + The subspace has the same domain as the original + construct. Missing data is inserted at unselected + locations, unless a non-zero halo has been + defined.""", + # indices valid modes Domain + "{{indices valid modes Domain}}": """ + Valid modes are: + + * ``'compress'`` This the default. + Unselected locations are removed to create the + subspace, unless they are required to ensure a + hyperrectangular result. + + * ``'envelope'`` + The subspace is the smallest subspace that + contains all of the selected locations. + + * ``'full'`` + The subspace has the same domain as the original + construct.""", } diff --git a/cf/domain.py b/cf/domain.py index 67d77737c0..34daaeb8d4 100644 --- a/cf/domain.py +++ b/cf/domain.py @@ -778,6 +778,10 @@ def indices(self, *mode, **kwargs): may still need to be inserted into the field construct's data. + **Halos** + + {{indices halos}} + .. versionadded:: 3.11.0 .. seealso:: `subspace`, `where`, `__getitem__`, @@ -786,57 +790,9 @@ def indices(self, *mode, **kwargs): :Parameters: mode: optional - Specify the mode of operation (``mode``) and a halo to - be added to the subspaced axes (``halo``) with - positional arguments in format ``mode``, ``halo``, or - ``mode, halo``, or with no positional arguments at - all. - - A mode of operation is given as a `str`, and a halo as - a non-negative `int` (or any object that can be - converted to one): - - ============== ====================================== - *mode* Description - ============== ====================================== - ```` If no positional arguments are - provided then assume the - ``'compress'`` basic mode of operation - with no halo added to the subspaced - axes. Note that this is equivalent to - ``'compress', 0``. - - ``mode`` Define the basic mode of operation - with no halo added to the subspaced - axes. One of: - - * ``'compress'``: Unselected locations - are removed to create the returned - subspace. Note that if a - multi-dimensional metadata construct - is being used to define the indices - then some unrequested locations may - also be selected. - - * ``'envelope'``: The subspace is the - smallest subspace that contains all - of the selected locations. - - Note that ``mode`` is equivalent to - ``mode, 0``. - - ``mode, halo`` Define a basic mode of operation (one - of ``'compress'``, ``'envelope'``, or - ``'full'``) and define a halo to be - added to the subspaced axes. Note that - ``mode, 0`` is equivalent to ``mode``. - - ``halo`` Assume the ``'compress'`` basic mode - of operation and define a halo to be - added to the subspaced axes. Note - that ``halo`` is equivalent to - ``'compress', halo``. - ============== ====================================== + {{indices mode options}} + + {{indices valid modes Domain}} kwargs: *optional* A keyword name is an identity of a metadata construct, @@ -904,7 +860,7 @@ def indices(self, *mode, **kwargs): halo = mode[1] else: raise ValueError( - "Can't provide more than two positional argument. " + "Can't provide more than two positional arguments. " f"Got: {', '.join(repr(x) for x in mode)}" ) @@ -1203,62 +1159,14 @@ def subspace(self, *mode, **kwargs): .. versionadded:: 3.11.0 - .. seealso:: `indices` + .. seealso:: `indices`, `cf.Field.subspace` :Parameters: mode: optional - Specify the mode of operation (``mode``) and a halo to - be added to the subspaced axes (``halo``) with - positional arguments in format ``mode``, ``halo``, or - ``mode, halo``, or with no positional arguments at - all. - - A mode of operation is given as a `str`, and a halo as - a non-negative `int` (or any object that can be - converted to one): - - ============== ====================================== - *mode* Description - ============== ====================================== - ```` If no positional arguments are - provided then assume the - ``'compress'`` basic mode of operation - with no halo added to the subspaced - axes. Note that this is equivalent to - ``'compress', 0``. - - ``mode`` Define the basic mode of operation - with no halo added to the subspaced - axes. One of: - - * ``'compress'``: Unselected locations - are removed to create the returned - subspace. Note that if a - multi-dimensional metadata construct - is being used to define the indices - then some unrequested locations may - also be selected. - - * ``'envelope'``: The subspace is the - smallest subspace that contains all - of the selected locations. - - Note that ``mode`` is equivalent to - ``mode, 0``. - - ``mode, halo`` Define a basic mode of operation (one - of ``'compress'``, ``'envelope'``, or - ``'full'``) and define a halo to be - added to the subspaced axes. Note that - ``mode, 0`` is equivalent to ``mode``. - - ``halo`` Assume the ``'compress'`` basic mode - of operation and define a halo to be - added to the subspaced axes. Note - that ``halo`` is equivalent to - ``'compress', halo``. - ============== ====================================== + {{indices mode options}} + + {{indices valid modes Domain}} kwargs: *optional* A keyword name is an identity of a metadata construct, and diff --git a/cf/field.py b/cf/field.py index 5e7349f9a6..cbb4d752c5 100644 --- a/cf/field.py +++ b/cf/field.py @@ -8862,27 +8862,23 @@ def indices(self, *mode, **kwargs): **Halos** - A halo on either "side" of an axis will be automatically - reduced if including the complete halo would extend the - subspace beyond the axis's extent. + {{indices halos}} If a non-zero halo has been defined then no ancillary masks will be created. + .. versionadded:: 1.0 + .. seealso:: `subspace`, `where`, `__getitem__`, `__setitem__`, `cf.Domain.indices` :Parameters: mode: optional - Specify the mode of operation (``mode``) and a halo to - be added to the subspaced axes (``halo``) with - positional arguments in format ``mode``, ``halo``, or - ``mode, halo``, or with no positional arguments at - all. - {{indices mode options}} + {{indices valid modes Field}} + kwargs: *optional* A keyword name is an identity of a metadata construct, and the keyword value provides a condition for @@ -13299,22 +13295,30 @@ def subspace(self): .. seealso:: `indices`, `squeeze`, `where`, `__getitem__` + **Halos** + + {{index halos}} + + .. note:: If a non-zero halo has been defined then no + ancillary masks will be created. + + .. versionadded:: 1.0 + + .. seealso:: `indices`, `where`, `__getitem__`, + `__setitem__`, `cf.Domain.subspace` + :Parameters: mode: optional - Specify the mode of operation (``mode``) and a halo to - be added to the subspaced axes (``halo``) with - positional arguments in format ``mode``, ``halo``, or - ``mode, halo``, or with no positional arguments at - all. + {{indices mode options}} + {{indices valid modes Field}} + In addition, an extra postional argument of ``'test'`` - is allowed. If provided, then the subspace is not - returned , rather `True` or `False` is returned + is allowed. When provided the subspace is not + returned, instead `True` or `False` is returned depending on whether or not it is possible for the - subspace to be created. - - {{indices mode options}} + requested subspace to be created. keyword parameters: *optional* A keyword name is an identity of a metadata construct, and From 8b52f4f83785600234f353867f51fada3c7207bd Mon Sep 17 00:00:00 2001 From: David Hassell Date: Thu, 11 Apr 2024 10:04:34 +0100 Subject: [PATCH 06/30] dev --- cf/docstring/docstring.py | 93 +++---- cf/domain.py | 63 ++--- cf/field.py | 44 +--- cf/mixin/fielddomain.py | 197 ++++++++++----- docs/attribute/cf.Domain.subspace.html | 165 ------------- docs/attribute/cf.Field.indices.html | 165 ------------- docs/method/cf.Field.subspace.html | 323 ------------------------- docs/source/class/cf.Domain.rst | 6 - 8 files changed, 215 insertions(+), 841 deletions(-) delete mode 100644 docs/attribute/cf.Domain.subspace.html delete mode 100644 docs/attribute/cf.Field.indices.html delete mode 100644 docs/method/cf.Field.subspace.html diff --git a/cf/docstring/docstring.py b/cf/docstring/docstring.py index 35716143d7..fd03a2ce8c 100644 --- a/cf/docstring/docstring.py +++ b/cf/docstring/docstring.py @@ -90,6 +90,18 @@ 20))``, ``f.indices(1, X=slice(11, 19))``, ``f.indices(2, X=slice(12, 18))``, etc. + The number of extra elements will be automatically reduced if + including full amount defined by the halo would extend the + subspace beyond the axis's extent.""", + # indices halos + "{{subspace halos}}": """ + If a halo is defined via a positional argument, then each + subspaced axis will be extended to include that many extra + elements at each "side" of the axis. For instance, + ``f.subspace(X=slice(10, 20))`` will give identical results to + each of ``f.subspace(0, X=slice(10, 20))``, ``f.subspace(1, + X=slice(11, 19))``, ``f.subspace(2, X=slice(12, 18))``, etc. + The number of extra elements will be automatically reduced if including full amount defined by the halo would extend the subspace beyond the axis's extent.""", @@ -626,6 +638,44 @@ "{{to_size: `int`, optional}}": """to_size: `int`, optional Pad the axis after so that the new axis has the given size.""", + # indices mode options + "{{mode: optional}}": """mode: optional + Specify the mode of operation (``mode``) and a halo to + be added to the subspaced axes (``halo``) with + positional arguments in format ``mode``, or ``halo``, + or ``mode, halo``, or with no positional arguments at + all. + + A mode of operation is given as a `str`, and a halo as + a non-negative `int` (or any object that can be + converted to one): + + ============== ====================================== + *mode* Description + ============== ====================================== + ```` If no positional arguments are + provided then assume the + ``'compress'`` mode of operation with + no halo added to the subspaced + axes. Note that this is equivalent to + either of ``'compress', 0`` and ``0``. + + ``mode`` Define the mode of operation with no + halo added to the subspaced axes. Note + that ``mode`` is equivalent to ``mode, + 0``. + + ``mode, halo`` Define a mode of operation, as well as + a halo to be added to the subspaced + axes. Note that ``mode, 0`` is + equivalent to ``mode``. + + ``halo`` Assume the ``'compress'`` mode of + operation and define a halo to be + added to the subspaced axes. Note + that ``halo`` is equivalent to + ``'compress', halo``. + ============== ======================================""", # ---------------------------------------------------------------- # Method description substitutions (4 levels of indentation) # ---------------------------------------------------------------- @@ -676,43 +726,6 @@ The removed CFA-netCDF file name substitution. If the substitution was not defined then an empty dictionary is returned.""", - # indices mode options - "{{indices mode options}}": """Specify the mode of operation (``mode``) and a halo to - be added to the subspaced axes (``halo``) with - positional arguments in format ``mode``, ``halo``, or - ``mode, halo``, or with no positional arguments at - all. - - A mode of operation is given as a `str`, and a halo as - a non-negative `int` (or any object that can be - converted to one): - - ============== ====================================== - *mode* Description - ============== ====================================== - ```` If no positional arguments are - provided then assume the - ``'compress'`` basic mode of operation - with no halo added to the subspaced - axes. Note that this is equivalent to - either of ``'compress', 0`` and ``0``. - - ``mode`` Define the basic mode of operation - with no halo added to the subspaced - axes. Note that ``mode`` is equivalent - to ``mode, 0``. - - ``mode, halo`` Define a basic mode of operation, as - well as a halo to be added to the - subspaced axes. Note that ``mode, 0`` - is equivalent to ``mode``. - - ``halo`` Assume the default basic mode of - operation and define a halo to be - added to the subspaced axes. Note - that ``halo`` is equivalent to - ``'compress', halo``. - ============== ======================================""", # indices valid modes Field "{{indices valid modes Field}}": """ Valid modes are: @@ -747,9 +760,5 @@ * ``'envelope'`` The subspace is the smallest subspace that - contains all of the selected locations. - - * ``'full'`` - The subspace has the same domain as the original - construct.""", + contains all of the selected locations.""", } diff --git a/cf/domain.py b/cf/domain.py index 34daaeb8d4..d1243fcd50 100644 --- a/cf/domain.py +++ b/cf/domain.py @@ -743,7 +743,7 @@ def indices(self, *mode, **kwargs): """Create indices that define a subspace of the domain construct. - The indices returned by this method be used to create the + The indices returned by this method may be used to create the subspace by passing them to the `subspace` method of the original domain construct. @@ -789,9 +789,7 @@ def indices(self, *mode, **kwargs): :Parameters: - mode: optional - {{indices mode options}} - + {{mode: optional}} {{indices valid modes Domain}} kwargs: *optional* @@ -846,34 +844,9 @@ def indices(self, *mode, **kwargs): : time(1) = [2019-01-01 00:00:00] """ - n_mode = len(mode) - if not n_mode: - halo = 0 - elif n_mode == 1: - try: - halo = int(mode[0]) - except ValueError: - halo = 0 - else: - mode = () - elif n_mode == 2: - halo = mode[1] - else: - raise ValueError( - "Can't provide more than two positional arguments. " - f"Got: {', '.join(repr(x) for x in mode)}" - ) - - if not mode or "compress" in mode: - mode = "compress" - elif "envelope" in mode: - mode = "envelope" - else: - raise ValueError(f"Invalid value for 'mode' argument: {mode[0]!r}") - # Get the indices for every domain axis in the domain, without # any auxiliary masks. - domain_indices = self._indices(mode, halo, None, False, kwargs) + domain_indices = self._indices(mode, None, False, kwargs) return domain_indices["indices"] @@ -1121,19 +1094,19 @@ def roll(self, axis, shift, inplace=False): return d def subspace(self, *mode, **kwargs): - """Create indices that define a subspace of the domain - construct. + """Create a subspace of the field construct. - The indices returned by this method be used to create the subspace - by passing them to the `subspace` method of the original domain - construct. + Creation of a new domain construct which spans a subspace of + the domain of an existing domain construct is achieved by + identifying indices based on the metadata constructs + (subspacing by metadata). The new domain construct is created + with the same properties as the original domain construct. - The subspace is defined by identifying indices based on the - metadata constructs. + **Subspacing by metadata** - Metadata constructs are selected conditions are specified on their - data. Indices for subspacing are then automatically inferred from - where the conditions are met. + Subspacing by metadata selects metadata constructs and + specifies conditions on their data. Indices for subspacing are + then automatically inferred from where the conditions are met. Metadata constructs and the conditions on their data are defined by keyword parameters. @@ -1143,6 +1116,9 @@ def subspace(self, *mode, **kwargs): * Multiple domain axes may be subspaced simultaneously, and it doesn't matter which order they are specified in. + * Subspace criteria may be provided for size 1 domain axes that + are not spanned by the field construct's data. + * Explicit indices may also be assigned to a domain axis identified by a metadata construct, with either a Python `slice` object, or a sequence of integers or booleans. @@ -1157,14 +1133,17 @@ def subspace(self, *mode, **kwargs): acting along orthogonal dimensions, some missing data may still need to be inserted into the field construct's data. + **Halos** + + {{subspace halos}} + .. versionadded:: 3.11.0 .. seealso:: `indices`, `cf.Field.subspace` :Parameters: - mode: optional - {{indices mode options}} + {{mode: optional}} {{indices valid modes Domain}} diff --git a/cf/field.py b/cf/field.py index cbb4d752c5..eab4169ef8 100644 --- a/cf/field.py +++ b/cf/field.py @@ -8874,8 +8874,7 @@ def indices(self, *mode, **kwargs): :Parameters: - mode: optional - {{indices mode options}} + {{mode: optional}} {{indices valid modes Field}} @@ -9012,40 +9011,10 @@ def indices(self, *mode, **kwargs): removed_at="4.0.0", ) # pragma: no cover - n_mode = len(mode) - if not n_mode: - halo = 0 - elif n_mode == 1: - try: - halo = int(mode[0]) - except ValueError: - halo = 0 - else: - mode = () - elif n_mode == 2: - halo = mode[1] - else: - raise ValueError( - "Can't provide more than two positional arguments. " - f"Got: {', '.join(repr(x) for x in mode)}" - ) - - if not mode or "compress" in mode: - mode = "compress" - elif "envelope" in mode: - mode = "envelope" - elif "full" in mode: - mode = "full" - else: - raise ValueError( - f"Invalid value for 'mode' positional argument: {mode[0]!r}" - ) - - data_axes = self.get_data_axes() - # Get the indices for every domain axis in the domain, # including any ancillary masks - domain_indices = self._indices(mode, halo, data_axes, True, kwargs) + data_axes = self.get_data_axes() + domain_indices = self._indices(mode, data_axes, True, kwargs) # Initialise the output indices with any ancillary masks. # Ensure that each ancillary mask is broadcastable to the @@ -13297,7 +13266,7 @@ def subspace(self): **Halos** - {{index halos}} + {{subspace halos}} .. note:: If a non-zero halo has been defined then no ancillary masks will be created. @@ -13309,11 +13278,10 @@ def subspace(self): :Parameters: - mode: optional - {{indices mode options}} + {{mode: optional}} {{indices valid modes Field}} - + In addition, an extra postional argument of ``'test'`` is allowed. When provided the subspace is not returned, instead `True` or `False` is returned diff --git a/cf/mixin/fielddomain.py b/cf/mixin/fielddomain.py index 88cb60a325..616ee555bc 100644 --- a/cf/mixin/fielddomain.py +++ b/cf/mixin/fielddomain.py @@ -204,7 +204,7 @@ def _equivalent_coordinate_references( return True - def _indices(self, mode, halo, data_axes, ancillary_mask, kwargs): + def _indices(self, mode, data_axes, ancillary_mask, kwargs): """Create indices that define a subspace of the field or domain construct. @@ -218,12 +218,8 @@ def _indices(self, mode, halo, data_axes, ancillary_mask, kwargs): :Parameters: - mode: `str` - The mode of operation. See the *mode* parameter of - `indices` for details. - - halo: `int` - The halo to be added to the subpsaced axes. See the + mode: `tuple` + The mode of operation and the halo size . See the *mode* parameter of `indices` for details. data_axes: sequence of `str`, or `None` @@ -259,11 +255,33 @@ def _indices(self, mode, halo, data_axes, ancillary_mask, kwargs): """ debug = is_log_level_debug(logger) - compress = mode == "compress" + # Parse mode + n_mode = len(mode) + if not n_mode: + mode = None + halo = 0 + elif n_mode == 1: + try: + halo = int(mode[0]) + except ValueError: + mode = mode[0] + halo = 0 + else: + mode = None + elif n_mode == 2: + mode, halo = mode + else: + raise ValueError( + "Can't provide more than two positional arguments. " + f"Got: {', '.join(repr(x) for x in mode)}" + ) + + compress = mode is None or mode == "compress" envelope = mode == "envelope" full = mode == "full" + if not (compress or envelope or full): + raise ValueError(f"Invalid mode of operation: {mode!r}") - # Parse halo if halo: try: halo = int(halo) @@ -343,6 +361,7 @@ def _indices(self, mode, halo, data_axes, ancillary_mask, kwargs): ) mask = {} + cyclic = {} for canonical_axes, axes_key_construct_value_id in parsed.items(): axes, keys, constructs, points, identities = tuple( @@ -469,12 +488,15 @@ def _indices(self, mode, halo, data_axes, ancillary_mask, kwargs): f"No indices found from: {identity}={value!r}" ) - index = slice(start, stop, 1) - +# index = slice(start, stop, 1) + index = list(range(start, stop)) + cyclic[axis] = True + if full: - d = self._Data(da.arange(size)) - d.cyclic(0) - ind = (d[index].array,) + # d = self._Data(da.arange(size)) + # d.cyclic(0) +# ind = (d[index].array,) + ind = (index,) index = slice(None) elif item is not None: @@ -651,7 +673,7 @@ def _point_not_in_cell(nodes_x, nodes_y, point): if ind is not None: mask_component_shape = [] masked_subspace_size = 1 - ind = np.array(ind) + ind = np.array(ind, copy=False) for i, (axis, start, stop) in enumerate( zip(canonical_axes, ind.min(axis=1), ind.max(axis=1)) @@ -702,21 +724,13 @@ def _point_not_in_cell(nodes_x, nodes_y, point): index = indices[axis] size = domain_axes[axis].get_size() - try: - index = normalize_cyclic_slice(index, size) - except IndexError: - # Non-cyclic index that is either a slice or a - # list/1-d array of int/bool. - # - # E.g. for halo=1 and size=5: - # slice(1, 3) -> [0, 1, 2, 3] - # slice(1, 4, 2) -> [0, 1, 3, 4] - # slice(2, 0, -1) -> [3, 2, 1, 0] - # [1, 2] -> [0, 1, 2, 3] - # [1, 3] -> [0, 1, 3, 4] - # [2, 1] -> [3, 2, 1, 0] - # [3, 1] -> [4, 3, 1, 0] - # [False, True, False, True, False] -> [0, 1, 3, 4] + # Extend the list at each end + if cyclic.get(axis): + mn = min(index) + mx = max(index) + left = range(max(*(0, mn - halo)), mn) + right = range(mx + 1, min(*(mx + 1 + halo, size))) + else: index = normalize_index((index,), (size,))[0] if isinstance(index, slice): step = index.step @@ -727,13 +741,13 @@ def _point_not_in_cell(nodes_x, nodes_y, point): # Convert dask to numpy, and bool to int. index = np.asanyarray(index) index = normalize_index((index,), (size,))[0] - + increasing = index[0] <= index[-1] # Convert 'index' to a list (which will # exclusively contain non-negative integers) index = index.tolist() - + # Extend the list at each end mn = min(index) mx = max(index) @@ -751,34 +765,97 @@ def _point_not_in_cell(nodes_x, nodes_y, point): right = range( mn - 1, max(*(mn - 1 - halo, -1)), -1 ) - - index[:0] = left - index.extend(right) - else: - # Cyclic slice index - # - # E.g. for halo=1 and size=5: - # slice(-1, 2) -> slice(-2, 3) - # slice(1, -2, -1) -> slice(2, -3, -1) - start = index.start - stop = index.stop - step = index.step - if abs(step) > 1: - raise ValueError( - "A cyclic slice can only have halos if it " - f"has step 1 or -1. Got {index!r}" - ) - - if step > 0: - # Increasing cyclic slice - start = max(start - halo, stop - size) - stop = min(stop + halo, size + start) - else: - # Decreasing cyclic slice - start = min(start + halo, size + stop) - stop = max(stop - halo, start - size) - - index = slice(start, stop, step) + + index[:0] = left + index.extend(right) + + + + +# try: +# index = normalize_cyclic_slice(index, size) +# except IndexError: +# # Non-cyclic index that is either a slice or a +# # list/1-d array of int/bool. +# # +# # E.g. for halo=1 and size=5: +# # slice(1, 3) -> [0, 1, 2, 3] +# # slice(1, 4, 2) -> [0, 1, 3, 4] +# # slice(2, 0, -1) -> [3, 2, 1, 0] +# # [1, 2] -> [0, 1, 2, 3] +# # [1, 3] -> [0, 1, 3, 4] +# # [2, 1] -> [3, 2, 1, 0] +# # [3, 1] -> [4, 3, 1, 0] +# # [False, True, False, True, False] -> [0, 1, 3, 4] +# index = normalize_index((index,), (size,))[0] +# if isinstance(index, slice): +# step = index.step +# increasing = step is None or step > 0 +# index = np.arange(size)[index] +# else: +# if is_dask_collection(index): +# # Convert dask to numpy, and bool to int. +# index = np.asanyarray(index) +# index = normalize_index((index,), (size,))[0] +# +# increasing = index[0] <= index[-1] +# +# # Convert 'index' to a list (which will +# # exclusively contain non-negative integers) +# index = index.tolist() +# +# # Extend the list at each end +# mn = min(index) +# mx = max(index) +# if increasing: +# # "Increasing" sequence: Put smaller +# # indices in the left hand halo, and +# # larger ones in the right hand halo +# left = range(max(*(0, mn - halo)), mn) +# right = range(mx + 1, min(*(mx + 1 + halo, size))) +# else: +# # "Decreasing" sequence: Put larger +# # indices in the left hand halo, and +# # smaller ones in the right hand halo +# left = range(min(*(size - 1, mx + halo)), mx, -1) +# right = range( +# mn - 1, max(*(mn - 1 - halo, -1)), -1 +# ) +# +# index[:0] = left +# index.extend(right) +# else: +# # Cyclic slice index +# # +# # E.g. for halo=1 and size=5: +# # slice(-1, 2) -> slice(-2, 3) +# # slice(1, -2, -1) -> slice(2, -3, -1) +# start = index.start +# stop = index.stop +# step = index.step +# if step not in (1, -1): +# # This restriction is due to the fact that +# # the extended index is a slice (rather +# # than a list of integers), and so we +# # can't represent the uneven spacing that +# # would be required. Note that cyclic +# # slices created by this method will +# # always have a step of 1. +# raise ValueError( +# "A cyclic slice index can only have halos if " +# f"it has step 1 or -1. Got {index!r}" +# ) +# +# if step > 0: +# # Increasing cyclic slice +# start = max(start - halo, stop - size) +# stop = min(stop + halo, size + start) +# else: +# # Decreasing cyclic slice +# start = min(start + halo, size + stop) +# stop = max(stop - halo, start - size) +# +# index = slice(start, stop, step) indices[axis] = index diff --git a/docs/attribute/cf.Domain.subspace.html b/docs/attribute/cf.Domain.subspace.html deleted file mode 100644 index 96308a8d0b..0000000000 --- a/docs/attribute/cf.Domain.subspace.html +++ /dev/null @@ -1,165 +0,0 @@ - - - - - - - - cf.Domain.subspace — Documentation - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - -
-
- - -
- -
-

cf.Domain.subspace

-
-
-Domain.subspace = <function Domain.subspace>[source]
-
- -
- - -
- -
-
-
-
- - - - - - - \ No newline at end of file diff --git a/docs/attribute/cf.Field.indices.html b/docs/attribute/cf.Field.indices.html deleted file mode 100644 index 95e2dce415..0000000000 --- a/docs/attribute/cf.Field.indices.html +++ /dev/null @@ -1,165 +0,0 @@ - - - - - - - - cf.Field.indices — Documentation - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - -
-
- - -
- -
-

cf.Field.indices

-
-
-Field.indices = <function Field.indices>[source]
-
- -
- - -
- -
-
-
-
- - - - - - - \ No newline at end of file diff --git a/docs/method/cf.Field.subspace.html b/docs/method/cf.Field.subspace.html deleted file mode 100644 index ce4b00f785..0000000000 --- a/docs/method/cf.Field.subspace.html +++ /dev/null @@ -1,323 +0,0 @@ - - - - - - - - cf.Field.subspace — Documentation - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - -
-
- - -
- -
-

cf.Field.subspace

-
-
-Field.subspace()
-

Create a subspace of the field construct.

-

Creation of a new field construct which spans a subspace of the -domain of an existing field construct is achieved either by -identifying indices based on the metadata constructs (subspacing -by metadata) or by indexing the field construct directly -(subspacing by index).

-

The subspacing operation, in either case, also subspaces any -metadata constructs of the field construct (e.g. coordinate -metadata constructs) which span any of the domain axis constructs -that are affected. The new field construct is created with the -same properties as the original field construct.

-

Subspacing by metadata

-

Subspacing by metadata, signified by the use of round brackets, -selects metadata constructs and specifies conditions on their -data. Indices for subspacing are then automatically inferred from -where the conditions are met.

-

Metadata constructs and the conditions on their data are defined -by keyword parameters.

-
    -
  • Any domain axes that have not been identified remain unchanged.

  • -
  • Multiple domain axes may be subspaced simultaneously, and it -doesn’t matter which order they are specified in.

  • -
  • Subspace criteria may be provided for size 1 domain axes that -are not spanned by the field construct’s data.

  • -
  • Explicit indices may also be assigned to a domain axis -identified by a metadata construct, with either a Python slice -object, or a sequence of integers or booleans.

  • -
  • For a dimension that is cyclic, a subspace defined by a slice or -by a Query instance is assumed to “wrap” around the edges of -the data.

  • -
  • Conditions may also be applied to multi-dimensional metadata -constructs. The “compress” mode is still the default mode (see -the positional arguments), but because the indices may not be -acting along orthogonal dimensions, some missing data may still -need to be inserted into the field construct’s data.

  • -
-

Subspacing by index

-

Subspacing by indexing, signified by the use of square brackets, -uses rules that are very similar to the numpy indexing rules, the -only differences being:

-
    -
  • An integer index i specified for a dimension reduces the size of -this dimension to unity, taking just the i-th element, but keeps -the dimension itself, so that the rank of the array is not -reduced.

  • -
  • When two or more dimensions’ indices are sequences of integers -then these indices work independently along each dimension -(similar to the way vector subscripts work in Fortran). This is -the same indexing behaviour as on a Variable object of the -netCDF4 package.

  • -
  • For a dimension that is cyclic, a range of indices specified by -a slice that spans the edges of the data (such as -2:3 or -3:-2:-1) is assumed to “wrap” around, rather then producing -a null result.

  • -
-
-

See also

-

indices, squeeze, where, __getitem__

-
-
-
Parameters
-
-
positional arguments: optional

There are three modes of operation, each of which provides -a different type of subspace:

- ---- - - - - - - - - - - - - - - - - - - - -

argument

Description

'compress'

This is the default mode. Unselected -locations are removed to create the -returned subspace. Note that if a -multi-dimensional metadata construct is -being used to define the indices then some -missing data may still be inserted at -unselected locations.

'envelope'

The returned subspace is the smallest that -contains all of the selected -indices. Missing data is inserted at -unselected locations within the envelope.

'full'

The returned subspace has the same domain -as the original field construct. Missing -data is inserted at unselected locations.

'test'

May be used on its own or in addition to -one of the other positional arguments. Do -not create a subspace, but return True -or False depending on whether or not it -is possible to create specified the -subspace.

-
-
keyword parameters: optional

A keyword name is an identity of a metadata construct, and -the keyword value provides a condition for inferring -indices that apply to the dimension (or dimensions) -spanned by the metadata construct’s data. Indices are -created that select every location for which the metadata -construct’s data satisfies the condition.

-
-
-
-
Returns
-
-
Field or bool

An independent field construct containing the subspace of -the original field. If the 'test' positional argument -has been set then return True or False depending on -whether or not it is possible to create specified -subspace.

-
-
-
-
-

Examples

-

There are further worked examples -in the tutorial.

-
>>> g = f.subspace(X=112.5)
->>> g = f.subspace(X=112.5, latitude=cf.gt(-60))
->>> g = f.subspace(latitude=cf.eq(-45) | cf.ge(20))
->>> g = f.subspace(X=[1, 2, 4], Y=slice(None, None, -1))
->>> g = f.subspace(X=cf.wi(-100, 200))
->>> g = f.subspace(X=slice(-2, 4))
->>> g = f.subspace(Y=[True, False, True, True, False])
->>> g = f.subspace(T=410.5)
->>> g = f.subspace(T=cf.dt('1960-04-16'))
->>> g = f.subspace(T=cf.wi(cf.dt('1962-11-01'),
-...                        cf.dt('1967-03-17 07:30')))
->>> g = f.subspace('compress', X=[1, 2, 4, 6])
->>> g = f.subspace('envelope', X=[1, 2, 4, 6])
->>> g = f.subspace('full', X=[1, 2, 4, 6])
->>> g = f.subspace(latitude=cf.wi(51, 53))
-
-
-
>>> g = f.subspace[::-1, 0]
->>> g = f.subspace[:, :, 1]
->>> g = f.subspace[:, 0]
->>> g = f.subspace[..., 6:3:-1, 3:6]
->>> g = f.subspace[0, [2, 3, 9], [4, 8]]
->>> g = t.subspace[0, :, -2]
->>> g = f.subspace[0, [2, 3, 9], [4, 8]]
->>> g = f.subspace[:, -2:3]
->>> g = f.subspace[:, 3:-2:-1]
->>> g = f.subspace[..., [True, False, True, True, False]]
-
-
-
- -
- - -
- -
-
-
-
- - - - - - - \ No newline at end of file diff --git a/docs/source/class/cf.Domain.rst b/docs/source/class/cf.Domain.rst index 4d09b2b047..3117abd118 100644 --- a/docs/source/class/cf.Domain.rst +++ b/docs/source/class/cf.Domain.rst @@ -217,12 +217,6 @@ Subspacing :template: method.rst ~cf.Domain.indices - -.. autosummary:: - :nosignatures: - :toctree: ../attribute/ - :template: attribute.rst - ~cf.Domain.subspace NetCDF From 6bbd81a6d9b69b1bc3d59bafd5f19002821a4521 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Thu, 11 Apr 2024 12:47:51 +0100 Subject: [PATCH 07/30] dev --- cf/mixin/fielddomain.py | 174 ++++++++++++++++++++-------------------- 1 file changed, 87 insertions(+), 87 deletions(-) diff --git a/cf/mixin/fielddomain.py b/cf/mixin/fielddomain.py index 616ee555bc..8af362e366 100644 --- a/cf/mixin/fielddomain.py +++ b/cf/mixin/fielddomain.py @@ -488,15 +488,15 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): f"No indices found from: {identity}={value!r}" ) -# index = slice(start, stop, 1) - index = list(range(start, stop)) + index = slice(start, stop, 1) +# index = list(range(start, stop)) cyclic[axis] = True if full: # d = self._Data(da.arange(size)) # d.cyclic(0) # ind = (d[index].array,) - ind = (index,) + ind = (np.arange(size)[np.arange(start, stop)],) index = slice(None) elif item is not None: @@ -772,90 +772,90 @@ def _point_not_in_cell(nodes_x, nodes_y, point): -# try: -# index = normalize_cyclic_slice(index, size) -# except IndexError: -# # Non-cyclic index that is either a slice or a -# # list/1-d array of int/bool. -# # -# # E.g. for halo=1 and size=5: -# # slice(1, 3) -> [0, 1, 2, 3] -# # slice(1, 4, 2) -> [0, 1, 3, 4] -# # slice(2, 0, -1) -> [3, 2, 1, 0] -# # [1, 2] -> [0, 1, 2, 3] -# # [1, 3] -> [0, 1, 3, 4] -# # [2, 1] -> [3, 2, 1, 0] -# # [3, 1] -> [4, 3, 1, 0] -# # [False, True, False, True, False] -> [0, 1, 3, 4] -# index = normalize_index((index,), (size,))[0] -# if isinstance(index, slice): -# step = index.step -# increasing = step is None or step > 0 -# index = np.arange(size)[index] -# else: -# if is_dask_collection(index): -# # Convert dask to numpy, and bool to int. -# index = np.asanyarray(index) -# index = normalize_index((index,), (size,))[0] -# -# increasing = index[0] <= index[-1] -# -# # Convert 'index' to a list (which will -# # exclusively contain non-negative integers) -# index = index.tolist() -# -# # Extend the list at each end -# mn = min(index) -# mx = max(index) -# if increasing: -# # "Increasing" sequence: Put smaller -# # indices in the left hand halo, and -# # larger ones in the right hand halo -# left = range(max(*(0, mn - halo)), mn) -# right = range(mx + 1, min(*(mx + 1 + halo, size))) -# else: -# # "Decreasing" sequence: Put larger -# # indices in the left hand halo, and -# # smaller ones in the right hand halo -# left = range(min(*(size - 1, mx + halo)), mx, -1) -# right = range( -# mn - 1, max(*(mn - 1 - halo, -1)), -1 -# ) -# -# index[:0] = left -# index.extend(right) -# else: -# # Cyclic slice index -# # -# # E.g. for halo=1 and size=5: -# # slice(-1, 2) -> slice(-2, 3) -# # slice(1, -2, -1) -> slice(2, -3, -1) -# start = index.start -# stop = index.stop -# step = index.step -# if step not in (1, -1): -# # This restriction is due to the fact that -# # the extended index is a slice (rather -# # than a list of integers), and so we -# # can't represent the uneven spacing that -# # would be required. Note that cyclic -# # slices created by this method will -# # always have a step of 1. -# raise ValueError( -# "A cyclic slice index can only have halos if " -# f"it has step 1 or -1. Got {index!r}" -# ) -# -# if step > 0: -# # Increasing cyclic slice -# start = max(start - halo, stop - size) -# stop = min(stop + halo, size + start) -# else: -# # Decreasing cyclic slice -# start = min(start + halo, size + stop) -# stop = max(stop - halo, start - size) -# -# index = slice(start, stop, step) + try: + index = normalize_cyclic_slice(index, size) + except IndexError: + # Non-cyclic index that is either a slice or a + # list/1-d array of int/bool. + # + # E.g. for halo=1 and size=5: + # slice(1, 3) -> [0, 1, 2, 3] + # slice(1, 4, 2) -> [0, 1, 3, 4] + # slice(2, 0, -1) -> [3, 2, 1, 0] + # [1, 2] -> [0, 1, 2, 3] + # [1, 3] -> [0, 1, 3, 4] + # [2, 1] -> [3, 2, 1, 0] + # [3, 1] -> [4, 3, 1, 0] + # [False, True, False, True, False] -> [0, 1, 3, 4] + index = normalize_index((index,), (size,))[0] + if isinstance(index, slice): + step = index.step + increasing = step is None or step > 0 + index = np.arange(size)[index] + else: + if is_dask_collection(index): + # Convert dask to numpy, and bool to int. + index = np.asanyarray(index) + index = normalize_index((index,), (size,))[0] + + increasing = index[0] <= index[-1] + + # Convert 'index' to a list (which will + # exclusively contain non-negative integers) + index = index.tolist() + + # Extend the list at each end + mn = min(index) + mx = max(index) + if increasing: + # "Increasing" sequence: Put smaller + # indices in the left hand halo, and + # larger ones in the right hand halo + left = range(max(*(0, mn - halo)), mn) + right = range(mx + 1, min(*(mx + 1 + halo, size))) + else: + # "Decreasing" sequence: Put larger + # indices in the left hand halo, and + # smaller ones in the right hand halo + left = range(min(*(size - 1, mx + halo)), mx, -1) + right = range( + mn - 1, max(*(mn - 1 - halo, -1)), -1 + ) + + index[:0] = left + index.extend(right) + else: + # Cyclic slice indexself.iscyclic(axis) + # + # E.g. for halo=1 and size=5: + # slice(-1, 2) -> slice(-2, 3) + # slice(1, -2, -1) -> slice(2, -3, -1) + start = index.start + stop = index.stop + step = index.step + if step not in (1, -1): + # This restriction is due to the fact that + # the extended index is a slice (rather + # than a list of integers), and so we + # can't represent the uneven spacing that + # would be required. Note that cyclic + # slices created by this method will + # always have a step of 1. + raise ValueError( + "A cyclic slice index can only have halos if " + f"it has step 1 or -1. Got {index!r}" + ) + + if step > 0: + # Increasing cyclic slice + start = max(start - halo, stop - size) + stop = min(stop + halo, size + start) + else: + # Decreasing cyclic slice + start = min(start + halo, size + stop) + stop = max(stop - halo, start - size) + + index = slice(start, stop, step) indices[axis] = index From 28ad356a9b5f8d433745065f9ba4db5f02ac26b0 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Thu, 11 Apr 2024 16:55:25 +0100 Subject: [PATCH 08/30] dev --- cf/docstring/docstring.py | 54 ++++++----- cf/field.py | 139 +++++++++++++------------- cf/mixin/fielddomain.py | 199 ++++++++++++++++---------------------- cf/test/test_Field.py | 11 +++ 4 files changed, 192 insertions(+), 211 deletions(-) diff --git a/cf/docstring/docstring.py b/cf/docstring/docstring.py index fd03a2ce8c..af2c05f800 100644 --- a/cf/docstring/docstring.py +++ b/cf/docstring/docstring.py @@ -92,7 +92,7 @@ The number of extra elements will be automatically reduced if including full amount defined by the halo would extend the - subspace beyond the axis's extent.""", + subspace beyond the axis limits.""", # indices halos "{{subspace halos}}": """ If a halo is defined via a positional argument, then each @@ -104,7 +104,7 @@ The number of extra elements will be automatically reduced if including full amount defined by the halo would extend the - subspace beyond the axis's extent.""", + subspace beyond the axis limits.""", # ---------------------------------------------------------------- # Method description substitutions (3 levels of indentation) # ---------------------------------------------------------------- @@ -657,23 +657,19 @@ provided then assume the ``'compress'`` mode of operation with no halo added to the subspaced - axes. Note that this is equivalent to - either of ``'compress', 0`` and ``0``. + axes. ``mode`` Define the mode of operation with no - halo added to the subspaced axes. Note - that ``mode`` is equivalent to ``mode, - 0``. + halo added to the subspaced axes. ``mode, halo`` Define a mode of operation, as well as a halo to be added to the subspaced - axes. Note that ``mode, 0`` is - equivalent to ``mode``. + axes. ``halo`` Assume the ``'compress'`` mode of operation and define a halo to be - added to the subspaced axes. Note - that ``halo`` is equivalent to + added to the subspaced axes. Note that + ``halo`` is equivalent to ``'compress', halo``. ============== ======================================""", # ---------------------------------------------------------------- @@ -731,34 +727,40 @@ Valid modes are: * ``'compress'`` This the default. + Unselected locations are removed to create the - subspace. Missing data may be inserted at - unselected locations, but only where it is - required to ensure a hyperrectangular result, and - not if a non-zero halo has been defined. + subspace. If the result is not hperrecctangular + then the minimum amount of unselected locations + required to make it so will also be specially + selected. Missing data is inserted at the + specially selected locations, unless a halo has + been defined (of any size, including 0). * ``'envelope'`` - The subspace is the smallest subspace that - contains all of the selected locations. Missing - data is inserted at unselected locations within - the envelope, unless a non-zero halo has been - defined. + The subspace is the smallest hperrecctangular + subspace that contains all of the selected + locations. Missing data is inserted at unselected + locations within the envelope, unless a halo has + been defined (of any size, including 0). * ``'full'`` The subspace has the same domain as the original construct. Missing data is inserted at unselected - locations, unless a non-zero halo has been - defined.""", + locations, unless a halo has been defined (of any + size, including 0).""", # indices valid modes Domain "{{indices valid modes Domain}}": """ Valid modes are: * ``'compress'`` This the default. Unselected locations are removed to create the - subspace, unless they are required to ensure a - hyperrectangular result. + subspace. If the result is not hperrecctangular + then the minimum amount of unselected locations + required to make it so will also be specially + selected. * ``'envelope'`` - The subspace is the smallest subspace that - contains all of the selected locations.""", + The subspace is the smallest hperrecctangular + subspace that contains all of the selected + locations.""", } diff --git a/cf/field.py b/cf/field.py index eab4169ef8..189832aeca 100644 --- a/cf/field.py +++ b/cf/field.py @@ -8864,8 +8864,8 @@ def indices(self, *mode, **kwargs): {{indices halos}} - If a non-zero halo has been defined then no ancillary masks - will be created. + If a halo has been defined (of any size, including 0), then no + ancillary masks will be created. .. versionadded:: 1.0 @@ -9066,92 +9066,92 @@ def set_data( ): """Set the field construct data. - .. versionadded:: 3.0.0 + .. versionadded:: 3.0.0 - .. seealso:: `data`, `del_data`, `get_data`, `has_data`, - `set_construct` + .. seealso:: `data`, `del_data`, `get_data`, `has_data`, + `set_construct` - :Parameters: + :Parameters: - data: `Data` - The data to be inserted. + data: `Data` + The data to be inserted. - {{data_like}} + {{data_like}} - axes: (sequence of) `str` or `int`, optional - Set the domain axes constructs that are spanned by the - data. If unset, and the *set_axes* parameter is True, then - an attempt will be made to assign existing domain axis - constructs to the data. + axes: (sequence of) `str` or `int`, optional + Set the domain axes constructs that are spanned by the + data. If unset, and the *set_axes* parameter is True, then + an attempt will be made to assign existing domain axis + constructs to the data. - The contents of the *axes* parameter is mapped to domain - axis constructs by translating each element into a domain - axis construct key via the `domain_axis` method. + The contents of the *axes* parameter is mapped to domain + axis constructs by translating each element into a domain + axis construct key via the `domain_axis` method. - *Parameter example:* - ``axes='domainaxis1'`` + *Parameter example:* + ``axes='domainaxis1'`` - *Parameter example:* - ``axes='X'`` + *Parameter example:* + ``axes='X'`` - *Parameter example:* - ``axes=['latitude']`` + *Parameter example:* + ``axes=['latitude']`` - *Parameter example:* - ``axes=['X', 'longitude']`` + *Parameter example:* + ``axes=['X', 'longitude']`` - *Parameter example:* - ``axes=[1, 0]`` + *Parameter example:* + ``axes=[1, 0]`` - set_axes: `bool`, optional - If False then do not set the domain axes constructs that - are spanned by the data, even if the *axes* parameter has - been set. By default the axes are set either according to - the *axes* parameter, or if any domain axis constructs - exist then an attempt will be made to assign existing - domain axis constructs to the data. + set_axes: `bool`, optional + If False then do not set the domain axes constructs that + are spanned by the data, even if the *axes* parameter has + been set. By default the axes are set either according to + the *axes* parameter, or if any domain axis constructs + exist then an attempt will be made to assign existing + domain axis constructs to the data. - If the *axes* parameter is `None` and no domain axis - constructs exist then no attempt is made to assign domain - axes constructs to the data, regardless of the value of - s *set_axes*. + If the *axes* parameter is `None` and no domain axis + constructs exist then no attempt is made to assign domain + axes constructs to the data, regardless of the value of + *set_axes*. - copy: `bool`, optional - If True then set a copy of the data. By default the data - are copied. + copy: `bool`, optional + If True then set a copy of the data. By default the data + are copied. - {{inplace: `bool`, optional (default True)}} + {{inplace: `bool`, optional (default True)}} - .. versionadded:: 3.7.0 + .. versionadded:: 3.7.0 - :Returns: + :Returns: - `None` or `Field` - If the operation was in-place then `None` is returned, - otherwise return a new `Field` instance containing the new - data. + `None` or `Field` + If the operation was in-place then `None` is returned, + otherwise return a new `Field` instance containing the new + data. - **Examples** + **Examples** - >>> f = cf.Field() - >>> f.set_data([1, 2, 3]) - >>> f.has_data() - True - >>> f.get_data() - - >>> f.data - - >>> f.del_data() - - >>> g = f.set_data([4, 5, 6], inplace=False) - >>> g.data - - >>> f.has_data() - False - >>> print(f.get_data(None)) - None - >>> print(f.del_data(None)) - None + >>> f = cf.Field() + >>> f.set_data([1, 2, 3]) + >>> f.has_data() + True + >>> f.get_data() + + >>> f.data + + >>> f.del_data() + + >>> g = f.set_data([4, 5, 6], inplace=False) + >>> g.data + + >>> f.has_data() + False + >>> print(f.get_data(None)) + None + >>> print(f.del_data(None)) + None """ data = self._Data(data, copy=False) @@ -13268,9 +13268,6 @@ def subspace(self): {{subspace halos}} - .. note:: If a non-zero halo has been defined then no - ancillary masks will be created. - .. versionadded:: 1.0 .. seealso:: `indices`, `where`, `__getitem__`, diff --git a/cf/mixin/fielddomain.py b/cf/mixin/fielddomain.py index 8af362e366..c35bfe992a 100644 --- a/cf/mixin/fielddomain.py +++ b/cf/mixin/fielddomain.py @@ -1,7 +1,6 @@ import logging from numbers import Integral -import dask.array as da import numpy as np from cfdm import is_log_level_debug, is_log_level_info from dask.array.slicing import normalize_index @@ -259,13 +258,13 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): n_mode = len(mode) if not n_mode: mode = None - halo = 0 + halo = None elif n_mode == 1: try: halo = int(mode[0]) except ValueError: mode = mode[0] - halo = 0 + halo = None else: mode = None elif n_mode == 2: @@ -282,7 +281,7 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): if not (compress or envelope or full): raise ValueError(f"Invalid mode of operation: {mode!r}") - if halo: + if halo is not None: try: halo = int(halo) except ValueError: @@ -361,7 +360,6 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): ) mask = {} - cyclic = {} for canonical_axes, axes_key_construct_value_id in parsed.items(): axes, keys, constructs, points, identities = tuple( @@ -417,18 +415,18 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): # [7,4,2], slice(0,4,2), # numpy.array([2,4,7]), # [True,False,True] - index = value if debug: logger.debug(" 1-d CASE 1:") # pragma: no cover index = value if envelope or full: + # Set ind size = domain_axes[axis].get_size() - # TODODASK: consider using dask.arange here - d = np.arange(size) # self._Data(range(size)) - ind = (d[value],) # .array,) - index = slice(None) + d = np.arange(size) + ind = (d[value],) + # Placeholder which will be overwritten later + index = None elif ( item is not None @@ -489,15 +487,24 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): ) index = slice(start, stop, 1) -# index = list(range(start, stop)) - cyclic[axis] = True - + if full: - # d = self._Data(da.arange(size)) - # d.cyclic(0) -# ind = (d[index].array,) - ind = (np.arange(size)[np.arange(start, stop)],) - index = slice(None) + # Set ind + try: + index = normalize_cyclic_slice(index, size) + except IndexError: + # Index is not a cyclic slice + ind = (np.arange(size)[index],) + else: + # Index is a cyclic slice + ind = ( + np.arange(size)[ + np.arange(index.start, index.stop) + ], + ) + + # Placeholder which will be overwritten later + index = None elif item is not None: # 1-d CASE 3: All other 1-d cases @@ -505,16 +512,18 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): logger.debug(" 1-d CASE 3:") # pragma: no cover index = item == value - index = index.data.to_dask_array() + index = index.to_dask_array() if envelope or full: + # Set ind index = np.asanyarray(index) if np.ma.isMA(index): ind = np.ma.where(index) else: ind = np.where(index) - index = slice(None) + # Placeholder which will be overwritten later + index = None else: raise ValueError( @@ -570,17 +579,22 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): ] # Find loctions that are True in all of the - # construct's matches + # constructs' matches item_match = item_matches.pop() for m in item_matches: item_match &= m - item_match = item_match.compute() + # Set ind + item_match = np.asanyarray(item_match) if np.ma.isMA(item_match): ind = np.ma.where(item_match) else: ind = np.where(item_match) + # Placeholders which will be overwritten later + for axis in canonical_axes: + indices[axis] = None + if debug: logger.debug( f" item_match = {item_match}\n" @@ -600,13 +614,13 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): if item.has_bounds() ] - # If there are exactly two 2-d contructs constructs, - # both with cell bounds and both with 'cf.contains' - # values, then do an extra check to remove any cells - # already selected for which the given value is in - # fact outside of the cell. This could happen if the - # cells are not rectangular (e.g. for curvilinear - # latitudes and longitudes arrays). + # If there are exactly two 2-d constructs, both with + # cell bounds and both with 'cf.contains' values, then + # do an extra check to remove any cells already + # selected for which the given value is in fact + # outside of the cell. This could happen if the cells + # are not rectilinear (e.g. for curvilinear latitudes + # and longitudes arrays). if n_items == constructs[0].ndim == len(bounds) == 2: point2 = [] for v, construct in zip(points, transposed_constructs): @@ -681,28 +695,35 @@ def _point_not_in_cell(nodes_x, nodes_y, point): if data_axes and axis not in data_axes: continue - if indices[axis] == slice(None): - if compress: - # Create a compressed index for this axis - size = stop - start + 1 - index = sorted(set(ind[i])) - elif envelope: - # Create an envelope index for this axis - stop += 1 - size = stop - start - index = slice(start, stop) - elif full: - # Create a full index for this axis - start = 0 - stop = domain_axes[axis].get_size() - size = stop - start - index = slice(None) - else: - raise ValueError( - "Must have mode full, envelope, or compress" - ) # pragma: no cover + if compress: + # Create a compressed index for this axis + size = stop - start + 1 + index = sorted(set(ind[i])) + elif envelope: + # Create an envelope index for this axis + stop += 1 + size = stop - start + index = slice(start, stop) + elif full: + # Create a full index for this axis + start = 0 + stop = domain_axes[axis].get_size() + size = stop - start + index = slice(None) + else: + raise ValueError( + "Must have mode full, envelope, or compress" + ) # pragma: no cover - indices[axis] = index + # Overwrite the placeholder value of None + if indices[axis] is not None: + raise ValueError( + "This error means that there is a bug: The " + "'indices' dictionary should contain None for " + "each axis with an 'ind'." + ) + + indices[axis] = index mask_component_shape.append(size) masked_subspace_size *= size @@ -710,6 +731,7 @@ def _point_not_in_cell(nodes_x, nodes_y, point): create_mask = ( ancillary_mask + and halo is None and data_axes and ind.shape[1] < masked_subspace_size ) @@ -717,66 +739,21 @@ def _point_not_in_cell(nodes_x, nodes_y, point): create_mask = False # -------------------------------------------------------- - # TODO + # Add a halo to the subspaced axes # -------------------------------------------------------- if halo: for axis in item_axes: index = indices[axis] size = domain_axes[axis].get_size() - # Extend the list at each end - if cyclic.get(axis): - mn = min(index) - mx = max(index) - left = range(max(*(0, mn - halo)), mn) - right = range(mx + 1, min(*(mx + 1 + halo, size))) - else: - index = normalize_index((index,), (size,))[0] - if isinstance(index, slice): - step = index.step - increasing = step is None or step > 0 - index = np.arange(size)[index] - else: - if is_dask_collection(index): - # Convert dask to numpy, and bool to int. - index = np.asanyarray(index) - index = normalize_index((index,), (size,))[0] - - increasing = index[0] <= index[-1] - - # Convert 'index' to a list (which will - # exclusively contain non-negative integers) - index = index.tolist() - - # Extend the list at each end - mn = min(index) - mx = max(index) - if increasing: - # "Increasing" sequence: Put smaller - # indices in the left hand halo, and - # larger ones in the right hand halo - left = range(max(*(0, mn - halo)), mn) - right = range(mx + 1, min(*(mx + 1 + halo, size))) - else: - # "Decreasing" sequence: Put larger - # indices in the left hand halo, and - # smaller ones in the right hand halo - left = range(min(*(size - 1, mx + halo)), mx, -1) - right = range( - mn - 1, max(*(mn - 1 - halo, -1)), -1 - ) - - index[:0] = left - index.extend(right) - - - - try: + # Is index a cyclic slice? index = normalize_cyclic_slice(index, size) except IndexError: - # Non-cyclic index that is either a slice or a - # list/1-d array of int/bool. + # Non-cyclic index + # + # Either a slice or a list/1-d array of + # int/bool. # # E.g. for halo=1 and size=5: # slice(1, 3) -> [0, 1, 2, 3] @@ -825,7 +802,7 @@ def _point_not_in_cell(nodes_x, nodes_y, point): index[:0] = left index.extend(right) else: - # Cyclic slice indexself.iscyclic(axis) + # Cyclic slice index # # E.g. for halo=1 and size=5: # slice(-1, 2) -> slice(-2, 3) @@ -838,9 +815,9 @@ def _point_not_in_cell(nodes_x, nodes_y, point): # the extended index is a slice (rather # than a list of integers), and so we # can't represent the uneven spacing that - # would be required. Note that cyclic - # slices created by this method will - # always have a step of 1. + # would be required if abs(step) != 1. + # Note that cyclic slices created by this + # method will always have a step of 1. raise ValueError( "A cyclic slice index can only have halos if " f"it has step 1 or -1. Got {index!r}" @@ -866,15 +843,9 @@ def _point_not_in_cell(nodes_x, nodes_y, point): ) # pragma: no cover if create_mask: - if halo: - logger.warn( - "Not applying ancillary masks to a subspace with " - "non-zero halos" - ) - else: - mask[canonical_axes] = _create_ancillary_mask_component( - mask_component_shape, ind, compress - ) + mask[canonical_axes] = _create_ancillary_mask_component( + mask_component_shape, ind, compress + ) indices = {"indices": indices, "mask": mask} diff --git a/cf/test/test_Field.py b/cf/test/test_Field.py index 8cf38c0966..06eb40b594 100644 --- a/cf/test/test_Field.py +++ b/cf/test/test_Field.py @@ -1646,6 +1646,17 @@ def test_Field_indices(self): (x == [120, 80, 40, 0, -40, -80, -120, -160, -200]).all() ) + # Halos: ancillary masking + index = cf.wi(90, 100) + indices = f.indices(longitude=index) + g = f[indices] + self.assertTrue(np.ma.is_masked(g.array)) + + for halo in (0, 1): + indices = f.indices(halo, longitude=index) + g = f[indices] + self.assertFalse(np.ma.is_masked(g.array)) + # Subspace has size 0 axis resulting from dask array index indices = f.indices(grid_latitude=cf.contains(-23.2)) with self.assertRaises(IndexError): From 62e8e622211de23153d88ca310115565d8d72c44 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Fri, 12 Apr 2024 09:51:31 +0100 Subject: [PATCH 09/30] dev --- cf/docstring/docstring.py | 48 ++++++++++++--------------------------- cf/domain.py | 11 +++++---- cf/field.py | 25 ++++++++++++++------ cf/functions.py | 3 ++- cf/mixin/fielddomain.py | 8 ++++++- 5 files changed, 47 insertions(+), 48 deletions(-) diff --git a/cf/docstring/docstring.py b/cf/docstring/docstring.py index af2c05f800..bae034f5d1 100644 --- a/cf/docstring/docstring.py +++ b/cf/docstring/docstring.py @@ -80,31 +80,14 @@ Whether `esmpy` logging is enabled or not is determined by `cf.regrid_logging`. If it is enabled then logging takes place after every call. By default logging is disabled.""", - # indices halos - "{{indices halos}}": """ - If a halo is defined via a positional argument, then the - returned index for each subspaced axis will be extended to - include that many extra elements at each "side" of the - axis. For instance, ``f.indices(X=slice(10, 20))`` will give - identical results to each of ``f.indices(0, X=slice(10, - 20))``, ``f.indices(1, X=slice(11, 19))``, ``f.indices(2, - X=slice(12, 18))``, etc. - - The number of extra elements will be automatically reduced if - including full amount defined by the halo would extend the - subspace beyond the axis limits.""", - # indices halos + # subspace halos "{{subspace halos}}": """ If a halo is defined via a positional argument, then each subspaced axis will be extended to include that many extra - elements at each "side" of the axis. For instance, - ``f.subspace(X=slice(10, 20))`` will give identical results to - each of ``f.subspace(0, X=slice(10, 20))``, ``f.subspace(1, - X=slice(11, 19))``, ``f.subspace(2, X=slice(12, 18))``, etc. - - The number of extra elements will be automatically reduced if - including full amount defined by the halo would extend the - subspace beyond the axis limits.""", + elements at each "side" of the axis. The number of extra + elements will be automatically reduced if including the full + amount defined by the halo would extend the subspace beyond + the axis limits.""", # ---------------------------------------------------------------- # Method description substitutions (3 levels of indentation) # ---------------------------------------------------------------- @@ -638,7 +621,7 @@ "{{to_size: `int`, optional}}": """to_size: `int`, optional Pad the axis after so that the new axis has the given size.""", - # indices mode options + # subspace mode options "{{mode: optional}}": """mode: optional Specify the mode of operation (``mode``) and a halo to be added to the subspaced axes (``halo``) with @@ -722,14 +705,12 @@ The removed CFA-netCDF file name substitution. If the substitution was not defined then an empty dictionary is returned.""", - # indices valid modes Field - "{{indices valid modes Field}}": """ - Valid modes are: + # subspace valid modes Field + "{{subspace valid modes Field}}": """Valid modes are: * ``'compress'`` This the default. - Unselected locations are removed to create the - subspace. If the result is not hperrecctangular + subspace. If the result is not hyperrectangular then the minimum amount of unselected locations required to make it so will also be specially selected. Missing data is inserted at the @@ -737,7 +718,7 @@ been defined (of any size, including 0). * ``'envelope'`` - The subspace is the smallest hperrecctangular + The subspace is the smallest hyperrectangular subspace that contains all of the selected locations. Missing data is inserted at unselected locations within the envelope, unless a halo has @@ -748,19 +729,18 @@ construct. Missing data is inserted at unselected locations, unless a halo has been defined (of any size, including 0).""", - # indices valid modes Domain - "{{indices valid modes Domain}}": """ - Valid modes are: + # subspace valid modes Domain + "{{subspace valid modes Domain}}": """Valid modes are: * ``'compress'`` This the default. Unselected locations are removed to create the - subspace. If the result is not hperrecctangular + subspace. If the result is not hyperrectangular then the minimum amount of unselected locations required to make it so will also be specially selected. * ``'envelope'`` - The subspace is the smallest hperrecctangular + The subspace is the smallest hyperrectangular subspace that contains all of the selected locations.""", } diff --git a/cf/domain.py b/cf/domain.py index d1243fcd50..db9e1e1199 100644 --- a/cf/domain.py +++ b/cf/domain.py @@ -780,7 +780,7 @@ def indices(self, *mode, **kwargs): **Halos** - {{indices halos}} + {{subspace halos}} .. versionadded:: 3.11.0 @@ -790,9 +790,10 @@ def indices(self, *mode, **kwargs): :Parameters: {{mode: optional}} - {{indices valid modes Domain}} - kwargs: *optional* + {{subspace valid modes Domain}} + + kwargs: optional A keyword name is an identity of a metadata construct, and the keyword value provides a condition for inferring indices that apply to the dimension (or @@ -1145,9 +1146,9 @@ def subspace(self, *mode, **kwargs): {{mode: optional}} - {{indices valid modes Domain}} + {{subspace valid modes Domain}} - kwargs: *optional* + kwargs: optional A keyword name is an identity of a metadata construct, and the keyword value provides a condition for inferring indices that apply to the dimension (or dimensions) diff --git a/cf/field.py b/cf/field.py index 189832aeca..c4a583053f 100644 --- a/cf/field.py +++ b/cf/field.py @@ -8862,7 +8862,12 @@ def indices(self, *mode, **kwargs): **Halos** - {{indices halos}} + {{subspace halos}} + + For instance, ``f.indices(X=slice(10, 20))`` will give + identical results to each of ``f.indices(0, X=slice(10, + 20))``, ``f.indices(1, X=slice(11, 19))``, ``f.indices(2, + X=slice(12, 18))``, etc. If a halo has been defined (of any size, including 0), then no ancillary masks will be created. @@ -8876,9 +8881,9 @@ def indices(self, *mode, **kwargs): {{mode: optional}} - {{indices valid modes Field}} + {{subspace valid modes Field}} - kwargs: *optional* + kwargs: optional A keyword name is an identity of a metadata construct, and the keyword value provides a condition for inferring indices that apply to the dimension (or @@ -9011,9 +9016,10 @@ def indices(self, *mode, **kwargs): removed_at="4.0.0", ) # pragma: no cover + data_axes = self.get_data_axes() + # Get the indices for every domain axis in the domain, # including any ancillary masks - data_axes = self.get_data_axes() domain_indices = self._indices(mode, data_axes, True, kwargs) # Initialise the output indices with any ancillary masks. @@ -13268,6 +13274,11 @@ def subspace(self): {{subspace halos}} + For instance, ``f.subspace(X=slice(10, 20))`` will give + identical results to each of ``f.subspace(0, X=slice(10, + 20))``, ``f.subspace(1, X=slice(11, 19))``, ``f.subspace(2, + X=slice(12, 18))``, etc. + .. versionadded:: 1.0 .. seealso:: `indices`, `where`, `__getitem__`, @@ -13277,15 +13288,15 @@ def subspace(self): {{mode: optional}} - {{indices valid modes Field}} + {{subspace valid modes Field}} In addition, an extra postional argument of ``'test'`` - is allowed. When provided the subspace is not + is allowed. When provided, the subspace is not returned, instead `True` or `False` is returned depending on whether or not it is possible for the requested subspace to be created. - keyword parameters: *optional* + keyword parameters: optional A keyword name is an identity of a metadata construct, and the keyword value provides a condition for inferring indices that apply to the dimension (or dimensions) diff --git a/cf/functions.py b/cf/functions.py index c06c402b56..002777edbf 100644 --- a/cf/functions.py +++ b/cf/functions.py @@ -1969,9 +1969,10 @@ def parse_indices(shape, indices, cyclic=False, keepdims=True): try: index = normalize_cyclic_slice(index, size) except IndexError: + # Non-cyclic slice pass else: - # Got a cyclic slice + # Cyclic slice start = index.start stop = index.stop step = index.step diff --git a/cf/mixin/fielddomain.py b/cf/mixin/fielddomain.py index c35bfe992a..56ab3821ed 100644 --- a/cf/mixin/fielddomain.py +++ b/cf/mixin/fielddomain.py @@ -499,7 +499,9 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): # Index is a cyclic slice ind = ( np.arange(size)[ - np.arange(index.start, index.stop) + np.arange( + index.start, index.stop, index.step + ) ], ) @@ -742,6 +744,10 @@ def _point_not_in_cell(nodes_x, nodes_y, point): # Add a halo to the subspaced axes # -------------------------------------------------------- if halo: + # Note: We're about to make 'indices' inconsistent + # with 'ind', but that's OK because we're not + # going to use 'ind' again since 'create_mask' + # is False. for axis in item_axes: index = indices[axis] size = domain_axes[axis].get_size() From 3c90399fd7daca096386bd9da840218316b53d87 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Fri, 12 Apr 2024 11:00:27 +0100 Subject: [PATCH 10/30] dev --- cf/mixin/fielddomain.py | 1 + 1 file changed, 1 insertion(+) diff --git a/cf/mixin/fielddomain.py b/cf/mixin/fielddomain.py index 56ab3821ed..12d9e9ce3f 100644 --- a/cf/mixin/fielddomain.py +++ b/cf/mixin/fielddomain.py @@ -689,6 +689,7 @@ def _point_not_in_cell(nodes_x, nodes_y, point): if ind is not None: mask_component_shape = [] masked_subspace_size = 1 + # TODONUMPY2: https://numpy.org/devdocs/numpy_2_0_migration_guide.html#adapting-to-changes-in-the-copy-keyword ind = np.array(ind, copy=False) for i, (axis, start, stop) in enumerate( From 0b0c3c3c5a65a565bc733544eda801aefbce1870 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Fri, 12 Apr 2024 11:46:01 +0100 Subject: [PATCH 11/30] dev --- Changelog.rst | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/Changelog.rst b/Changelog.rst index 3a0b2560c6..c76f1f3520 100644 --- a/Changelog.rst +++ b/Changelog.rst @@ -3,6 +3,15 @@ version NEXT **2024-??-??** +* Allow a halo to be added by `cf.Field.indices` and + `cf.Field.subspace` + (https://github.com/NCAS-CMS/cf-python/issues/759) + +version NEXT +------------ + +**2024-??-??** + * Added spherical regridding to discrete sampling geometry destination grids (https://github.com/NCAS-CMS/cf-python/issues/716) * Added 3-d spherical regridding to `cf.Field.regrids`, and the option From bde2e2843d7045dd8aa007749b0199bfce4dfd7e Mon Sep 17 00:00:00 2001 From: David Hassell Date: Mon, 15 Apr 2024 13:18:10 +0100 Subject: [PATCH 12/30] dev --- cf/docstring/docstring.py | 7 +- cf/domain.py | 3 - cf/functions.py | 55 ++++++++---- cf/mixin/fielddomain.py | 175 +++++++++++++++++++++++++------------- cf/test/test_Field.py | 70 +++++++++------ cf/test/test_functions.py | 30 +++++-- 6 files changed, 223 insertions(+), 117 deletions(-) diff --git a/cf/docstring/docstring.py b/cf/docstring/docstring.py index bae034f5d1..5cc0dd0d8e 100644 --- a/cf/docstring/docstring.py +++ b/cf/docstring/docstring.py @@ -639,8 +639,7 @@ ```` If no positional arguments are provided then assume the ``'compress'`` mode of operation with - no halo added to the subspaced - axes. + no halo added to the subspaced axes. ``mode`` Define the mode of operation with no halo added to the subspaced axes. @@ -651,9 +650,7 @@ ``halo`` Assume the ``'compress'`` mode of operation and define a halo to be - added to the subspaced axes. Note that - ``halo`` is equivalent to - ``'compress', halo``. + added to the subspaced axes. ============== ======================================""", # ---------------------------------------------------------------- # Method description substitutions (4 levels of indentation) diff --git a/cf/domain.py b/cf/domain.py index db9e1e1199..2824d07ebf 100644 --- a/cf/domain.py +++ b/cf/domain.py @@ -1117,9 +1117,6 @@ def subspace(self, *mode, **kwargs): * Multiple domain axes may be subspaced simultaneously, and it doesn't matter which order they are specified in. - * Subspace criteria may be provided for size 1 domain axes that - are not spanned by the field construct's data. - * Explicit indices may also be assigned to a domain axis identified by a metadata construct, with either a Python `slice` object, or a sequence of integers or booleans. diff --git a/cf/functions.py b/cf/functions.py index 002777edbf..bdcce6dbc6 100644 --- a/cf/functions.py +++ b/cf/functions.py @@ -1967,7 +1967,7 @@ def parse_indices(shape, indices, cyclic=False, keepdims=True): if cyclic and isinstance(index, slice): # Check for a cyclic slice try: - index = normalize_cyclic_slice(index, size) + index = normalize_slice(index, size, cyclic=True) except IndexError: # Non-cyclic slice pass @@ -2030,12 +2030,11 @@ def parse_indices(shape, indices, cyclic=False, keepdims=True): return parsed_indices, roll -def normalize_cyclic_slice(index, size): - """Normalise a cyclic slice. +def normalize_slice(index, size, cyclic=False): + """Normalise a slice. - If the *index* is found to be a cyclic slice then a normalised - copy if it is returned. If a non-cyclic slice *index* is found - then an `IndexError` is raised. + If *index* is not a slice, or *cyclic* is True and *index* is not + a cylcic slice, then an `IndexError` is raised. .. versionadded:: NEXTRELEASE @@ -2047,32 +2046,48 @@ def normalize_cyclic_slice(index, size): size: `int` The size of the axis to which the *index* applies. + cyclic: `bool` + If True then normalise a cyclic slice, and raise an + exception if the *index* is a not a cyclic slice. + :Returns: `slice` - The normalised cyclic slice. + The normalised slice. **Examples** - >>> cf.normalize_cyclic_slice(slice(-2, 3), 8) + >>> cf.normalize_slice(slice(1, 4), 8) + slice(1, 4, 1) + >>> cf.normalize_slice(slice(None), 8) + slice(0, 8, 1) + >>> cf.normalize_slice(slice(6, None, -1), 8) + slice(6, None, -1) + >>> cf.normalize_slice(slice(-2, 4), 8) + slice(6, 4, 1) + + >>> cf.normalize_slice([1, 2], 8) + IndexError: [1, 2] is not a slice + + >>> cf.normalize_slice(slice(-2, 3), 8, cyclic=True) slice(-2, 3, 1) - >>> cf.normalize_cyclic_slice(slice(6, 3), 8) + >>> cf.normalize_slice(slice(6, 3), 8, cyclic=True) slice(-2, 3, 1) - >>> cf.normalize_cyclic_slice(slice(2, -3, -1), 8) + >>> cf.normalize_slice(slice(2, -3, -1), 8, cyclic=True) slice(2, -3, -1) - >>> cf.normalize_cyclic_slice(slice(2, 5, -1), 8) + >>> cf.normalize_slice(slice(2, 5, -1), 8, cyclic=True) slice(2, -3, -1) - >>> cf.normalize_cyclic_slice(slice(1, 6), 8) + >>> cf.normalize_slice(slice(1, 6), 8, cyclic=True) IndexError: slice(1, 6, None) is not a cyclic slice - >>> cf.normalize_cyclic_slice([1, 2, 3, 4], 8) + >>> cf.normalize_slice([1, 2, 3, 4], 8, cyclic=True) IndexError: [1, 2, 3, 4] is not a cyclic slice """ if not isinstance(index, slice): step = 0 - else: + elif cyclic: start = index.start stop = index.stop step = index.step @@ -2108,13 +2123,23 @@ def normalize_cyclic_slice(index, size): # 3:6:-1 => 3:-4:-1 # 3:9:-1 => 3:-1:-1 stop -= size + else: + start, stop, step = index.indices(size) + if step < 0 and stop < 0: + stop = None + + # Return the normalized non-cyclic slice + return slice(start, stop, step) if not ( (step > 0 and start < 0 and stop > 0) or (step < 0 and start > 0 and stop < 0) ): - raise IndexError(f"{index!r} is not a cyclic slice") + raise IndexError( + f"{index!r} is not a {'cyclic ' if cyclic else ''}slice" + ) + # Return the normalized cyclic slice return slice(start, stop, step) diff --git a/cf/mixin/fielddomain.py b/cf/mixin/fielddomain.py index 12d9e9ce3f..2046f81669 100644 --- a/cf/mixin/fielddomain.py +++ b/cf/mixin/fielddomain.py @@ -16,7 +16,7 @@ from ..functions import ( _DEPRECATION_ERROR_KWARGS, bounds_combination_mode, - normalize_cyclic_slice, + normalize_slice, ) from ..query import Query from ..units import Units @@ -254,7 +254,7 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): """ debug = is_log_level_debug(logger) - # Parse mode + # Parse mode and halo n_mode = len(mode) if not n_mode: mode = None @@ -423,8 +423,7 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): if envelope or full: # Set ind size = domain_axes[axis].get_size() - d = np.arange(size) - ind = (d[value],) + ind = (np.arange(size)[value],) # Placeholder which will be overwritten later index = None @@ -491,7 +490,7 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): if full: # Set ind try: - index = normalize_cyclic_slice(index, size) + index = normalize_slice(index, size, cyclic=True) except IndexError: # Index is not a cyclic slice ind = (np.arange(size)[index],) @@ -747,73 +746,43 @@ def _point_not_in_cell(nodes_x, nodes_y, point): if halo: # Note: We're about to make 'indices' inconsistent # with 'ind', but that's OK because we're not - # going to use 'ind' again since 'create_mask' - # is False. + # going to use 'ind' again as 'create_mask' is + # False. for axis in item_axes: index = indices[axis] size = domain_axes[axis].get_size() try: # Is index a cyclic slice? - index = normalize_cyclic_slice(index, size) + index = normalize_slice(index, size, cyclic=True) except IndexError: - # Non-cyclic index - # - # Either a slice or a list/1-d array of - # int/bool. - # - # E.g. for halo=1 and size=5: - # slice(1, 3) -> [0, 1, 2, 3] - # slice(1, 4, 2) -> [0, 1, 3, 4] - # slice(2, 0, -1) -> [3, 2, 1, 0] - # [1, 2] -> [0, 1, 2, 3] - # [1, 3] -> [0, 1, 3, 4] - # [2, 1] -> [3, 2, 1, 0] - # [3, 1] -> [4, 3, 1, 0] - # [False, True, False, True, False] -> [0, 1, 3, 4] - index = normalize_index((index,), (size,))[0] - if isinstance(index, slice): - step = index.step - increasing = step is None or step > 0 - index = np.arange(size)[index] - else: - if is_dask_collection(index): - # Convert dask to numpy, and bool to int. - index = np.asanyarray(index) - index = normalize_index((index,), (size,))[0] - - increasing = index[0] <= index[-1] - - # Convert 'index' to a list (which will - # exclusively contain non-negative integers) - index = index.tolist() - - # Extend the list at each end - mn = min(index) - mx = max(index) - if increasing: - # "Increasing" sequence: Put smaller - # indices in the left hand halo, and - # larger ones in the right hand halo - left = range(max(*(0, mn - halo)), mn) - right = range(mx + 1, min(*(mx + 1 + halo, size))) + try: + index = normalize_slice(index, size) + except IndexError: + # Index is not a slice + cyclic = False else: - # "Decreasing" sequence: Put larger - # indices in the left hand halo, and - # smaller ones in the right hand halo - left = range(min(*(size - 1, mx + halo)), mx, -1) - right = range( - mn - 1, max(*(mn - 1 - halo, -1)), -1 + # Index is a non-cyclic slice, but if the + # axis is cyclic then it could become a + # cyclic slice once the halo is added. + cyclic = self.iscyclic(axis) + else: + # Index is a cyclic slice + cyclic = self.iscyclic(axis) + if not cyclic: + raise IndexError( + "Can't take a cyclic slice of a non-cyclic " + "axis" ) - index[:0] = left - index.extend(right) - else: - # Cyclic slice index + if cyclic: + # Cyclic slice, or potentially cyclic slice. # # E.g. for halo=1 and size=5: - # slice(-1, 2) -> slice(-2, 3) - # slice(1, -2, -1) -> slice(2, -3, -1) + # slice(0, 2) -> slice(-1, 3) + # slice(-1, 2) -> slice(-2, 3) + # slice(1, None, -1) -> slice(2, -2, -1) + # slice(1, -2, -1) -> slice(2, -3, -1) start = index.start stop = index.stop step = index.step @@ -825,11 +794,14 @@ def _point_not_in_cell(nodes_x, nodes_y, point): # would be required if abs(step) != 1. # Note that cyclic slices created by this # method will always have a step of 1. - raise ValueError( + raise IndexError( "A cyclic slice index can only have halos if " f"it has step 1 or -1. Got {index!r}" ) + if step < 0 and stop is None: + stop = -1 + if step > 0: # Increasing cyclic slice start = max(start - halo, stop - size) @@ -840,7 +812,88 @@ def _point_not_in_cell(nodes_x, nodes_y, point): stop = max(stop - halo, start - size) index = slice(start, stop, step) + else: + # A list/1-d array of int/bool, or a + # non-cyclic slice that can't beconme cyclic. + # + # E.g. for halo=1 and size=5: + # slice(1, 3) -> [0, 1, 2, 3] + # slice(1, 4, 2) -> [0, 1, 3, 4] + # slice(2, 0, -1) -> [3, 2, 1, 0] + # slice(2, 0, -1) -> [3, 2, 1, 0] + # [1, 2] -> [0, 1, 2, 3] + # [1, 3] -> [0, 1, 3, 4] + # [2, 1] -> [3, 2, 1, 0] + # [3, 1] -> [4, 3, 1, 0] + # [1, 3, 2] -> [0, 1, 3, 2, 1] + # [False, True, False, True, False] -> [0, 1, 3, 4] + if isinstance(index, slice): + index = np.arange(size)[index] + else: + if is_dask_collection(index): + index = np.asanyarray(index) + + index = normalize_index((index,), (size,))[0] + + # Find the left-most and right-most elements + # ('iL' and iR') of the sequence of positive + # integers, and whether the sequence is + # increasing or decreasing at each end + # ('increasing_L' and 'increasing_R') + # + # For instance: + # + # ------------ -- -- ------------ ------------ + # index iL iR increasing_L increasing_R + # ------------ -- -- ------------ ------------ + # [1, 2, 3, 4] 1 4 True True + # [4, 3, 2, 1] 4 1 False False + # [2, 1, 3, 4] 2 4 False True + # [1, 2, 4, 3] 1 3 True False + # [2, 2, 3, 4] 2 4 True True + # [2, 2, 4, 3] 2 3 True False + # [3, 3, 4, 4] 3 4 True True + # [4, 4, 3, 3] 4 3 False False + # [10] 10 10 True True + # ------------ -- -- ------------ ------------ + n_index = index.size + if n_index == 1: + iL = index[0] + iR = iL + increasing_L = True + increasing_R = True + elif n_index > 1: + iL = index[0] + iR = index[-1] + increasing_L = iL <= index[np.argmax(index != iL)] + increasing_R = ( + iR >= index[-1 - np.argmax(index[::-1] != iR)] + ) + else: + raise IndexError( + "Can't add a halo to a zero-sized index: " + f"{index}" + ) + + # Extend the list at each end, but not + # exceeding the axis limits. + if increasing_L: + left = range(max(*(0, iL - halo)), iL) + else: + left = range(min(*(size - 1, iL + halo)), iL, -1) + + if increasing_R: + right = range(iR + 1, min(*(iR + 1 + halo, size))) + else: + right = range( + iR - 1, max(*(iR - 1 - halo, -1)), -1 + ) + + index = index.tolist() + index[:0] = left + index.extend(right) + # Reset the returned index indices[axis] = index # Create an ancillary mask for these axes diff --git a/cf/test/test_Field.py b/cf/test/test_Field.py index 06eb40b594..bd9c239f4d 100644 --- a/cf/test/test_Field.py +++ b/cf/test/test_Field.py @@ -1561,26 +1561,25 @@ def test_Field_indices(self): self.assertEqual(g.construct("aux_x").array, 160) self.assertEqual(g.construct("aux_y").array, 3) - # Halos - x0 = f.dimension_coordinate("X").array - for index in (cf.wi(70, 200), slice(2, 6), [2, 3, 4, 5]): - indices = f.indices(0, grid_longitude=index) - g = f[indices] - self.assertEqual(g.shape, (1, 10, 4)) - x = g.dimension_coordinate("X").array - self.assertTrue((x == [80, 120, 160, 200]).all()) + # Halos: monotonic increasing sequence + index = [2, 3, 4, 5] + indices = f.indices(0, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 4)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == [80, 120, 160, 200]).all()) - indices = f.indices(1, grid_longitude=index) - g = f[indices] - self.assertEqual(g.shape, (1, 10, 6)) - x = g.dimension_coordinate("X").array - self.assertTrue((x == [40, 80, 120, 160, 200, 240]).all()) + indices = f.indices(1, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 6)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == [40, 80, 120, 160, 200, 240]).all()) - indices = f.indices(999, grid_longitude=index) - g = f[indices] - self.assertEqual(g.shape, (1, 10, 9)) - x = g.dimension_coordinate("X").array - self.assertTrue((x == x0).all()) + indices = f.indices(999, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 9)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == f.dimension_coordinate("X").array).all()) # Halos: non-monotonic sequence index = [2, 3, 4, 1] @@ -1594,13 +1593,36 @@ def test_Field_indices(self): g = f[indices] self.assertEqual(g.shape, (1, 10, 6)) x = g.dimension_coordinate("X").array - self.assertTrue((x == [200, 80, 120, 160, 40, 0]).all()) + self.assertTrue((x == [40, 80, 120, 160, 40, 0]).all()) - indices = f.indices(2, grid_longitude=index) - g = f[indices] - self.assertEqual(g.shape, (1, 10, 7)) - x = g.dimension_coordinate("X").array - self.assertTrue((x == [240, 200, 80, 120, 160, 40, 0]).all()) + for halo in (2, 999): + indices = f.indices(halo, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 7)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == [0, 40, 80, 120, 160, 40, 0]).all()) + + # Halos: cyclic slice increasing + for index in (cf.wi(70, 200), slice(2, 6)): + indices = f.indices(0, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 4)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == [80, 120, 160, 200]).all()) + + indices = f.indices(1, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 6)) + x = g.dimension_coordinate("X").array + self.assertTrue((x == [40, 80, 120, 160, 200, 240]).all()) + + indices = f.indices(999, grid_longitude=index) + g = f[indices] + self.assertEqual(g.shape, (1, 10, 9)) + x = g.dimension_coordinate("X").array + self.assertTrue( + (x == [-120, -80, -40, 0, 40, 80, 120, 160, 200]).all() + ) # Halos: cyclic slice increasing index = cf.wi(-170, 40) diff --git a/cf/test/test_functions.py b/cf/test/test_functions.py index 74596d0b93..daff8a485e 100644 --- a/cf/test/test_functions.py +++ b/cf/test/test_functions.py @@ -354,29 +354,41 @@ def test_size(self): def test_CFA(self): self.assertEqual(cf.CFA(), cf.__cfa_version__) - def test_normalize_cyclic_slice(self): + def test_normalize_slice(self): + self.assertEqual(cf.normalize_slice(slice(1, 4), 8), slice(1, 4, 1)) + self.assertEqual(cf.normalize_slice(slice(None), 8), slice(0, 8, 1)) + self.assertEqual( + cf.normalize_slice(slice(6, None, -1), 8), slice(6, None, -1) + ) + self.assertEqual(cf.normalize_slice(slice(-2, 4), 8), slice(6, 4, 1)) + # Cyclic slices self.assertEqual( - cf.normalize_cyclic_slice(slice(-2, 3), 8), slice(-2, 3, 1) + cf.normalize_slice(slice(-2, 3), 8, cyclic=True), slice(-2, 3, 1) ) self.assertEqual( - cf.normalize_cyclic_slice(slice(6, 3), 8), slice(-2, 3, 1) + cf.normalize_slice(slice(6, 3), 8, cyclic=True), slice(-2, 3, 1) ) self.assertEqual( - cf.normalize_cyclic_slice(slice(6, 3, 2), 8), slice(-2, 3, 2) + cf.normalize_slice(slice(6, 3, 2), 8, cyclic=True), slice(-2, 3, 2) ) self.assertEqual( - cf.normalize_cyclic_slice(slice(2, -3, -1), 8), slice(2, -3, -1) + cf.normalize_slice(slice(2, -3, -1), 8, cyclic=True), + slice(2, -3, -1), ) self.assertEqual( - cf.normalize_cyclic_slice(slice(2, 5, -1), 8), slice(2, -3, -1) + cf.normalize_slice(slice(2, 5, -1), 8, cyclic=True), + slice(2, -3, -1), ) self.assertEqual( - cf.normalize_cyclic_slice(slice(2, 5, -2), 8), slice(2, -3, -2) + cf.normalize_slice(slice(2, 5, -2), 8, cyclic=True), + slice(2, -3, -2), ) - # Non-cylic slices + with self.assertRaises(IndexError): + cf.normalize_slice([1, 2], 8) + for index in ( slice(1, 6), slice(6, 1, -1), @@ -386,7 +398,7 @@ def test_normalize_cyclic_slice(self): 5, ): with self.assertRaises(IndexError): - cf.normalize_cyclic_slice(index, 8) + cf.normalize_slice(index, 8, cyclic=True) if __name__ == "__main__": From e25bdc84f392ceb572d9c2f54683219dadeb4584 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Fri, 19 Apr 2024 09:29:55 +0100 Subject: [PATCH 13/30] 1-d index performance --- cf/field.py | 4 +++- cf/mixin/fielddomain.py | 17 ++++++++++++++--- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/cf/field.py b/cf/field.py index c4a583053f..8887469fa1 100644 --- a/cf/field.py +++ b/cf/field.py @@ -9052,7 +9052,9 @@ def indices(self, *mode, **kwargs): # spanned by the data if len(axis_indices) > len(data_axes): for axis, index in axis_indices.items(): - if axis in data_axes or index == slice(None): + if axis in data_axes or ( + isinstance(index, slice) and index == slice(None) + ): continue import dask.array as da diff --git a/cf/mixin/fielddomain.py b/cf/mixin/fielddomain.py index 2046f81669..6ddb275846 100644 --- a/cf/mixin/fielddomain.py +++ b/cf/mixin/fielddomain.py @@ -513,7 +513,15 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): logger.debug(" 1-d CASE 3:") # pragma: no cover index = item == value - index = index.to_dask_array() + + # Performance: Convert the 1-d 'index' to a numpy + # array of bool. + # + # This is beacuse Dask can be *very* slow at + # instantiation time when the 'index' is a Dask + # array, in which case contents of 'index' are + # unknown. + index = np.asanyarray(index) if envelope or full: # Set ind @@ -525,7 +533,10 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): # Placeholder which will be overwritten later index = None - + else: + # Convert bool to int, to save memory. + size = domain_axes[axis].get_size() + index = normalize_index(index, (size,))[0] else: raise ValueError( "Must specify a domain axis construct or a " @@ -833,7 +844,7 @@ def _point_not_in_cell(nodes_x, nodes_y, point): if is_dask_collection(index): index = np.asanyarray(index) - index = normalize_index((index,), (size,))[0] + index = normalize_index(index, (size,))[0] # Find the left-most and right-most elements # ('iL' and iR') of the sequence of positive From c9fc315d0779d07727595f31f50063da93774fdf Mon Sep 17 00:00:00 2001 From: David Hassell Date: Mon, 22 Apr 2024 16:42:09 +0100 Subject: [PATCH 14/30] Clarify docstring Co-authored-by: Sadie L. Bartholomew --- cf/docstring/docstring.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cf/docstring/docstring.py b/cf/docstring/docstring.py index 5cc0dd0d8e..d7b85279cb 100644 --- a/cf/docstring/docstring.py +++ b/cf/docstring/docstring.py @@ -636,7 +636,7 @@ ============== ====================================== *mode* Description ============== ====================================== - ```` If no positional arguments are + Not provided If no positional arguments are provided then assume the ``'compress'`` mode of operation with no halo added to the subspaced axes. From 64863272aad3179a0b44a0d572ee3bfdd77b20ee Mon Sep 17 00:00:00 2001 From: David Hassell Date: Mon, 22 Apr 2024 16:44:04 +0100 Subject: [PATCH 15/30] Typo Co-authored-by: Sadie L. Bartholomew --- cf/docstring/docstring.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cf/docstring/docstring.py b/cf/docstring/docstring.py index d7b85279cb..1f8cbbad8f 100644 --- a/cf/docstring/docstring.py +++ b/cf/docstring/docstring.py @@ -705,7 +705,7 @@ # subspace valid modes Field "{{subspace valid modes Field}}": """Valid modes are: - * ``'compress'`` This the default. + * ``'compress'`` This is the default. Unselected locations are removed to create the subspace. If the result is not hyperrectangular then the minimum amount of unselected locations From 472a570d6dfc9aa50df97d956e720a8c38e19267 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Mon, 22 Apr 2024 16:44:17 +0100 Subject: [PATCH 16/30] Typo Co-authored-by: Sadie L. Bartholomew --- cf/docstring/docstring.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cf/docstring/docstring.py b/cf/docstring/docstring.py index 1f8cbbad8f..3df5c6221e 100644 --- a/cf/docstring/docstring.py +++ b/cf/docstring/docstring.py @@ -729,7 +729,7 @@ # subspace valid modes Domain "{{subspace valid modes Domain}}": """Valid modes are: - * ``'compress'`` This the default. + * ``'compress'`` This is the default. Unselected locations are removed to create the subspace. If the result is not hyperrectangular then the minimum amount of unselected locations From 6c7fa9cb8aea275c1b5046f66fd81b20278602c6 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Mon, 22 Apr 2024 16:47:00 +0100 Subject: [PATCH 17/30] Typos Co-authored-by: Sadie L. Bartholomew --- cf/mixin/fielddomain.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cf/mixin/fielddomain.py b/cf/mixin/fielddomain.py index 6ddb275846..486d3a115d 100644 --- a/cf/mixin/fielddomain.py +++ b/cf/mixin/fielddomain.py @@ -218,7 +218,7 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): :Parameters: mode: `tuple` - The mode of operation and the halo size . See the + The mode of operation and the halo size. See the *mode* parameter of `indices` for details. data_axes: sequence of `str`, or `None` @@ -825,7 +825,7 @@ def _point_not_in_cell(nodes_x, nodes_y, point): index = slice(start, stop, step) else: # A list/1-d array of int/bool, or a - # non-cyclic slice that can't beconme cyclic. + # non-cyclic slice that can't become cyclic. # # E.g. for halo=1 and size=5: # slice(1, 3) -> [0, 1, 2, 3] From d70b77ae3e28bcd9a47b6ca8c0a3d075aa491bf9 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Mon, 22 Apr 2024 16:47:52 +0100 Subject: [PATCH 18/30] Update changelog Co-authored-by: Sadie L. Bartholomew --- Changelog.rst | 6 ------ 1 file changed, 6 deletions(-) diff --git a/Changelog.rst b/Changelog.rst index c76f1f3520..f4a10d80ae 100644 --- a/Changelog.rst +++ b/Changelog.rst @@ -6,12 +6,6 @@ version NEXT * Allow a halo to be added by `cf.Field.indices` and `cf.Field.subspace` (https://github.com/NCAS-CMS/cf-python/issues/759) - -version NEXT ------------- - -**2024-??-??** - * Added spherical regridding to discrete sampling geometry destination grids (https://github.com/NCAS-CMS/cf-python/issues/716) * Added 3-d spherical regridding to `cf.Field.regrids`, and the option From 47f051357a43d128b095fe54a40f87e8c22db92b Mon Sep 17 00:00:00 2001 From: David Hassell Date: Mon, 22 Apr 2024 16:49:45 +0100 Subject: [PATCH 19/30] Typo Co-authored-by: Sadie L. Bartholomew --- cf/field.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cf/field.py b/cf/field.py index 8887469fa1..f1d1469630 100644 --- a/cf/field.py +++ b/cf/field.py @@ -13292,7 +13292,7 @@ def subspace(self): {{subspace valid modes Field}} - In addition, an extra postional argument of ``'test'`` + In addition, an extra positional argument of ``'test'`` is allowed. When provided, the subspace is not returned, instead `True` or `False` is returned depending on whether or not it is possible for the From 252b644dc11d3be30d66c22011587e79d4ca6378 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Mon, 22 Apr 2024 16:50:15 +0100 Subject: [PATCH 20/30] Typo Co-authored-by: Sadie L. Bartholomew --- cf/functions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cf/functions.py b/cf/functions.py index bdcce6dbc6..a7f67bc3d3 100644 --- a/cf/functions.py +++ b/cf/functions.py @@ -2034,7 +2034,7 @@ def normalize_slice(index, size, cyclic=False): """Normalise a slice. If *index* is not a slice, or *cyclic* is True and *index* is not - a cylcic slice, then an `IndexError` is raised. + a cyclic slice, then an `IndexError` is raised. .. versionadded:: NEXTRELEASE From 15c34d126901a06e97b9afc6cf5413767955bbac Mon Sep 17 00:00:00 2001 From: David Hassell Date: Tue, 23 Apr 2024 09:32:47 +0100 Subject: [PATCH 21/30] Update cf/docstring/docstring.py Co-authored-by: Sadie L. Bartholomew --- cf/docstring/docstring.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cf/docstring/docstring.py b/cf/docstring/docstring.py index 3df5c6221e..c0014149fb 100644 --- a/cf/docstring/docstring.py +++ b/cf/docstring/docstring.py @@ -81,9 +81,9 @@ `cf.regrid_logging`. If it is enabled then logging takes place after every call. By default logging is disabled.""", # subspace halos - "{{subspace halos}}": """ - If a halo is defined via a positional argument, then each - subspaced axis will be extended to include that many extra + "{{subspace halos}}": """If a halo is defined via + a positional argument, then each subspaced axis + will be extended to include that many extra elements at each "side" of the axis. The number of extra elements will be automatically reduced if including the full amount defined by the halo would extend the subspace beyond From c1dc603bb585006156704c009e247cc4b1d7f277 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Tue, 23 Apr 2024 09:33:23 +0100 Subject: [PATCH 22/30] Update cf/docstring/docstring.py Co-authored-by: Sadie L. Bartholomew --- cf/docstring/docstring.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cf/docstring/docstring.py b/cf/docstring/docstring.py index c0014149fb..761c9f33a0 100644 --- a/cf/docstring/docstring.py +++ b/cf/docstring/docstring.py @@ -636,7 +636,7 @@ ============== ====================================== *mode* Description ============== ====================================== - Not provided If no positional arguments are + Not provided If no positional arguments are provided then assume the ``'compress'`` mode of operation with no halo added to the subspaced axes. From 6121429f4c4e76b2d3507bfafef7459b92957f82 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Tue, 23 Apr 2024 09:36:16 +0100 Subject: [PATCH 23/30] format --- cf/docstring/docstring.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/cf/docstring/docstring.py b/cf/docstring/docstring.py index 761c9f33a0..b904b6d010 100644 --- a/cf/docstring/docstring.py +++ b/cf/docstring/docstring.py @@ -81,9 +81,8 @@ `cf.regrid_logging`. If it is enabled then logging takes place after every call. By default logging is disabled.""", # subspace halos - "{{subspace halos}}": """If a halo is defined via - a positional argument, then each subspaced axis - will be extended to include that many extra + "{{subspace halos}}": """If a halo is defined via a positional argument, then each + subspaced axis will be extended to include that many extra elements at each "side" of the axis. The number of extra elements will be automatically reduced if including the full amount defined by the halo would extend the subspace beyond From 3e1e3c739d7a0b3538d5d817df195f8a8c2edc76 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Tue, 23 Apr 2024 09:45:02 +0100 Subject: [PATCH 24/30] subspace docstring --- cf/field.py | 2 - cf/subspacefield.py | 118 ++++++++++++++++++++++++++++++++------------ 2 files changed, 87 insertions(+), 33 deletions(-) diff --git a/cf/field.py b/cf/field.py index f1d1469630..fd6ea4f370 100644 --- a/cf/field.py +++ b/cf/field.py @@ -13270,8 +13270,6 @@ def subspace(self): ``3:-2:-1``) is assumed to "wrap" around, rather then producing a null result. - .. seealso:: `indices`, `squeeze`, `where`, `__getitem__` - **Halos** {{subspace halos}} diff --git a/cf/subspacefield.py b/cf/subspacefield.py index 6e6a554ba7..3505d756ea 100644 --- a/cf/subspacefield.py +++ b/cf/subspacefield.py @@ -2,7 +2,7 @@ class SubspaceField(mixin.Subspace): - """Create a subspace of a field construct. + """Create a subspace of the field construct. Creation of a new field construct which spans a subspace of the domain of an existing field construct is achieved either by @@ -70,37 +70,86 @@ class SubspaceField(mixin.Subspace): ``3:-2:-1``) is assumed to "wrap" around, rather then producing a null result. - .. seealso:: `cf.Field.indices`, `cf.Field.squeeze`, - `cf.Field.where` + **Halos** + + If a halo is defined via a positional argument, then each + subspaced axis will be extended to include that many extra + elements at each "side" of the axis. The number of extra elements + will be automatically reduced if including the full amount defined + by the halo would extend the subspace beyond the axis limits. + + For instance, ``f.subspace(X=slice(10, 20))`` will give identical + results to each of ``f.subspace(0, X=slice(10, 20))``, + ``f.subspace(1, X=slice(11, 19))``, ``f.subspace(2, X=slice(12, + 18))``, etc. + + .. seealso:: `cf.Field.indices`, `cf.Field.where`, + `cf.Field.__getitem__`, `cf.Field.__setitem__`, + `cf.Domain.subspace` :Parameters: - positional arguments: *optional* - There are three modes of operation, each of which provides - a different type of subspace: - - ============== ========================================== - *argument* Description - ============== ========================================== - ``'compress'`` This is the default mode. Unselected - locations are removed to create the - returned subspace. Note that if a - multi-dimensional metadata construct is - being used to define the indices then some - missing data may still be inserted at - unselected locations. - - ``'envelope'`` The returned subspace is the smallest that - contains all of the selected - indices. Missing data is inserted at - unselected locations within the envelope. - - ``'full'`` The returned subspace has the same domain - as the original field construct. Missing - data is inserted at unselected locations. - ============== ========================================== - - keyword parameters: *optional* + mode: optional + Specify the mode of operation (``mode``) and a halo to be + added to the subspaced axes (``halo``) with positional + arguments in format ``mode``, or ``halo``, or ``mode, + halo``, or with no positional arguments at all. + + A mode of operation is given as a `str`, and a halo as a + non-negative `int` (or any object that can be converted to + one): + + ============== ====================================== + *mode* Description + ============== ====================================== + Not provided If no positional arguments are + provided then assume the + ``'compress'`` mode of operation with + no halo added to the subspaced axes. + + ``mode`` Define the mode of operation with no + halo added to the subspaced axes. + + ``mode, halo`` Define a mode of operation, as well as + a halo to be added to the subspaced + axes. + + ``halo`` Assume the ``'compress'`` mode of + operation and define a halo to be + added to the subspaced axes. + ============== ====================================== + + Valid modes are: + + * ``'compress'`` This is the default. + Unselected locations are removed to create the + subspace. If the result is not hyperrectangular then + the minimum amount of unselected locations required + to make it so will also be specially + selected. Missing data is inserted at the specially + selected locations, unless a halo has been defined + (of any size, including 0). + + * ``'envelope'`` + The subspace is the smallest hyperrectangular + subspace that contains all of the selected + locations. Missing data is inserted at unselected + locations within the envelope, unless a halo has been + defined (of any size, including 0). + + * ``'full'`` + The subspace has the same domain as the original + construct. Missing data is inserted at unselected + locations, unless a halo has been defined (of any + size, including 0). + + In addition, an extra positional argument of ``'test'`` is + allowed. When provided, the subspace is not returned, + instead `True` or `False` is returned depending on whether + or not it is possible for the requested subspace to be + created. + + keyword parameters: optional A keyword name is an identity of a metadata construct, and the keyword value provides a condition for inferring indices that apply to the dimension (or dimensions) @@ -110,9 +159,12 @@ class SubspaceField(mixin.Subspace): :Returns: - `Field` + `Field` or `bool` An independent field construct containing the subspace of - the original field. + the original field. If the ``'test'`` positional argument + has been set then return `True` or `False` depending on + whether or not it is possible to create specified + subspace. **Examples** @@ -190,6 +242,10 @@ def __call__(self, *args, **kwargs): acting along orthogonal dimensions, some missing data may still need to be inserted into the field construct's data. + **Halos** + + {{subspace halos}} + .. seealso:: `cf.Field.indices`, `cf.Field.squeeze`, `cf.Field.where` From 3a092e646fee2fbd7ec63ce211a2eb0ab283c0c3 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Tue, 23 Apr 2024 09:54:01 +0100 Subject: [PATCH 25/30] size 0 halo note --- cf/docstring/docstring.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/cf/docstring/docstring.py b/cf/docstring/docstring.py index b904b6d010..69b2df70ac 100644 --- a/cf/docstring/docstring.py +++ b/cf/docstring/docstring.py @@ -635,7 +635,7 @@ ============== ====================================== *mode* Description ============== ====================================== - Not provided If no positional arguments are + Not provided If no positional arguments are provided then assume the ``'compress'`` mode of operation with no halo added to the subspaced axes. @@ -724,7 +724,14 @@ The subspace has the same domain as the original construct. Missing data is inserted at unselected locations, unless a halo has been defined (of any - size, including 0).""", + size, including 0). + + .. note:: Setting a halo size of `0` differs from not + not defining a halo at all. The shape of the + returned field will always be the same, but + in the former case missing data will not be + inserted at unselected locations (if any) + within the output domain.""", # subspace valid modes Domain "{{subspace valid modes Domain}}": """Valid modes are: From cb535c45235dd525d9414d270875974ab7f1a361 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Tue, 23 Apr 2024 10:14:16 +0100 Subject: [PATCH 26/30] subspace mode -> config --- cf/docstring/docstring.py | 14 +++++++------- cf/domain.py | 20 ++++++++++---------- cf/field.py | 10 +++++----- cf/mixin/fielddomain.py | 22 +++++++++++----------- cf/subspacefield.py | 23 ++++++++++++----------- 5 files changed, 45 insertions(+), 44 deletions(-) diff --git a/cf/docstring/docstring.py b/cf/docstring/docstring.py index 69b2df70ac..10d9d5c87a 100644 --- a/cf/docstring/docstring.py +++ b/cf/docstring/docstring.py @@ -620,13 +620,13 @@ "{{to_size: `int`, optional}}": """to_size: `int`, optional Pad the axis after so that the new axis has the given size.""", - # subspace mode options - "{{mode: optional}}": """mode: optional - Specify the mode of operation (``mode``) and a halo to - be added to the subspaced axes (``halo``) with - positional arguments in format ``mode``, or ``halo``, - or ``mode, halo``, or with no positional arguments at - all. + # subspace config options + "{{config: optional}}": """config: optional + Configure the subspace by specifying the mode of + operation (``mode``) and any halo to be added to the + subspaced axes (``halo``), with positional arguments + in the format ``mode``, or ``halo``, or ``mode, + halo``, or with no positional arguments at all. A mode of operation is given as a `str`, and a halo as a non-negative `int` (or any object that can be diff --git a/cf/domain.py b/cf/domain.py index 2824d07ebf..1fbdb73e1c 100644 --- a/cf/domain.py +++ b/cf/domain.py @@ -739,7 +739,7 @@ def identities(self): return out - def indices(self, *mode, **kwargs): + def indices(self, *config, **kwargs): """Create indices that define a subspace of the domain construct. @@ -789,7 +789,7 @@ def indices(self, *mode, **kwargs): :Parameters: - {{mode: optional}} + {{config: optional}} {{subspace valid modes Domain}} @@ -847,7 +847,7 @@ def indices(self, *mode, **kwargs): """ # Get the indices for every domain axis in the domain, without # any auxiliary masks. - domain_indices = self._indices(mode, None, False, kwargs) + domain_indices = self._indices(config, None, False, kwargs) return domain_indices["indices"] @@ -1094,7 +1094,7 @@ def roll(self, axis, shift, inplace=False): return d - def subspace(self, *mode, **kwargs): + def subspace(self, *config, **kwargs): """Create a subspace of the field construct. Creation of a new domain construct which spans a subspace of @@ -1141,7 +1141,7 @@ def subspace(self, *mode, **kwargs): :Parameters: - {{mode: optional}} + {{config: optional}} {{subspace valid modes Domain}} @@ -1186,19 +1186,19 @@ def subspace(self, *mode, **kwargs): """ test = False - if "test" in mode: - mode = list(mode) - mode.remove("test") + if "test" in config: + config = list(config) + config.remove("test") test = True - if not mode and not kwargs: + if not config and not kwargs: if test: return True return self.copy() try: - indices = self.indices(*mode, **kwargs) + indices = self.indices(*config, **kwargs) except ValueError as error: if test: return False diff --git a/cf/field.py b/cf/field.py index fd6ea4f370..12c8c885f6 100644 --- a/cf/field.py +++ b/cf/field.py @@ -8809,7 +8809,7 @@ def insert_dimension( inplace=inplace, ) - def indices(self, *mode, **kwargs): + def indices(self, *config, **kwargs): """Create indices that define a subspace of the field construct. The subspace is defined by identifying indices based on the @@ -8879,7 +8879,7 @@ def indices(self, *mode, **kwargs): :Parameters: - {{mode: optional}} + {{config: optional}} {{subspace valid modes Field}} @@ -9006,7 +9006,7 @@ def indices(self, *mode, **kwargs): [-- -- -- -- -- -- 270.6 273.0 270.6]]] """ - if "exact" in mode: + if "exact" in config: _DEPRECATION_ERROR_ARG( self, "indices", @@ -9020,7 +9020,7 @@ def indices(self, *mode, **kwargs): # Get the indices for every domain axis in the domain, # including any ancillary masks - domain_indices = self._indices(mode, data_axes, True, kwargs) + domain_indices = self._indices(config, data_axes, True, kwargs) # Initialise the output indices with any ancillary masks. # Ensure that each ancillary mask is broadcastable to the @@ -13286,7 +13286,7 @@ def subspace(self): :Parameters: - {{mode: optional}} + {{config: optional}} {{subspace valid modes Field}} diff --git a/cf/mixin/fielddomain.py b/cf/mixin/fielddomain.py index 486d3a115d..a187d5df46 100644 --- a/cf/mixin/fielddomain.py +++ b/cf/mixin/fielddomain.py @@ -203,7 +203,7 @@ def _equivalent_coordinate_references( return True - def _indices(self, mode, data_axes, ancillary_mask, kwargs): + def _indices(self, config, data_axes, ancillary_mask, kwargs): """Create indices that define a subspace of the field or domain construct. @@ -217,9 +217,9 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): :Parameters: - mode: `tuple` + config: `tuple` The mode of operation and the halo size. See the - *mode* parameter of `indices` for details. + *config* parameter of `indices` for details. data_axes: sequence of `str`, or `None` The domain axis identifiers of the data axes, or @@ -255,24 +255,24 @@ def _indices(self, mode, data_axes, ancillary_mask, kwargs): debug = is_log_level_debug(logger) # Parse mode and halo - n_mode = len(mode) - if not n_mode: + n_config = len(config) + if not n_config: mode = None halo = None - elif n_mode == 1: + elif n_config == 1: try: - halo = int(mode[0]) + halo = int(config[0]) except ValueError: - mode = mode[0] + mode = config[0] halo = None else: mode = None - elif n_mode == 2: - mode, halo = mode + elif n_config == 2: + mode, halo = config else: raise ValueError( "Can't provide more than two positional arguments. " - f"Got: {', '.join(repr(x) for x in mode)}" + f"Got: {', '.join(repr(x) for x in config)}" ) compress = mode is None or mode == "compress" diff --git a/cf/subspacefield.py b/cf/subspacefield.py index 3505d756ea..926c5627f8 100644 --- a/cf/subspacefield.py +++ b/cf/subspacefield.py @@ -89,11 +89,12 @@ class SubspaceField(mixin.Subspace): :Parameters: - mode: optional - Specify the mode of operation (``mode``) and a halo to be - added to the subspaced axes (``halo``) with positional - arguments in format ``mode``, or ``halo``, or ``mode, - halo``, or with no positional arguments at all. + config: optional + Configure the subspace by specifying the mode of operation + (``mode``) and any halo to be added to the subspaced axes + (``halo``), with positional arguments in the format + ``mode``, or ``halo``, or ``mode, halo``, or with no + positional arguments at all. A mode of operation is given as a `str`, and a halo as a non-negative `int` (or any object that can be converted to @@ -201,7 +202,7 @@ class SubspaceField(mixin.Subspace): __slots__ = [] - def __call__(self, *args, **kwargs): + def __call__(self, *config, **kwargs): """Create a subspace of a field construct. Creation of a new field construct which spans a subspace of the @@ -324,19 +325,19 @@ def __call__(self, *args, **kwargs): field = self.variable test = False - if "test" in args: - args = list(args) - args.remove("test") + if "test" in config: + config = list(config) + config.remove("test") test = True - if not args and not kwargs: + if not config and not kwargs: if test: return True return field.copy() try: - indices = field.indices(*args, **kwargs) + indices = field.indices(*config, **kwargs) out = field[indices] except (ValueError, IndexError) as error: if test: From f038dd8e60c84911d9fbf3990e469eba988f533b Mon Sep 17 00:00:00 2001 From: David Hassell Date: Tue, 23 Apr 2024 10:50:41 +0100 Subject: [PATCH 27/30] warn when reducing halo size --- cf/mixin/fielddomain.py | 59 ++++++++++++++++++++++++++++++++++------- 1 file changed, 49 insertions(+), 10 deletions(-) diff --git a/cf/mixin/fielddomain.py b/cf/mixin/fielddomain.py index a187d5df46..9cd3a1a680 100644 --- a/cf/mixin/fielddomain.py +++ b/cf/mixin/fielddomain.py @@ -759,6 +759,7 @@ def _point_not_in_cell(nodes_x, nodes_y, point): # with 'ind', but that's OK because we're not # going to use 'ind' again as 'create_mask' is # False. + reduced_halo = False for axis in item_axes: index = indices[axis] size = domain_axes[axis].get_size() @@ -815,12 +816,27 @@ def _point_not_in_cell(nodes_x, nodes_y, point): if step > 0: # Increasing cyclic slice - start = max(start - halo, stop - size) - stop = min(stop + halo, size + start) + start = start - halo + if start < stop - size: + start = stop - size + reduced_halo = True + + stop = stop + halo + if stop > size + start: + stop = size + start + reduced_halo = True + else: # Decreasing cyclic slice - start = min(start + halo, size + stop) - stop = max(stop - halo, start - size) + start = start + halo + if start > size + stop: + start = size + stop + reduced_halo = True + + stop = stop - halo + if stop < start - size: + stop = start - size + reduced_halo = True index = slice(start, stop, step) else: @@ -889,16 +905,34 @@ def _point_not_in_cell(nodes_x, nodes_y, point): # Extend the list at each end, but not # exceeding the axis limits. if increasing_L: - left = range(max(*(0, iL - halo)), iL) + start = iL - halo + if start < 0: + start = 0 + reduced_halo = True + + left = range(start, iL) else: - left = range(min(*(size - 1, iL + halo)), iL, -1) + start = iL + halo + if start > size - 1: + start = size - 1 + reduced_halo = True + + left = range(start, iL, -1) if increasing_R: - right = range(iR + 1, min(*(iR + 1 + halo, size))) + stop = iR + 1 + halo + if stop > size: + stop = size + reduced_halo = True + + right = range(iR + 1, stop) else: - right = range( - iR - 1, max(*(iR - 1 - halo, -1)), -1 - ) + stop = iR - 1 - halo + if stop < -1: + stop = -1 + reduced_halo = True + + right = range(iR - 1, stop, -1) index = index.tolist() index[:0] = left @@ -907,6 +941,11 @@ def _point_not_in_cell(nodes_x, nodes_y, point): # Reset the returned index indices[axis] = index + if reduced_halo: + logger.warning( + "Halo reduced to keep subspace within axis limits" + ) + # Create an ancillary mask for these axes if debug: logger.debug( From 0a8b8a067f066e771f526b85ce8fabdfef8fc955 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Tue, 23 Apr 2024 10:59:46 +0100 Subject: [PATCH 28/30] test indices API --- cf/test/test_Field.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/cf/test/test_Field.py b/cf/test/test_Field.py index bd9c239f4d..eef39732e3 100644 --- a/cf/test/test_Field.py +++ b/cf/test/test_Field.py @@ -1679,6 +1679,14 @@ def test_Field_indices(self): g = f[indices] self.assertFalse(np.ma.is_masked(g.array)) + # Test API with 0/1/2 arguments + kwargs = {"grid_latitude": [1]} + i = f.indices(**kwargs) + j = f.indices(0, **kwargs) + k = f.indices("compress", 0, **kwargs) + self.assertEqual(i, j) + self.assertEqual(i, k) + # Subspace has size 0 axis resulting from dask array index indices = f.indices(grid_latitude=cf.contains(-23.2)) with self.assertRaises(IndexError): From 31ec154bd7a2e518782c9018ea2cafe06c5fad41 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Tue, 23 Apr 2024 14:32:21 +0100 Subject: [PATCH 29/30] susbapce API test --- cf/test/test_Field.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/cf/test/test_Field.py b/cf/test/test_Field.py index eef39732e3..4dc2cab40d 100644 --- a/cf/test/test_Field.py +++ b/cf/test/test_Field.py @@ -2749,6 +2749,14 @@ def test_Field_subspace(self): h = f.subspace(grid_longitude=np.float64(20)) self.assertTrue(g.equals(h)) + # Test API with 0/1/2 arguments + kwargs = {"grid_latitude": [1]} + i = f.subspace(**kwargs) + j = f.subspace(0, **kwargs) + k = f.subspace("compress", 0, **kwargs) + self.assertEqual(i, j) + self.assertEqual(i, k) + def test_Field_auxiliary_to_dimension_to_auxiliary(self): f = cf.example_field(0) nd = len(f.dimension_coordinates()) From 7d9c32faa7c3e6fb8845d6b18ddb127248a00175 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Wed, 24 Apr 2024 15:41:09 +0100 Subject: [PATCH 30/30] add missing docs --- cf/subspacefield.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/cf/subspacefield.py b/cf/subspacefield.py index 926c5627f8..b9e1cf25d5 100644 --- a/cf/subspacefield.py +++ b/cf/subspacefield.py @@ -144,6 +144,13 @@ class SubspaceField(mixin.Subspace): locations, unless a halo has been defined (of any size, including 0). + .. note:: Setting a halo size of `0` differs from not not + defining a halo at all. The shape of the + returned field will always be the same, but in + the former case missing data will not be + inserted at unselected locations (if any) within + the output domain. + In addition, an extra positional argument of ``'test'`` is allowed. When provided, the subspace is not returned, instead `True` or `False` is returned depending on whether