From 094a60a3f329b74ce3323e7667cf25b35e18d8ea Mon Sep 17 00:00:00 2001 From: David Hassell Date: Wed, 20 Sep 2023 23:06:52 +0100 Subject: [PATCH] dev --- Changelog.rst | 4 +- cf/__init__.py | 2 +- cf/data/dask_regrid.py | 19 +- cf/data/data.py | 25 +- cf/field.py | 2966 ++++++++++++++++++------------------ cf/regrid/regrid.py | 165 +- cf/test/test_general.py | 23 - cf/test/test_read_write.py | 5 +- cf/weights.py | 1428 +++++++++++++++++ 9 files changed, 3085 insertions(+), 1552 deletions(-) create mode 100644 cf/weights.py diff --git a/Changelog.rst b/Changelog.rst index c216a15d7a..fb81561a34 100644 --- a/Changelog.rst +++ b/Changelog.rst @@ -256,8 +256,8 @@ version 3.11.0 * Fix for `cf.aggregate` failures when a datum or coordinate conversion parameter has an array value (https://github.com/NCAS-CMS/cf-python/issues/230) -* Allow for regridding using a destination field featuring size 1 dimension(s) - (https://github.com/NCAS-CMS/cf-python/issues/250) +* Allow for regridding using a destination field featuring size 1 + dimension(s) (https://github.com/NCAS-CMS/cf-python/issues/250) * Fix bug that sometimes caused `cf.Field.autocyclic` to fail when setting a construct that is cyclic and has a defined period * Fix bug that sometimes caused a failure when reading PP extra data diff --git a/cf/__init__.py b/cf/__init__.py index 3cbcd13cf2..f22835065b 100644 --- a/cf/__init__.py +++ b/cf/__init__.py @@ -75,7 +75,7 @@ __Conventions__ = "CF-1.11" __date__ = "2023-??-??" -__version__ = "3.16.0" +__version__ = "3.17.0" _requires = ( "numpy", diff --git a/cf/data/dask_regrid.py b/cf/data/dask_regrid.py index 4a60d261d3..a1a9727561 100644 --- a/cf/data/dask_regrid.py +++ b/cf/data/dask_regrid.py @@ -287,18 +287,25 @@ def regrid( a = a.reshape(non_regrid_shape + tuple(dst_shape)) n_dst_axes = len(dst_shape) + if n_src_axes == n_dst_axes: pass elif n_src_axes == 1 and n_dst_axes > 1: # E.g. UGRID regridded to regular lat-lon; changes # 'axis_order' from [0,2,1] to [0,3,1,2] - r = axis_order[-1] - axis_order = [i + n_dst_axes - 1 if i > r else i for i in axis_order] - axis_order[-1:] = range(r, r + n_dst_axes) - elif n_dst_axes == 1 and n_src_axes > 1: + raxis = axis_order[-1] + axis_order = [ + i if i <= raxis else i + n_dst_axes - 1 for i in axis_order + ] + axis_order[-1:] = range(raxis, raxis + n_dst_axes) + elif n_src_axes == 2 and n_dst_axes == 1: # E.g. regular lat-lon regridded to UGRID; changes - # 'axis_order' from [0,3,2,1] to [0,2,1] - pass # TODOUGRID + # 'axis_order' from [0,2,4,5,1,3] to [0,2,3,4,1], or + # [0,2,4,5,3,1] to [0,1,3,4,2] + raxis0, raxis = axis_order[-2:] + print(axis_order) + axis_order = [i if i <= raxis else i - 1 for i in axis_order[:-1]] + print(axis_order) else: raise ValueError("TODOUGRID") diff --git a/cf/data/data.py b/cf/data/data.py index 2377e2686b..90b340ab16 100644 --- a/cf/data/data.py +++ b/cf/data/data.py @@ -3734,6 +3734,7 @@ def _regrid( from .dask_regrid import regrid, regrid_weights shape = self.shape + ndim = self.ndim src_shape = tuple(shape[i] for i in regrid_axes) if src_shape != operator.src_shape: raise ValueError( @@ -3755,12 +3756,17 @@ def _regrid( # Define the regridded chunksizes regridded_chunks = [] # The 'chunks' parameter to `map_blocks` + drop_axis = [] # The 'drop_axis' parameter to `map_blocks` new_axis = [] # The 'new_axis' parameter to `map_blocks` n = 0 for i, c in enumerate(dx.chunks): if i in regridded_sizes: sizes = regridded_sizes[i] n_sizes = len(sizes) + if not n_sizes: + drop_axis.append(i) + continue + regridded_chunks.extend(sizes) if n_sizes > 1: new_axis.extend(range(n + 1, n + n_sizes)) @@ -3770,17 +3776,21 @@ def _regrid( n += 1 + # Update the axis identifiers. + # + # This is necessary when regridding changes the number of data + # dimensions (e.g. as happens when regridding a mesh topology + # axis to/from separate lat and lon axes). if new_axis: - # Update the axis identifiers. - # - # This is necessary when regridding changes the number of - # data dimensions (e.g. as happens when regridding a mesh - # topology axis to separate lat and lon axes). axes = list(self._axes) for i in new_axis: axes.insert(i, new_axis_identifier(tuple(axes))) - self._axes = tuple(axes) + self._axes = axes + elif drop_axis: + axes = self._axes + axes = [axes[i] for i in range(ndim) if i not in drop_axis] + self._axes = axes # Set the output data type if method in ("nearest_dtos", "nearest_stod"): @@ -3788,7 +3798,7 @@ def _regrid( else: dst_dtype = float - non_regrid_axes = [i for i in range(self.ndim) if i not in regrid_axes] + non_regrid_axes = [i for i in range(ndim) if i not in regrid_axes] src_mask = operator.src_mask if src_mask is not None: @@ -3813,6 +3823,7 @@ def _regrid( weights_dst_mask=weights_dst_mask, ref_src_mask=src_mask, chunks=regridded_chunks, + drop_axis=drop_axis, new_axis=new_axis, meta=np.array((), dtype=dst_dtype), ) diff --git a/cf/field.py b/cf/field.py index b20a41b42d..50b1fd891b 100644 --- a/cf/field.py +++ b/cf/field.py @@ -1834,1433 +1834,1431 @@ def _equivalent_construct_data( # ---------------------------------------------------------------- # Worker functions for weights # ---------------------------------------------------------------- - def _weights_area_XY( - self, - comp, - weights_axes, - auto=False, - measure=False, - radius=None, - methods=False, - ): - """Calculate area weights from X and Y dimension coordinate - constructs. - - :Parameters: - - measure: `bool` - If True then make sure that the weights represent true - cell sizes. - - methods: `bool`, optional - If True then add a description of the method used to - create the weights to the *comp* dictionary, as opposed to - the actual weights. - - :Returns: - - `bool` or `None` - - """ - xkey, xcoord = self.dimension_coordinate( - "X", item=True, default=(None, None) - ) - ykey, ycoord = self.dimension_coordinate( - "Y", item=True, default=(None, None) - ) - - if xcoord is None or ycoord is None: - if auto: - return - - raise ValueError( - "No unique 'X' and 'Y' dimension coordinate constructs for " - "calculating area weights" - ) - - if xcoord.Units.equivalent( - Units("radians") - ) and ycoord.Units.equivalent(Units("radians")): - pass - elif xcoord.Units.equivalent( - Units("metres") - ) and ycoord.Units.equivalent(Units("metres")): - radius = None - else: - if auto: - return - - raise ValueError( - "Insufficient coordinate constructs for calculating " - "area weights" - ) - - xaxis = self.get_data_axes(xkey)[0] - yaxis = self.get_data_axes(ykey)[0] - - for axis in (xaxis, yaxis): - if axis in weights_axes: - if auto: - return - - raise ValueError( - "Multiple weights specifications for " - f"{self.constructs.domain_axis_identity(axis)!r} axis" - ) - - if measure and radius is not None: - radius = self.radius(default=radius) - - if measure or xcoord.size > 1: - if not xcoord.has_bounds(): - if auto: - return - - raise ValueError( - "Can't create area weights: No bounds for " - f"{xcoord.identity()!r} axis" - ) - - if methods: - comp[(xaxis,)] = "linear " + xcoord.identity() - else: - cells = xcoord.cellsize - if xcoord.Units.equivalent(Units("radians")): - cells.Units = _units_radians - if measure: - cells *= radius - cells.override_units(radius.Units, inplace=True) - else: - cells.Units = Units("metres") - - comp[(xaxis,)] = cells - - weights_axes.add(xaxis) - - if measure or ycoord.size > 1: - if not ycoord.has_bounds(): - if auto: - return - - raise ValueError( - "Can't create area weights: No bounds for " - f"{ycoord.identity()!r} axis" - ) - - if ycoord.Units.equivalent(Units("radians")): - ycoord = ycoord.clip(-90, 90, units=Units("degrees")) - ycoord.sin(inplace=True) - - if methods: - comp[(yaxis,)] = "linear sine " + ycoord.identity() - else: - cells = ycoord.cellsize - if measure: - cells *= radius - - comp[(yaxis,)] = cells - else: - if methods: - comp[(yaxis,)] = "linear " + ycoord.identity() - else: - cells = ycoord.cellsize - comp[(yaxis,)] = cells - - weights_axes.add(yaxis) - - return True - - def _weights_data( - self, - w, - comp, - weights_axes, - axes=None, - data=False, - components=False, - methods=False, - ): - """Creates weights for the field construct's data array. - - :Parameters: - - methods: `bool`, optional - If True then add a description of the method used to - create the weights to the *comp* dictionary, as - opposed to the actual weights. - - """ - # -------------------------------------------------------- - # Data weights - # -------------------------------------------------------- - field_data_axes = self.get_data_axes() - - if axes is not None: - if isinstance(axes, (str, int)): - axes = (axes,) - - if len(axes) != w.ndim: - raise ValueError( - "'axes' parameter must provide an axis identifier " - f"for each weights data dimension. Got {axes!r} for " - f"{w.ndim} dimension(s)." - ) - - iaxes = [ - field_data_axes.index(self.domain_axis(axis, key=True)) - for axis in axes - ] - - if data: - for i in range(self.ndim): - if i not in iaxes: - w = w.insert_dimension(position=i) - iaxes.insert(i, i) - - w = w.transpose(iaxes) - - if w.ndim > 0: - while w.shape[0] == 1: - w = w.squeeze(0) - - else: - iaxes = list(range(self.ndim - w.ndim, self.ndim)) - - if not (components or methods): - if not self._is_broadcastable(w.shape): - raise ValueError( - f"The 'Data' weights (shape {w.shape}) are not " - "broadcastable to the field construct's data " - f"(shape {self.shape})." - ) - - axes0 = field_data_axes[self.ndim - w.ndim :] - else: - axes0 = [field_data_axes[i] for i in iaxes] - - for axis0 in axes0: - if axis0 in weights_axes: - raise ValueError( - "Multiple weights specified for " - f"{self.constructs.domain_axis_identity(axis0)!r} axis" - ) - - if methods: - comp[tuple(axes0)] = "custom data" - else: - comp[tuple(axes0)] = w - - weights_axes.update(axes0) - - return True - - def _weights_field(self, fields, comp, weights_axes, methods=False): - """Creates a weights field.""" - s = self.analyse_items() - - domain_axes = self.domain_axes(todict=True) - for w in fields: - t = w.analyse_items() - # TODO CHECK this with org - - if t["undefined_axes"]: - w_domain_axes_1 = w.domain_axes( - filter_by_size=(1,), todict=True - ) - if set(w_domain_axes_1).intersection(t["undefined_axes"]): - raise ValueError( - f"Weights field {w} has at least one undefined " - f"domain axes: {w_domain_axes_1}." - ) - - w = w.squeeze() - - w_domain_axes = w.domain_axes(todict=True) - - axis1_to_axis0 = {} - - coordinate_references = self.coordinate_references(todict=True) - w_coordinate_references = w.coordinate_references(todict=True) - - for axis1 in w.get_data_axes(): - identity = t["axis_to_id"].get(axis1, None) - if identity is None: - raise ValueError( - "Weights field has unmatched, size > 1 " - f"{w.constructs.domain_axis_identity(axis1)!r} axis" - ) - - axis0 = s["id_to_axis"].get(identity, None) - if axis0 is None: - raise ValueError( - f"Weights field has unmatched, size > 1 {identity!r} " - "axis" - ) - - w_axis_size = w_domain_axes[axis1].get_size() - self_axis_size = domain_axes[axis0].get_size() - - if w_axis_size != self_axis_size: - raise ValueError( - f"Weights field has incorrectly sized {identity!r} " - f"axis ({w_axis_size} != {self_axis_size})" - ) - - axis1_to_axis0[axis1] = axis0 - - # Check that the defining coordinate data arrays are - # compatible - key0 = s["axis_to_coord"][axis0] - key1 = t["axis_to_coord"][axis1] - - if not self._equivalent_construct_data( - w, key0=key0, key1=key1, s=s, t=t - ): - raise ValueError( - f"Weights field has incompatible {identity!r} " - "coordinates" - ) - - # Still here? Then the defining coordinates have - # equivalent data arrays - - # If the defining coordinates are attached to - # coordinate references then check that those - # coordinate references are equivalent - refs0 = [ - key - for key, ref in coordinate_references.items() - if key0 in ref.coordinates() - ] - refs1 = [ - key - for key, ref in w_coordinate_references.items() - if key1 in ref.coordinates() - ] - - nrefs = len(refs0) - if nrefs > 1 or nrefs != len(refs1): - # The defining coordinate are associated with - # different numbers of coordinate references - equivalent_refs = False - elif not nrefs: - # Neither defining coordinate is associated with a - # coordinate reference - equivalent_refs = True - else: - # Each defining coordinate is associated with - # exactly one coordinate reference - equivalent_refs = self._equivalent_coordinate_references( - w, key0=refs0[0], key1=refs1[0], s=s, t=t - ) - - if not equivalent_refs: - raise ValueError( - "Input weights field has an incompatible " - "coordinate reference" - ) - - axes0 = tuple( - [axis1_to_axis0[axis1] for axis1 in w.get_data_axes()] - ) - - for axis0 in axes0: - if axis0 in weights_axes: - raise ValueError( - "Multiple weights specified for " - f"{self.constructs.domain_axis_identity(axis0)!r} " - "axis" - ) - - comp[tuple(axes0)] = w.data - - weights_axes.update(axes0) - - return True - - def _weights_field_scalar(self, methods=False): - """Return a scalar field of weights with long_name ``'weight'``. - - :Returns: - - `Field` - - """ - data = Data(1.0, "1") - - f = type(self)() - f.set_data(data, copy=False) - f.long_name = "weight" - f.comment = f"Weights for {self!r}" - - return f - - def _weights_geometry_area( - self, - domain_axis, - comp, - weights_axes, - auto=False, - measure=False, - radius=None, - great_circle=False, - return_areas=False, - methods=False, - ): - """Creates area weights for polygon geometry cells. - - .. versionadded:: 3.2.0 - - :Parameters: - - domain_axis : `str` or `None` - - measure: `bool` - If True then make sure that the weights represent true - cell sizes. - - methods: `bool`, optional - If True then add a description of the method used to - create the weights to the *comp* dictionary, as opposed to - the actual weights. - - :Returns: - - `bool` or `Data` - - """ - axis, aux_X, aux_Y, aux_Z = self._weights_yyy( - domain_axis, "polygon", methods=methods, auto=auto - ) - - if axis is None: - if auto: - return False - - if domain_axis is None: - raise ValueError("No polygon geometries") - - raise ValueError( - "No polygon geometries for " - f"{self.constructs.domain_axis_identity(domain_axis)!r} axis" - ) - - if axis in weights_axes: - if auto: - return False - - raise ValueError( - "Multiple weights specifications for " - f"{self.constructs.domain_axis_identity(axis)!r} axis" - ) - - ugrid = self.domain_topologies(filter_by_axis=(axis,), todict=True) - - # Check for interior rings - if ugrid: - interior_ring = None - else: - interior_ring_X = aux_X.get_interior_ring(None) - interior_ring_Y = aux_Y.get_interior_ring(None) - if interior_ring_X is None and interior_ring_Y is None: - interior_ring = None - elif interior_ring_X is None: - raise ValueError( - "Can't find weights: X coordinates have missing " - "interior ring variable" - ) - elif interior_ring_Y is None: - raise ValueError( - "Can't find weights: Y coordinates have missing " - "interior ring variable" - ) - elif not interior_ring_X.data.equals(interior_ring_Y.data): - raise ValueError( - "Can't find weights: Interior ring variables for X and Y " - "coordinates have different data values" - ) - else: - interior_ring = interior_ring_X.data - if interior_ring.shape != aux_X.bounds.shape[:-1]: - raise ValueError( - "Can't find weights: Interior ring variables have " - f"incorrect shape. Got {interior_ring.shape}, " - f"expected {aux_X.bounds.shape[:-1]}" - # TODOUGRID: is this not tested in set_interior_ring? - ) - - x = aux_X.bounds.data - y = aux_Y.bounds.data - spherical = False - - if x.Units.equivalent(_units_metres) and y.Units.equivalent( - _units_metres - ): - # -------------------------------------------------------- - # Plane polygons defined by straight lines. - # - # The area is given by the shoelace formula: - # https://en.wikipedia.org/wiki/Shoelace_formula - # - # Do this in preference to weights based on spherical - # polygons, which require the great circle assumption. - # -------------------------------------------------------- - if methods: - comp[(axis,)] = "area plane polygon geometry" - return True - - y.Units = x.Units - - all_areas = (x[..., :-1] * y[..., 1:]).sum(-1, squeeze=True) - ( - x[..., 1:] * y[..., :-1] - ).sum(-1, squeeze=True) - - # TODO: replace this loop! - for i, (parts_x, parts_y) in enumerate(zip(x, y)): - for j, (nodes_x, nodes_y) in enumerate(zip(parts_x, parts_y)): - nodes_x = nodes_x.compressed() - nodes_y = nodes_y.compressed() - - if (nodes_x.size and nodes_x[0] != nodes_x[-1]) or ( - nodes_y.size and nodes_y[0] != nodes_y[-1] - ): - # First and last nodes of this polygon - # part are different => need to account - # for the "last" edge of the polygon that - # joins the first and last points. - all_areas[i, j] += x[-1] * y[0] - x[0] * y[-1] - - all_areas = all_areas.abs() * 0.5 - - elif x.Units.equivalent(_units_radians) and y.Units.equivalent( - _units_radians - ): - # -------------------------------------------------------- - # Spherical polygons defined by great circles - # - # The area is given by the sum of the interior angles - # minus (N-2)pi, where N is the number of sides: - # https://en.wikipedia.org/wiki/Spherical_trigonometry - # -------------------------------------------------------- - spherical = True - - if not great_circle: - raise ValueError( - "Must set great_circle=True to allow the derivation of " - "area weights from spherical polygons composed from " - "great circle segments." - ) - - if methods: - comp[(axis,)] = "area spherical polygon geometry" - return True - - x.Units = _units_radians - y.Units = _units_radians - x = x.array - y = y.array - - cols = np.ma.count(x, axis=1) - n_cells = cols.size - if (cols != np.ma.count(y, axis=1)).any(): - raise ValueError( - "Can't create area weights for " - f"{self.constructs.domain_axis_identity(axis)!r} axis: " - f"{aux_X!r} and {aux_Y!r} have inconsistent bounds " - "specifications" - ) - - col0 = cols[0] - if (cols == col0).all(): - # All cells have the same number of nodes, so we can - # use an integer and a slice in place of two integer - # arrays, which is much faster. - cols = col0 - rows = slice(None) - else: - rows = np.arange(x.shape[0]) - - if not ugrid: - # Ascertain whether or not the first and last boundary - # vertices are coincident in all cells. This can never - # be the case for UGRID cells. - cols_m1 = cols - 1 - duplicate = self._Data(x[:, 0]).isclose( - x[rows, cols_m1] - ) & self._Data(y[:, 0]).isclose(y[rows, cols_m1]) - - empty_col_x = np.ma.masked_all((n_cells, 1), dtype=x.dtype) - empty_col_y = np.ma.masked_all((n_cells, 1), dtype=y.dtype) - - if ugrid or not duplicate.any(): - # The first and last boundary vertices are different - # in all cells, i.e. No. nodes = No. edges. - # - # E.g. for triangle cells: - # [[4, 5, 6]] -> [[6, 4, 5, 6, 4]] - # [[4, 5, 6, --]] -> [[6, 4, 5, 6, 4, --]] - cols_p1 = cols + 1 - x = np.ma.concatenate((empty_col_x, x, empty_col_x), axis=1) - y = np.ma.concatenate((empty_col_y, y, empty_col_y), axis=1) - x[:, 0] = x[rows, cols] - y[:, 0] = y[rows, cols] - x[rows, cols_p1] = x[rows, 1] - y[rows, cols_p1] = y[rows, 1] - x = self._Data(x, units=_units_radians) - y = self._Data(y, units=_units_radians) - - # The number of edges of each polygon - N = cols - elif duplicate.all(): - # The first and last boundary vertices coincide in all - # cells, i.e. No. nodes = No. edges + 1. - # - # E.g. for triangle cells: - # [[4, 5, 6, 4]] -> [[6, 4, 5, 6, 4]] - # [[4, 5, 6, 4, --]] -> [[6, 4, 5, 6, 4, --]] - x = np.ma.concatenate((empty_col_x, x), axis=1) - y = np.ma.concatenate((empty_col_y, y), axis=1) - x[:, 0] = x[rows, cols_m1] - y[:, 0] = y[rows, cols_m1] - x = self._Data(x, units=_units_radians) - y = self._Data(y, units=_units_radians) - - # The number of edges of each polygon - N = cols - 1 - else: - raise ValueError( - "Can't calculate spherical polygon cell areas when " - "the first and last boundary vertices coincide in " - "some cells but not all" - ) - - # Find the interior angles of each polygon - interior_angles = self._weights_interior_angles( - x, y, interior_ring - ) - - # # Find the number of edges of each polygon - # N = interior_angles.count(-1, keepdims=False) - - # Find the polgon areas - all_areas = interior_angles.sum(-1, squeeze=True) - (N - 2) * np.pi - - # It's worth checking for negative areas, as this has been - # known to trap problems. TODOUGRID: is it? can we do this - # check on the weights later? Is it worth it even then? - if all_areas.min() < 0: - if ugrid: - raise ValueError( - "A UGRID spherical polygon cell has negative area" - ) - else: - raise ValueError( - "A spherical polygon geometry cell part has " - "negative area" - ) - else: - return False - - if ugrid: - areas = all_areas - else: - if interior_ring is not None: - # Change the sign of areas for geometry polygons that - # are interior rings - all_areas.where(interior_ring, -all_areas, inplace=True) - - # Sum the areas of each geometry polygon part to get the - # total area of each cell - areas = all_areas.sum(-1, squeeze=True) - - areas.override_units(_units_1, inplace=True) - - if measure and spherical: - # Multiply by radius squared, accounting for any Z - # coordinates, to get the actual area - if aux_Z is None: - # No Z coordinates - r = radius - else: - z = aux_Z.get_data(None, _fill_value=False) - if z is None: - # No Z coordinates - r = radius - else: - if not z.Units.equivalent(_units_metres): - raise ValueError( - "Z coordinates must have units equivalent to " - f"metres for area calculations. Got {z.Units!r}" - ) - - positive = aux_Z.get_property("positive", None) - if positive is None: - raise ValueError( - "Z coordinate 'positive' property is not defined" - ) - - if positive.lower() == "up": - r = radius + z - elif positive.lower() == "down": - r = radius - z - else: - raise ValueError( - "Bad value of Z coordinate 'positive' property: " - f"{positive!r}." - ) - - areas *= r**2 - - if return_areas: - return areas - - comp[(axis,)] = areas - - weights_axes.add(axis) - - return True - - def _weights_geometry_line( - self, - domain_axis, - comp, - weights_axes, - auto=False, - measure=False, - radius=None, - great_circle=False, - methods=False, - ): - """Creates line-length weights for line geometries. - - .. versionadded:: 3.2.0 - - :Parameters: - - measure: `bool` - If True then make sure that the weights represent true - cell sizes. - - methods: `bool`, optional - If True then add a description of the method used to - create the weights to the *comp* dictionary, as opposed to - the actual weights. - - """ - axis, aux_X, aux_Y, aux_Z = self._weights_yyy( - domain_axis, "line", methods=methods, auto=auto - ) - - if axis is None: - if auto: - return False - - if domain_axis is None: - raise ValueError("No line geometries") - - raise ValueError( - "No line geometries for " - f"{self.constructs.domain_axis_identity(domain_axis)!r} axis" - ) - - if axis in weights_axes: - if auto: - return False - - raise ValueError( - "Multiple weights specifications for " - f"{self.constructs.domain_axis_identity(axis)!r} axis" - ) - - x = aux_X.bounds.data - y = aux_Y.bounds.data - - if x.Units.equivalent(_units_metres) and y.Units.equivalent( - _units_metres - ): - # ---------------------------------------------------- - # Plane lines. - # - # Each line segment is the simple cartesian distance - # between two adjacent nodes. - # ---------------------------------------------------- - if methods: - comp[(axis,)] = "linear plane line geometry" - return True - - y.Units = x.Units - - delta_x = x.diff(axis=-1) - delta_y = y.diff(axis=-1) - - all_lengths = (delta_x**2 + delta_y**2) ** 0.5 - all_lengths = all_lengths.sum(-1, squeeze=True) - - elif x.Units.equivalent(_units_radians) and y.Units.equivalent( - _units_radians - ): - # ---------------------------------------------------- - # Spherical lines. - # - # Each line segment is a great circle arc between two - # adjacent nodes. - # - # The length of the great circle arc is the the central - # angle multiplied by the radius. The central angle is - # calculated with a special case of the Vincenty formula: - # https://en.wikipedia.org/wiki/Great-circle_distance - # ---------------------------------------------------- - if not great_circle: - raise ValueError( - "Must set great_circle=True to allow the derivation " - "of line-length weights from great circle segments." - ) - - if methods: - comp[(axis,)] = "linear spherical line geometry" - return True - - x.Units = _units_radians - y.Units = _units_radians - - central_angles = self._weights_central_angles(x, y) - - # It's worth checking for negative areas, as this has been - # known to trap problems. TODOUGRID: is it? can we do this - # check on the weights later? Is it worth it even then? - if central_angles.min() < 0: - raise ValueError( - "A spherical line geometry segment has " - f"negative length: {central_angles.min() * radius!r}" - ) - - all_lengths = central_angles.sum(-1, squeeze=True) - - if measure: - all_lengths *= radius - else: - return False - - # Sum the lengths of each part to get the total length of - # each cell - lengths = all_lengths.sum(-1, squeeze=True) - - comp[(axis,)] = lengths - - weights_axes.add(axis) - - return True - - def _weights_geometry_volume( - self, - comp, - weights_axes, - auto=False, - measure=False, - radius=None, - great_circle=False, - methods=False, - ): - """Creates volume weights for polygon geometry cells. - - .. versionadded:: 3.2.0 - - :Parameters: - - measure: `bool` - If True then make sure that the weights represent true - cell sizes. - - methods: `bool`, optional - If True then add a description of the method used to - create the weights to the *comp* dictionary, as opposed to - the actual weights. - - """ - axis, aux_X, aux_Y, aux_Z = self._weights_yyy( - "polygon", methods=methods, auto=auto - ) - - if axis is None and auto: - return False - - if axis in weights_axes: - if auto: - return False - - raise ValueError( - "Multiple weights specifications for " - f"{self.constructs.domain_axis_identity(axis)!r} axis" - ) - - x = aux_X.bounds.data - y = aux_Y.bounds.data - z = aux_Z.bounds.data - - if not z.Units.equivalent(_units_metres): - if auto: - return False - - raise ValueError( - "Z coordinate bounds must have units equivalent to " - f"metres for volume calculations. Got {z.Units!r}." - ) - - if not methods: - # Initialise cell volumes as the cell areas - volumes = self._weights_geometry_area( - comp, - weights_axes, - auto=auto, - measure=measure, - radius=radius, - great_circle=great_circle, - methods=False, - return_areas=True, - ) - - if measure: - delta_z = abs(z[..., 1] - z[..., 0]) - delta_z.squeeze(axis=-1, inplace=True) - - if x.Units.equivalent(_units_metres) and y.Units.equivalent( - _units_metres - ): - # ---------------------------------------------------- - # Plane polygons defined by straight lines. - # - # Do this in preference to weights based on spherical - # polygons, which require the great circle assumption. - # ---------------------------------------------------- - if methods: - comp[(axis,)] = "volume plane polygon geometry" - return True - - if measure: - volumes *= delta_z - - elif x.Units.equivalent(_units_radians) and y.Units.equivalent( - _units_radians - ): - # ---------------------------------------------------- - # Spherical polygons defined by great circles - # - # The area of such a spherical polygon is given by the - # sum of the interior angles minus (N-2)pi, where N is - # the number of sides (Todhunter): - # - # https://en.wikipedia.org/wiki/Spherical_trigonometry#Area_and_spherical_excess - # - # The interior angle of a side is calculated with a - # special case of the Vincenty formula: - # https://en.wikipedia.org/wiki/Great-circle_distance - # ---------------------------------------------------- - if not great_circle: - raise ValueError( - "Must set great_circle=True to allow the derivation " - "of volume weights from spherical polygons composed " - "from great circle segments." - ) - - if methods: - comp[(axis,)] = "volume spherical polygon geometry" - return True - - if measure: - r = radius - - # actual_volume = - # [actual_area/(4*pi*r**2)] - # * (4/3)*pi*[(r+delta_z)**3 - r**3)] - volumes *= ( - delta_z**3 / (3 * r**2) + delta_z**2 / r + delta_z - ) - else: - raise ValueError( - "X and Y coordinate bounds must both have units " - "equivalent to either metres, for plane polygon, or radians, " - "for spherical polygon, volume calculations. Got " - f"{x.Units!r} and {y.Units!r}." - ) - - comp[(axis,)] = volumes - - weights_axes.add(axis) - - return True - - def _weights_interior_angles(self, lon, lat, interior_ring=None): - """Find the interior angles at spherical polygon vertices. - - The interior angle at a vertex is that formed by two adjacent - sides. Each interior angle is therefore a function of three - vertices - the vertex at which the interior angle is required - and the two adjacent vertices either side of it. - - **Method** - - Bevis, M., Cambareri, G. Computing the area of a spherical - polygon of arbitrary shape. Math Geol 19, 335–346 - (1987). - - http://bemlar.ism.ac.jp/zhuang/Refs/Refs/bevis1987mathgeol.pdf - - Abstract: An algorithm for determining the area of a spherical - polygon of arbitrary shape is presented. The kernel of the - problem is to compute the interior angle at each vertex of - the spherical polygon; a well-known relationship between - the area of a spherical polygon and the sum of its - interior angles then may be exploited. The algorithm has - been implemented as a FORTRAN subroutine and a listing is - provided. - - .. versionadded:: TODOUGRIDVER - - .. seealso:: `_weights_central_angles` - - :Parameters: - - lon: `Data` - Longitudes. Must have units of radians, which is not - checked. - - lat: `Data` - Latitudes. Must have units of radians, which is not - checked. - - interior_ring: `Data`, optional - The interior ring flags for geometry cells. - - :Returns: - - `Data` - The interior angles of each spherical polygon, in - units of radians. - - """ - # P denotes a vertex at which the interior angle is required, - # A denotes the adjacent point clockwise from P, and B denotes - # the adjacent point anticlockwise from P. - - lon_A = lon[..., :-2] - lon_P = lon[..., 1:-1] - lon_B = lon[..., 2:] - - lon_A_minus_lon_P = lon_A - lon_P - lon_B_minus_lon_P = lon_B - lon_P - - cos_lat = lat.cos() - cos_lat_A = cos_lat[..., :-2] - cos_lat_P = cos_lat[..., 1:-1] - cos_lat_B = cos_lat[..., 2:] - - sin_lat = lat.sin() - sin_lat_A = sin_lat[..., :-2] - sin_lat_P = sin_lat[..., 1:-1] - sin_lat_B = sin_lat[..., 2:] - - y = lon_A_minus_lon_P.sin() * cos_lat_A - x = ( - sin_lat_A * cos_lat_P - - cos_lat_A * sin_lat_P * lon_A_minus_lon_P.cos() - ) - lat_A_primed = self._Data.arctan2(y, x) - - y = lon_B_minus_lon_P.sin() * cos_lat_B - x = ( - sin_lat_B * cos_lat_P - - cos_lat_B * sin_lat_P * lon_B_minus_lon_P.cos() - ) - lat_B_primed = self._Data.arctan2(y, x) - - # The CF vertices here are, in general, given in anticlockwise - # order, so we do alpha_P = lat_B_primed - lat_A_primed, - # rather than the "alpha_P = lat_A_primed - lat_B_primed" - # given in Bevis and Cambareri, which assumes clockwise order. - alpha_P = lat_B_primed - lat_A_primed - - if interior_ring is not None: - # However, interior rings *are* given in clockwise order - # in CF, so we need to negate alpha_P in these cases. - alpha_P.where(interior_ring, -alpha_P, inplace=True) - - alpha_P = alpha_P.where(lt(0), alpha_P + 2 * np.pi, alpha_P) - return alpha_P - - def _weights_central_angles(self, lon, lat): - r"""Find the central angles for spherical great circle line segments. - - The central angle of two points on the sphere is the angle - subtended from the centre of the sphere by the two points on - its surface points. It is calculated with a special case of - the Vincenty formula that is accurate for all spherical - distances: - - \Delta \sigma =\arctan {\frac {\sqrt {\left(\cos \phi - _{2}\sin(\Delta \lambda )\right)^{2}+\left(\cos \phi _{1}\sin - \phi _{2}-\sin \phi _{1}\cos \phi _{2}\cos(\Delta \lambda - )\right)^{2}}}{\sin \phi _{1}\sin \phi _{2}+\cos \phi _{1}\cos - \phi _{2}\cos(\Delta \lambda )}} - - The quadrant for \Delta \sigma should be governed by the signs - of the numerator and denominator of the right hand side using - the atan2 function. - - (https://en.wikipedia.org/wiki/Great-circle_distance) - - :Parameters: - - lon: `Data` - Longitudes. Must have units of radians, which is not - checked. - - lat: `Data` - Latitudes. Must have units of radians, which is not - checked. - - :Returns: - - `Data` - The central angles in units of radians. - - """ - # A and B denote adjacent vertices that define each line - # segment - - delta_lon = lon.diff(axis=-1) - - cos_lat = lat.cos() - sin_lat = lat.sin() - - cos_lat_A = cos_lat[..., :-1] - cos_lat_B = cos_lat[..., 1:] - - sin_lat_A = sin_lat[..., :-1] - sin_lat_B = sin_lat[..., 1:] - - cos_delta_lon = delta_lon.cos() - sin_delta_lon = delta_lon.sin() - - y = ( - (cos_lat_B * sin_delta_lon) ** 2 - + (cos_lat_A * sin_lat_B - sin_lat_A * cos_lat_B * cos_delta_lon) - ** 2 - ) ** 0.5 - x = sin_lat_A * sin_lat_B + cos_lat_A * cos_lat_B * cos_delta_lon - - delta_sigma = self._Data.arctan2(y, x) - delta_sigma.override_units(_units_1, inplace=True) - return delta_sigma - - def _weights_linear( - self, - axis, - comp, - weights_axes, - auto=False, - measure=False, - methods=False, - ): - """1-d linear weights from dimension coordinate constructs. - - :Parameters: - - axis: `str` - - comp: `dict` - - weights_axes: `set` - - auto: `bool` - - measure: `bool` - If True then make sure that the weights represent true - cell sizes. - - methods: `bool`, optional - If True then add a description of the method used to - create the weights to the *comp* dictionary, as opposed to - the actual weights. - - :Returns: - - `bool` - - """ - da_key = self.domain_axis(axis, key=True, default=None) - if da_key is None: - if auto: - return False - - raise ValueError( - "Can't create weights: Can't find domain axis " - f"matching {axis!r}" - ) - - dim = self.dimension_coordinate(filter_by_axis=(da_key,), default=None) - if dim is None: - if auto: - return False - - raise ValueError( - f"Can't create linear weights for {axis!r} axis: Can't find " - "dimension coordinate construct." - ) - - if not measure and dim.size == 1: - return False - - if da_key in weights_axes: - if auto: - return False - - raise ValueError( - f"Can't create linear weights for {axis!r} axis: Multiple " - "axis specifications" - ) - - if not dim.has_bounds(): - # Dimension coordinate has no bounds - if auto: - return False - - raise ValueError( - f"Can't create linear weights for {axis!r} axis: No bounds" - ) - else: - # Bounds exist - if methods: - comp[ - (da_key,) - ] = "linear " + self.constructs.domain_axis_identity(da_key) - else: - comp[(da_key,)] = dim.cellsize - - weights_axes.add(da_key) - - return True - - def _weights_measure( - self, measure, comp, weights_axes, methods=False, auto=False - ): - """Cell measure weights. - - :Parameters: - - methods: `bool`, optional - If True then add a description of the method used to - create the weights to the *comp* dictionary, as opposed to - the actual weights. - - :Returns: - - `bool` - - """ - m = self.cell_measures(filter_by_measure=(measure,), todict=True) - len_m = len(m) - - if not len_m: - if measure == "area": - return False - - if auto: - return - - raise ValueError( - f"Can't find weights: No {measure!r} cell measure" - ) - - elif len_m > 1: - if auto: - return False - - raise ValueError( - f"Can't find weights: Multiple {measure!r} cell measures" - ) - - key, clm = m.popitem() - - clm_axes0 = self.get_data_axes(key) - - clm_axes = tuple( - [axis for axis, n in zip(clm_axes0, clm.data.shape) if n > 1] - ) - - for axis in clm_axes: - if axis in weights_axes: - if auto: - return False - - raise ValueError( - "Multiple weights specifications for " - f"{self.constructs.domain_axis_identity(axis)!r} axis" - ) - - clm = clm.get_data(_fill_value=False).copy() - if clm_axes != clm_axes0: - iaxes = [clm_axes0.index(axis) for axis in clm_axes] - clm.squeeze(iaxes, inplace=True) - - if methods: - comp[tuple(clm_axes)] = measure + " cell measure" - else: - comp[tuple(clm_axes)] = clm - - weights_axes.update(clm_axes) - - return True - - def _weights_scale(self, w, scale): - """Scale the weights so that they are <= scale. - - :Parameters: - - w: `Data` - The weights to be scaled. - - scale: number or `None` - The maximum value of the scaled weights. If `None` then no - scaling is applied. - - :Returns: - - `Data` - The scaled weights. - - """ - if scale is None: - return w - - if scale <= 0: - raise ValueError( - "Can't set 'scale' parameter to a negative number. " - f"Got {scale!r}" - ) - - w = w / w.max() - if scale != 1: - w = w * scale - - return w - - def _weights_yyy( - self, domain_axis, geometry_type, methods=False, auto=False - ): - """Checks whether weights can be created for given coordinates. - - .. versionadded:: 3.2.0 - - :Parameters: - - domain_axis : `str` or `None` - - geometry_type: `str` - Either ``'polygon'`` or ``'line'``. - - auto: `bool` - - :Returns: - - `tuple` - - """ - aux_X = None - aux_Y = None - aux_Z = None - x_axis = None - y_axis = None - z_axis = None - - auxiliary_coordinates_1d = self.auxiliary_coordinates( - filter_by_naxes=(1,), todict=True - ) - - for key, aux in auxiliary_coordinates_1d.items(): - aux_axis = self.get_data_axes(key)[0] - dt = self.domain_topology(default=None, filter_by_axis=(aux_axis,)) - if dt is None and aux.get_geometry(None) != geometry_type: - continue - - # Still here? Then this auxiliary coordinate has either - # UGRID or geometry cells. - if aux.X: - aux_X = aux.copy() - x_axis = aux_axis # self.get_data_axes(key)[0] - if domain_axis is not None and x_axis != domain_axis: - aux_X = None - continue - elif aux.Y: - aux_Y = aux.copy() - y_axis = aux_axis # self.get_data_axes(key)[0] - if domain_axis is not None and y_axis != domain_axis: - aux_Y = None - continue - elif aux.Z: - aux_Z = aux.copy() - z_axis = aux_axis # self.get_data_axes(key)[0] - if domain_axis is not None and z_axis != domain_axis: - aux_Z = None - continue - - if aux_X is None or aux_Y is None: - if auto: - return (None,) * 4 - - raise ValueError( - "Can't create weights: Need both X and Y nodes to " - f"calculate {geometry_type} cell weights" - ) - - if x_axis != y_axis: - if auto: - return (None,) * 4 - - raise ValueError( - "Can't create weights: X and Y cells span different " - "domain axes" - ) - - axis = x_axis - - if aux_X.get_bounds(None) is None or aux_Y.get_bounds(None) is None: - # Not both X and Y coordinates have bounds - if auto: - return (None,) * 4 - - raise ValueError( - "Can't create weights: " - "Not both X and Y coordinates have bounds" - ) - - if aux_X.bounds.shape != aux_Y.bounds.shape: - if auto: - return (None,) * 4 - - raise ValueError( - "Can't create weights: UGRID/geometry X and Y coordinate " - "bounds must have the same shape. " - f"Got {aux_X.bounds.shape} and {aux_Y.bounds.shape}" - ) - - if aux_Z is None: - for key, aux in auxiliary_coordinates_1d.items(): - if aux.Z: - aux_Z = aux.copy() - z_axis = self.get_data_axes(key)[0] - - # Check Z coordinates - if aux_Z is not None: - if z_axis != x_axis: - if auto: - return (None,) * 4 - - raise ValueError( - "Z coordinates span different domain axis to X and Y " - "geometry coordinates" - ) - - return axis, aux_X, aux_Y, aux_Z + # def _weights_area_XY( + # self, + # comp, + # weights_axes, + # auto=False, + # measure=False, + # radius=None, + # methods=False, + # ): + # """Calculate area weights from X and Y dimension coordinate + # constructs. + # + # :Parameters: + # + # measure: `bool` + # If True then make sure that the weights represent true + # cell sizes. + # + # methods: `bool`, optional + # If True then add a description of the method used to + # create the weights to the *comp* dictionary, as opposed to + # the actual weights. + # + # :Returns: + # + # `bool` or `None` + # + # """ + # xkey, xcoord = self.dimension_coordinate( + # "X", item=True, default=(None, None) + # ) + # ykey, ycoord = self.dimension_coordinate( + # "Y", item=True, default=(None, None) + # ) + # + # if xcoord is None or ycoord is None: + # if auto: + # return + # + # raise ValueError( + # "No unique 'X' and 'Y' dimension coordinate constructs for " + # "calculating area weights" + # ) + # + # if xcoord.Units.equivalent( + # Units("radians") + # ) and ycoord.Units.equivalent(Units("radians")): + # pass + # elif xcoord.Units.equivalent( + # Units("metres") + # ) and ycoord.Units.equivalent(Units("metres")): + # radius = None + # else: + # if auto: + # return + # + # raise ValueError( + # "Insufficient coordinate constructs for calculating " + # "area weights" + # ) + # + # xaxis = self.get_data_axes(xkey)[0] + # yaxis = self.get_data_axes(ykey)[0] + # + # for axis in (xaxis, yaxis): + # if axis in weights_axes: + # if auto: + # return + # + # raise ValueError( + # "Multiple weights specifications for " + # f"{self.constructs.domain_axis_identity(axis)!r} axis" + # ) + # + # if measure and radius is not None: + # radius = self.radius(default=radius) + # + # if measure or xcoord.size > 1: + # if not xcoord.has_bounds(): + # if auto: + # return + # + # raise ValueError( + # "Can't create area weights: No bounds for " + # f"{xcoord.identity()!r} axis" + # ) + # + # if methods: + # comp[(xaxis,)] = "linear " + xcoord.identity() + # else: + # cells = xcoord.cellsize + # if xcoord.Units.equivalent(Units("radians")): + # cells.Units = _units_radians + # if measure: + # cells *= radius + # cells.override_units(radius.Units, inplace=True) + # else: + # cells.Units = Units("metres") + # + # comp[(xaxis,)] = cells + # + # weights_axes.add(xaxis) + # + # if measure or ycoord.size > 1: + # if not ycoord.has_bounds(): + # if auto: + # return + # + # raise ValueError( + # "Can't create area weights: No bounds for " + # f"{ycoord.identity()!r} axis" + # ) + # + # if ycoord.Units.equivalent(Units("radians")): + # ycoord = ycoord.clip(-90, 90, units=Units("degrees")) + # ycoord.sin(inplace=True) + # + # if methods: + # comp[(yaxis,)] = "linear sine " + ycoord.identity() + # else: + # cells = ycoord.cellsize + # if measure: + # cells *= radius + # + # comp[(yaxis,)] = cells + # else: + # if methods: + # comp[(yaxis,)] = "linear " + ycoord.identity() + # else: + # cells = ycoord.cellsize + # comp[(yaxis,)] = cells + # + # weights_axes.add(yaxis) + # + # return True + # + # def _weights_data( + # self, + # w, + # comp, + # weights_axes, + # axes=None, + # data=False, + # components=False, + # methods=False, + # ): + # """Creates weights for the field construct's data array. + # + # :Parameters: + # + # methods: `bool`, optional + # If True then add a description of the method used to + # create the weights to the *comp* dictionary, as + # opposed to the actual weights. + # + # """ + # # -------------------------------------------------------- + # # Data weights + # # -------------------------------------------------------- + # field_data_axes = self.get_data_axes() + # + # if axes is not None: + # if isinstance(axes, (str, int)): + # axes = (axes,) + # + # if len(axes) != w.ndim: + # raise ValueError( + # "'axes' parameter must provide an axis identifier " + # f"for each weights data dimension. Got {axes!r} for " + # f"{w.ndim} dimension(s)." + # ) + # + # iaxes = [ + # field_data_axes.index(self.domain_axis(axis, key=True)) + # for axis in axes + # ] + # + # if data: + # for i in range(self.ndim): + # if i not in iaxes: + # w = w.insert_dimension(position=i) + # iaxes.insert(i, i) + # + # w = w.transpose(iaxes) + # + # if w.ndim > 0: + # while w.shape[0] == 1: + # w = w.squeeze(0) + # + # else: + # iaxes = list(range(self.ndim - w.ndim, self.ndim)) + # + # if not (components or methods): + # if not self._is_broadcastable(w.shape): + # raise ValueError( + # f"The 'Data' weights (shape {w.shape}) are not " + # "broadcastable to the field construct's data " + # f"(shape {self.shape})." + # ) + # + # axes0 = field_data_axes[self.ndim - w.ndim :] + # else: + # axes0 = [field_data_axes[i] for i in iaxes] + # + # for axis0 in axes0: + # if axis0 in weights_axes: + # raise ValueError( + # "Multiple weights specified for " + # f"{self.constructs.domain_axis_identity(axis0)!r} axis" + # ) + # + # if methods: + # comp[tuple(axes0)] = "custom data" + # else: + # comp[tuple(axes0)] = w + # + # weights_axes.update(axes0) + # + # return True + # + # def _weights_field(self, fields, comp, weights_axes, methods=False): + # """Creates a weights field.""" + # s = self.analyse_items() + # + # domain_axes = self.domain_axes(todict=True) + # for w in fields: + # t = w.analyse_items() + # # TODO CHECK this with org + # + # if t["undefined_axes"]: + # w_domain_axes_1 = w.domain_axes( + # filter_by_size=(1,), todict=True + # ) + # if set(w_domain_axes_1).intersection(t["undefined_axes"]): + # raise ValueError( + # f"Weights field {w} has at least one undefined " + # f"domain axes: {w_domain_axes_1}." + # ) + # + # w = w.squeeze() + # + # w_domain_axes = w.domain_axes(todict=True) + # + # axis1_to_axis0 = {} + # + # coordinate_references = self.coordinate_references(todict=True) + # w_coordinate_references = w.coordinate_references(todict=True) + # + # for axis1 in w.get_data_axes(): + # identity = t["axis_to_id"].get(axis1, None) + # if identity is None: + # raise ValueError( + # "Weights field has unmatched, size > 1 " + # f"{w.constructs.domain_axis_identity(axis1)!r} axis" + # ) + # + # axis0 = s["id_to_axis"].get(identity, None) + # if axis0 is None: + # raise ValueError( + # f"Weights field has unmatched, size > 1 {identity!r} " + # "axis" + # ) + # + # w_axis_size = w_domain_axes[axis1].get_size() + # self_axis_size = domain_axes[axis0].get_size() + # + # if w_axis_size != self_axis_size: + # raise ValueError( + # f"Weights field has incorrectly sized {identity!r} " + # f"axis ({w_axis_size} != {self_axis_size})" + # ) + # + # axis1_to_axis0[axis1] = axis0 + # + # # Check that the defining coordinate data arrays are + # # compatible + # key0 = s["axis_to_coord"][axis0] + # key1 = t["axis_to_coord"][axis1] + # + # if not self._equivalent_construct_data( + # w, key0=key0, key1=key1, s=s, t=t + # ): + # raise ValueError( + # f"Weights field has incompatible {identity!r} " + # "coordinates" + # ) + # + # # Still here? Then the defining coordinates have + # # equivalent data arrays + # + # # If the defining coordinates are attached to + # # coordinate references then check that those + # # coordinate references are equivalent + # refs0 = [ + # key + # for key, ref in coordinate_references.items() + # if key0 in ref.coordinates() + # ] + # refs1 = [ + # key + # for key, ref in w_coordinate_references.items() + # if key1 in ref.coordinates() + # ] + # + # nrefs = len(refs0) + # if nrefs > 1 or nrefs != len(refs1): + # # The defining coordinate are associated with + # # different numbers of coordinate references + # equivalent_refs = False + # elif not nrefs: + # # Neither defining coordinate is associated with a + # # coordinate reference + # equivalent_refs = True + # else: + # # Each defining coordinate is associated with + # # exactly one coordinate reference + # equivalent_refs = self._equivalent_coordinate_references( + # w, key0=refs0[0], key1=refs1[0], s=s, t=t + # ) + # + # if not equivalent_refs: + # raise ValueError( + # "Input weights field has an incompatible " + # "coordinate reference" + # ) + # + # axes0 = tuple( + # [axis1_to_axis0[axis1] for axis1 in w.get_data_axes()] + # ) + # + # for axis0 in axes0: + # if axis0 in weights_axes: + # raise ValueError( + # "Multiple weights specified for " + # f"{self.constructs.domain_axis_identity(axis0)!r} " + # "axis" + # ) + # + # comp[tuple(axes0)] = w.data + # + # weights_axes.update(axes0) + # + # return True + # + # def _weights_field_scalar(self, methods=False): + # """Return a scalar field of weights with long_name ``'weight'``. + # + # :Returns: + # + # `Field` + # + # """ + # data = Data(1.0, "1") + # + # f = type(self)() + # f.set_data(data, copy=False) + # f.long_name = "weight" + # f.comment = f"Weights for {self!r}" + # + # return f + # + # def _weights_geometry_area( + # self, + # domain_axis, + # comp, + # weights_axes, + # auto=False, + # measure=False, + # radius=None, + # great_circle=False, + # return_areas=False, + # methods=False, + # ): + # """Creates area weights for polygon geometry cells. + # + # .. versionadded:: 3.2.0 + # + # :Parameters: + # + # domain_axis : `str` or `None` + # + # measure: `bool` + # If True then make sure that the weights represent true + # cell sizes. + # + # methods: `bool`, optional + # If True then add a description of the method used to + # create the weights to the *comp* dictionary, as opposed to + # the actual weights. + # + # :Returns: + # + # `bool` or `Data` + # + # """ + # axis, aux_X, aux_Y, aux_Z = self._weights_yyy( + # domain_axis, "polygon", methods=methods, auto=auto + # ) + # + # if axis is None: + # if auto: + # return False + # + # if domain_axis is None: + # raise ValueError("No polygon geometries") + # + # raise ValueError( + # "No polygon geometries for " + # f"{self.constructs.domain_axis_identity(domain_axis)!r} axis" + # ) + # + # if axis in weights_axes: + # if auto: + # return False + # + # raise ValueError( + # "Multiple weights specifications for " + # f"{self.constructs.domain_axis_identity(axis)!r} axis" + # ) + # + # ugrid = self.domain_topologies(filter_by_axis=(axis,), todict=True) + # + # # Check for interior rings + # if ugrid: + # interior_ring = None + # else: + # interior_ring_X = aux_X.get_interior_ring(None) + # interior_ring_Y = aux_Y.get_interior_ring(None) + # if interior_ring_X is None and interior_ring_Y is None: + # interior_ring = None + # elif interior_ring_X is None: + # raise ValueError( + # "Can't find weights: X coordinates have missing " + # "interior ring variable" + # ) + # elif interior_ring_Y is None: + # raise ValueError( + # "Can't find weights: Y coordinates have missing " + # "interior ring variable" + # ) + # elif not interior_ring_X.data.equals(interior_ring_Y.data): + # raise ValueError( + # "Can't find weights: Interior ring variables for X and Y " + # "coordinates have different data values" + # ) + # else: + # interior_ring = interior_ring_X.data + # if interior_ring.shape != aux_X.bounds.shape[:-1]: + # raise ValueError( + # "Can't find weights: Interior ring variables have " + # f"incorrect shape. Got {interior_ring.shape}, " + # f"expected {aux_X.bounds.shape[:-1]}" + # # TODOUGRID: is this not tested in set_interior_ring? + # ) + # + # x = aux_X.bounds.data + # y = aux_Y.bounds.data + # spherical = False + # + # if x.Units.equivalent(_units_metres) and y.Units.equivalent( + # _units_metres + # ): + # # -------------------------------------------------------- + # # Plane polygons defined by straight lines. + # # + # # The area is given by the shoelace formula: + # # https://en.wikipedia.org/wiki/Shoelace_formula + # # + # # Do this in preference to weights based on spherical + # # polygons, which require the great circle assumption. + # # -------------------------------------------------------- + # if methods: + # comp[(axis,)] = "area plane polygon geometry" + # return True + # + # y.Units = x.Units + # + # all_areas = (x[..., :-1] * y[..., 1:]).sum(-1, squeeze=True) - ( + # x[..., 1:] * y[..., :-1] + # ).sum(-1, squeeze=True) + # + # # TODO: replace this loop! + # for i, (parts_x, parts_y) in enumerate(zip(x, y)): + # for j, (nodes_x, nodes_y) in enumerate(zip(parts_x, parts_y)): + # nodes_x = nodes_x.compressed() + # nodes_y = nodes_y.compressed() + # + # if (nodes_x.size and nodes_x[0] != nodes_x[-1]) or ( + # nodes_y.size and nodes_y[0] != nodes_y[-1] + # ): + # # First and last nodes of this polygon + # # part are different => need to account + # # for the "last" edge of the polygon that + # # joins the first and last points. + # all_areas[i, j] += x[-1] * y[0] - x[0] * y[-1] + # + # all_areas = all_areas.abs() * 0.5 + # + # elif x.Units.equivalent(_units_radians) and y.Units.equivalent( + # _units_radians + # ): + # # -------------------------------------------------------- + # # Spherical polygons defined by great circles + # # + # # The area is given by the sum of the interior angles + # # minus (N-2)pi, where N is the number of sides: + # # https://en.wikipedia.org/wiki/Spherical_trigonometry + # # -------------------------------------------------------- + # spherical = True + # + # if not great_circle: + # raise ValueError( + # "Must set great_circle=True to allow the derivation of " + # "area weights from spherical polygons composed from " + # "great circle segments." + # ) + # + # if methods: + # comp[(axis,)] = "area spherical polygon geometry" + # return True + # + # x.Units = _units_radians + # y.Units = _units_radians + # x = x.array + # y = y.array + # + # cols = np.ma.count(x, axis=1) + # n_cells = cols.size + # if (cols != np.ma.count(y, axis=1)).any(): + # raise ValueError( + # "Can't create area weights for " + # f"{self.constructs.domain_axis_identity(axis)!r} axis: " + # f"{aux_X!r} and {aux_Y!r} have inconsistent bounds " + # "specifications" + # ) + # + # col0 = cols[0] + # if (cols == col0).all(): + # # All cells have the same number of nodes, so we can + # # use an integer and a slice in place of two integer + # # arrays, which is much faster. + # cols = col0 + # rows = slice(None) + # else: + # rows = np.arange(x.shape[0]) + # + # if not ugrid: + # # Ascertain whether or not the first and last boundary + # # vertices are coincident in all cells. Note that this + # # can never be the case for UGRID cells. + # cols_m1 = cols - 1 + # duplicate = self._Data(x[:, 0]).isclose( + # x[rows, cols_m1] + # ) & self._Data(y[:, 0]).isclose(y[rows, cols_m1]) + # + # empty_col_x = np.ma.masked_all((n_cells, 1), dtype=x.dtype) + # empty_col_y = np.ma.masked_all((n_cells, 1), dtype=y.dtype) + # + # # Pad out the boundary vertiex array for each cell with + # # first and last bounds neighbours + # if ugrid or not duplicate.any(): + # # The first and last boundary vertices are different + # # in all cells, i.e. No. nodes = No. edges. + # # + # # E.g. for triangle cells: + # # [[4, 5, 6]] -> [[6, 4, 5, 6, 4]] + # # [[4, 5, 6, --]] -> [[6, 4, 5, 6, 4, --]] + # cols_p1 = cols + 1 + # x = np.ma.concatenate((empty_col_x, x, empty_col_x), axis=1) + # y = np.ma.concatenate((empty_col_y, y, empty_col_y), axis=1) + # x[:, 0] = x[rows, cols] + # y[:, 0] = y[rows, cols] + # x[rows, cols_p1] = x[rows, 1] + # y[rows, cols_p1] = y[rows, 1] + # x = self._Data(x, units=_units_radians) + # y = self._Data(y, units=_units_radians) + # + # # The number of edges of each polygon + # N = cols + # elif duplicate.all(): + # # The first and last boundary vertices coincide in all + # # cells, i.e. No. nodes = No. edges + 1. + # # + # # E.g. for triangle cells: + # # [[4, 5, 6, 4]] -> [[6, 4, 5, 6, 4]] + # # [[4, 5, 6, 4, --]] -> [[6, 4, 5, 6, 4, --]] + # x = np.ma.concatenate((empty_col_x, x), axis=1) + # y = np.ma.concatenate((empty_col_y, y), axis=1) + # x[:, 0] = x[rows, cols_m1] + # y[:, 0] = y[rows, cols_m1] + # x = self._Data(x, units=_units_radians) + # y = self._Data(y, units=_units_radians) + # + # # The number of edges of each polygon + # N = cols - 1 + # else: + # raise ValueError( + # "Can't calculate spherical polygon cell areas when " + # "the first and last boundary vertices coincide in " + # "some cells but not all" + # ) + # + # # Find the interior angles of each polygon + # interior_angles = self._weights_interior_angles( + # x, y, interior_ring + # ) + # + # # Find the polgon areas + # all_areas = interior_angles.sum(-1, squeeze=True) - (N - 2) * np.pi + # + # # It's worth checking for negative areas, as this has been + # # known to trap problems. TODOUGRID: is it? can we do this + # # check on the weights later? Is it worth it even then? + # if all_areas.min() < 0: + # if ugrid: + # raise ValueError( + # "A UGRID spherical polygon cell has negative area" + # ) + # else: + # raise ValueError( + # "A spherical polygon geometry cell part has " + # "negative area" + # ) + # else: + # return False + # + # if ugrid: + # areas = all_areas + # else: + # if interior_ring is not None: + # # Change the sign of areas for geometry polygons that + # # are interior rings + # all_areas.where(interior_ring, -all_areas, inplace=True) + # + # # Sum the areas of each geometry polygon part to get the + # # total area of each cell + # areas = all_areas.sum(-1, squeeze=True) + # + # areas.override_units(_units_1, inplace=True) + # + # if measure and spherical: + # # Multiply by radius squared, accounting for any Z + # # coordinates, to get the actual area + # if aux_Z is None: + # # No Z coordinates + # r = radius + # else: + # z = aux_Z.get_data(None, _fill_value=False) + # if z is None: + # # No Z coordinates + # r = radius + # else: + # if not z.Units.equivalent(_units_metres): + # raise ValueError( + # "Z coordinates must have units equivalent to " + # f"metres for area calculations. Got {z.Units!r}" + # ) + # + # positive = aux_Z.get_property("positive", None) + # if positive is None: + # raise ValueError( + # "Z coordinate 'positive' property is not defined" + # ) + # + # if positive.lower() == "up": + # r = radius + z + # elif positive.lower() == "down": + # r = radius - z + # else: + # raise ValueError( + # "Bad value of Z coordinate 'positive' property: " + # f"{positive!r}." + # ) + # + # areas *= r**2 + # + # if return_areas: + # return areas + # + # comp[(axis,)] = areas + # + # weights_axes.add(axis) + # + # return True + # + # def _weights_geometry_line( + # self, + # domain_axis, + # comp, + # weights_axes, + # auto=False, + # measure=False, + # radius=None, + # great_circle=False, + # methods=False, + # ): + # """Creates line-length weights for line geometries. + # + # .. versionadded:: 3.2.0 + # + # :Parameters: + # + # measure: `bool` + # If True then make sure that the weights represent true + # cell sizes. + # + # methods: `bool`, optional + # If True then add a description of the method used to + # create the weights to the *comp* dictionary, as opposed to + # the actual weights. + # + # """ + # axis, aux_X, aux_Y, aux_Z = self._weights_yyy( + # domain_axis, "line", methods=methods, auto=auto + # ) + # + # if axis is None: + # if auto: + # return False + # + # if domain_axis is None: + # raise ValueError("No line geometries") + # + # raise ValueError( + # "No line geometries for " + # f"{self.constructs.domain_axis_identity(domain_axis)!r} axis" + # ) + # + # if axis in weights_axes: + # if auto: + # return False + # + # raise ValueError( + # "Multiple weights specifications for " + # f"{self.constructs.domain_axis_identity(axis)!r} axis" + # ) + # + # x = aux_X.bounds.data + # y = aux_Y.bounds.data + # + # if x.Units.equivalent(_units_metres) and y.Units.equivalent( + # _units_metres + # ): + # # ---------------------------------------------------- + # # Plane lines. + # # + # # Each line segment is the simple cartesian distance + # # between two adjacent nodes. + # # ---------------------------------------------------- + # if methods: + # comp[(axis,)] = "linear plane line geometry" + # return True + # + # y.Units = x.Units + # + # delta_x = x.diff(axis=-1) + # delta_y = y.diff(axis=-1) + # + # all_lengths = (delta_x**2 + delta_y**2) ** 0.5 + # all_lengths = all_lengths.sum(-1, squeeze=True) + # + # elif x.Units.equivalent(_units_radians) and y.Units.equivalent( + # _units_radians + # ): + # # ---------------------------------------------------- + # # Spherical lines. + # # + # # Each line segment is a great circle arc between two + # # adjacent nodes. + # # + # # The length of the great circle arc is the the central + # # angle multiplied by the radius. The central angle is + # # calculated with a special case of the Vincenty formula: + # # https://en.wikipedia.org/wiki/Great-circle_distance + # # ---------------------------------------------------- + # if not great_circle: + # raise ValueError( + # "Must set great_circle=True to allow the derivation " + # "of line-length weights from great circle segments." + # ) + # + # if methods: + # comp[(axis,)] = "linear spherical line geometry" + # return True + # + # x.Units = _units_radians + # y.Units = _units_radians + # + # central_angles = self._weights_central_angles(x, y) + # + # # It's worth checking for negative areas, as this has been + # # known to trap problems. TODOUGRID: is it? can we do this + # # check on the weights later? Is it worth it even then? + # if central_angles.min() < 0: + # raise ValueError( + # "A spherical line geometry segment has " + # f"negative length: {central_angles.min() * radius!r}" + # ) + # + # all_lengths = central_angles.sum(-1, squeeze=True) + # + # if measure: + # all_lengths *= radius + # else: + # return False + # + # # Sum the lengths of each part to get the total length of + # # each cell + # lengths = all_lengths.sum(-1, squeeze=True) + # + # comp[(axis,)] = lengths + # + # weights_axes.add(axis) + # + # return True + # + # def _weights_geometry_volume( + # self, + # comp, + # weights_axes, + # auto=False, + # measure=False, + # radius=None, + # great_circle=False, + # methods=False, + # ): + # """Creates volume weights for polygon geometry cells. + # + # .. versionadded:: 3.2.0 + # + # :Parameters: + # + # measure: `bool` + # If True then make sure that the weights represent true + # cell sizes. + # + # methods: `bool`, optional + # If True then add a description of the method used to + # create the weights to the *comp* dictionary, as opposed to + # the actual weights. + # + # """ + # axis, aux_X, aux_Y, aux_Z = self._weights_yyy( + # "polygon", methods=methods, auto=auto + # ) + # + # if axis is None and auto: + # return False + # + # if axis in weights_axes: + # if auto: + # return False + # + # raise ValueError( + # "Multiple weights specifications for " + # f"{self.constructs.domain_axis_identity(axis)!r} axis" + # ) + # + # x = aux_X.bounds.data + # y = aux_Y.bounds.data + # z = aux_Z.bounds.data + # + # if not z.Units.equivalent(_units_metres): + # if auto: + # return False + # + # raise ValueError( + # "Z coordinate bounds must have units equivalent to " + # f"metres for volume calculations. Got {z.Units!r}." + # ) + # + # if not methods: + # # Initialise cell volumes as the cell areas + # volumes = self._weights_geometry_area( + # comp, + # weights_axes, + # auto=auto, + # measure=measure, + # radius=radius, + # great_circle=great_circle, + # methods=False, + # return_areas=True, + # ) + # + # if measure: + # delta_z = abs(z[..., 1] - z[..., 0]) + # delta_z.squeeze(axis=-1, inplace=True) + # + # if x.Units.equivalent(_units_metres) and y.Units.equivalent( + # _units_metres + # ): + # # ---------------------------------------------------- + # # Plane polygons defined by straight lines. + # # + # # Do this in preference to weights based on spherical + # # polygons, which require the great circle assumption. + # # ---------------------------------------------------- + # if methods: + # comp[(axis,)] = "volume plane polygon geometry" + # return True + # + # if measure: + # volumes *= delta_z + # + # elif x.Units.equivalent(_units_radians) and y.Units.equivalent( + # _units_radians + # ): + # # ---------------------------------------------------- + # # Spherical polygons defined by great circles + # # + # # The area of such a spherical polygon is given by the + # # sum of the interior angles minus (N-2)pi, where N is + # # the number of sides (Todhunter): + # # + # # https://en.wikipedia.org/wiki/Spherical_trigonometry#Area_and_spherical_excess + # # + # # The interior angle of a side is calculated with a + # # special case of the Vincenty formula: + # # https://en.wikipedia.org/wiki/Great-circle_distance + # # ---------------------------------------------------- + # if not great_circle: + # raise ValueError( + # "Must set great_circle=True to allow the derivation " + # "of volume weights from spherical polygons composed " + # "from great circle segments." + # ) + # + # if methods: + # comp[(axis,)] = "volume spherical polygon geometry" + # return True + # + # if measure: + # r = radius + # + # # actual_volume = + # # [actual_area/(4*pi*r**2)] + # # * (4/3)*pi*[(r+delta_z)**3 - r**3)] + # volumes *= ( + # delta_z**3 / (3 * r**2) + delta_z**2 / r + delta_z + # ) + # else: + # raise ValueError( + # "X and Y coordinate bounds must both have units " + # "equivalent to either metres, for plane polygon, or radians, " + # "for spherical polygon, volume calculations. Got " + # f"{x.Units!r} and {y.Units!r}." + # ) + # + # comp[(axis,)] = volumes + # + # weights_axes.add(axis) + # + # return True + # + # def _weights_interior_angles(self, lon, lat, interior_ring=None): + # """Find the interior angles at spherical polygon vertices. + # + # The interior angle at a vertex is that formed by two adjacent + # sides. Each interior angle is therefore a function of three + # vertices - the vertex at which the interior angle is required + # and the two adjacent vertices either side of it. + # + # **Method** + # + # Bevis, M., Cambareri, G. Computing the area of a spherical + # polygon of arbitrary shape. Math Geol 19, 335–346 (1987). + # + # http://bemlar.ism.ac.jp/zhuang/Refs/Refs/bevis1987mathgeol.pdf + # + # Abstract: An algorithm for determining the area of a spherical + # polygon of arbitrary shape is presented. The kernel of the + # problem is to compute the interior angle at each vertex of + # the spherical polygon; a well-known relationship between the + # area of a spherical polygon and the sum of its interior + # angles then may be exploited. The algorithm has been + # implemented as a FORTRAN subroutine and a listing is + # provided. + # + # .. versionadded:: UGRIDVER + # + # .. seealso:: `_weights_central_angles` + # + # :Parameters: + # + # lon: `Data` + # Longitudes. Must have units of radians, which is not + # checked. + # + # lat: `Data` + # Latitudes. Must have units of radians, which is not + # checked. + # + # interior_ring: `Data`, optional + # The interior ring flags for geometry cells. + # + # :Returns: + # + # `Data` + # The interior angles of each spherical polygon, in + # units of radians. + # + # """ + # # P denotes a vertex at which the interior angle is required, + # # A denotes the adjacent point clockwise from P, and B denotes + # # the adjacent point anticlockwise from P. + # + # lon_A = lon[..., :-2] + # lon_P = lon[..., 1:-1] + # lon_B = lon[..., 2:] + # + # lon_A_minus_lon_P = lon_A - lon_P + # lon_B_minus_lon_P = lon_B - lon_P + # + # cos_lat = lat.cos() + # cos_lat_A = cos_lat[..., :-2] + # cos_lat_P = cos_lat[..., 1:-1] + # cos_lat_B = cos_lat[..., 2:] + # + # sin_lat = lat.sin() + # sin_lat_A = sin_lat[..., :-2] + # sin_lat_P = sin_lat[..., 1:-1] + # sin_lat_B = sin_lat[..., 2:] + # + # y = lon_A_minus_lon_P.sin() * cos_lat_A + # x = ( + # sin_lat_A * cos_lat_P + # - cos_lat_A * sin_lat_P * lon_A_minus_lon_P.cos() + # ) + # lat_A_primed = self._Data.arctan2(y, x) + # + # y = lon_B_minus_lon_P.sin() * cos_lat_B + # x = ( + # sin_lat_B * cos_lat_P + # - cos_lat_B * sin_lat_P * lon_B_minus_lon_P.cos() + # ) + # lat_B_primed = self._Data.arctan2(y, x) + # + # # The CF vertices here are, in general, given in anticlockwise + # # order, so we do alpha_P = lat_B_primed - lat_A_primed, + # # rather than the "alpha_P = lat_A_primed - lat_B_primed" + # # given in Bevis and Cambareri, which assumes clockwise order. + # alpha_P = lat_B_primed - lat_A_primed + # + # if interior_ring is not None: + # # However, interior rings *are* given in clockwise order + # # in CF, so we need to negate alpha_P in these cases. + # alpha_P.where(interior_ring, -alpha_P, inplace=True) + # + # alpha_P = alpha_P.where(lt(0), alpha_P + 2 * np.pi, alpha_P) + # return alpha_P + # + # def _weights_central_angles(self, lon, lat): + # r"""Find the central angles for spherical great circle line segments. + # + # The central angle of two points on the sphere is the angle + # subtended from the centre of the sphere by the two points on + # its surface points. It is calculated with a special case of + # the Vincenty formula that is accurate for all spherical + # distances: + # + # \Delta \sigma =\arctan {\frac {\sqrt {\left(\cos \phi + # _{2}\sin(\Delta \lambda )\right)^{2}+\left(\cos \phi _{1}\sin + # \phi _{2}-\sin \phi _{1}\cos \phi _{2}\cos(\Delta \lambda + # )\right)^{2}}}{\sin \phi _{1}\sin \phi _{2}+\cos \phi _{1}\cos + # \phi _{2}\cos(\Delta \lambda )}} + # + # The quadrant for \Delta \sigma should be governed by the signs + # of the numerator and denominator of the right hand side using + # the atan2 function. + # + # (https://en.wikipedia.org/wiki/Great-circle_distance) + # + # :Parameters: + # + # lon: `Data` + # Longitudes. Must have units of radians, which is not + # checked. + # + # lat: `Data` + # Latitudes. Must have units of radians, which is not + # checked. + # + # :Returns: + # + # `Data` + # The central angles in units of radians. + # + # """ + # # A and B denote adjacent vertices that define each line + # # segment + # + # delta_lon = lon.diff(axis=-1) + # + # cos_lat = lat.cos() + # sin_lat = lat.sin() + # + # cos_lat_A = cos_lat[..., :-1] + # cos_lat_B = cos_lat[..., 1:] + # + # sin_lat_A = sin_lat[..., :-1] + # sin_lat_B = sin_lat[..., 1:] + # + # cos_delta_lon = delta_lon.cos() + # sin_delta_lon = delta_lon.sin() + # + # y = ( + # (cos_lat_B * sin_delta_lon) ** 2 + # + (cos_lat_A * sin_lat_B - sin_lat_A * cos_lat_B * cos_delta_lon) + # ** 2 + # ) ** 0.5 + # x = sin_lat_A * sin_lat_B + cos_lat_A * cos_lat_B * cos_delta_lon + # + # delta_sigma = self._Data.arctan2(y, x) + # delta_sigma.override_units(_units_1, inplace=True) + # return delta_sigma + # + # def _weights_linear( + # self, + # axis, + # comp, + # weights_axes, + # auto=False, + # measure=False, + # methods=False, + # ): + # """1-d linear weights from dimension coordinate constructs. + # + # :Parameters: + # + # axis: `str` + # + # comp: `dict` + # + # weights_axes: `set` + # + # auto: `bool` + # + # measure: `bool` + # If True then make sure that the weights represent true + # cell sizes. + # + # methods: `bool`, optional + # If True then add a description of the method used to + # create the weights to the *comp* dictionary, as opposed to + # the actual weights. + # + # :Returns: + # + # `bool` + # + # """ + # da_key = self.domain_axis(axis, key=True, default=None) + # if da_key is None: + # if auto: + # return False + # + # raise ValueError( + # "Can't create weights: Can't find domain axis " + # f"matching {axis!r}" + # ) + # + # dim = self.dimension_coordinate(filter_by_axis=(da_key,), default=None) + # if dim is None: + # if auto: + # return False + # + # raise ValueError( + # f"Can't create linear weights for {axis!r} axis: Can't find " + # "dimension coordinate construct." + # ) + # + # if not measure and dim.size == 1: + # return False + # + # if da_key in weights_axes: + # if auto: + # return False + # + # raise ValueError( + # f"Can't create linear weights for {axis!r} axis: Multiple " + # "axis specifications" + # ) + # + # if not dim.has_bounds(): + # # Dimension coordinate has no bounds + # if auto: + # return False + # + # raise ValueError( + # f"Can't create linear weights for {axis!r} axis: No bounds" + # ) + # else: + # # Bounds exist + # if methods: + # comp[ + # (da_key,) + # ] = "linear " + self.constructs.domain_axis_identity(da_key) + # else: + # comp[(da_key,)] = dim.cellsize + # + # weights_axes.add(da_key) + # + # return True + # + # def _weights_measure( + # self, measure, comp, weights_axes, methods=False, auto=False + # ): + # """Cell measure weights. + # + # :Parameters: + # + # methods: `bool`, optional + # If True then add a description of the method used to + # create the weights to the *comp* dictionary, as opposed to + # the actual weights. + # + # :Returns: + # + # `bool` + # + # """ + # m = self.cell_measures(filter_by_measure=(measure,), todict=True) + # len_m = len(m) + # + # if not len_m: + # if measure == "area": + # return False + # + # if auto: + # return + # + # raise ValueError( + # f"Can't find weights: No {measure!r} cell measure" + # ) + # + # elif len_m > 1: + # if auto: + # return False + # + # raise ValueError( + # f"Can't find weights: Multiple {measure!r} cell measures" + # ) + # + # key, clm = m.popitem() + # + # clm_axes0 = self.get_data_axes(key) + # + # clm_axes = tuple( + # [axis for axis, n in zip(clm_axes0, clm.data.shape) if n > 1] + # ) + # + # for axis in clm_axes: + # if axis in weights_axes: + # if auto: + # return False + # + # raise ValueError( + # "Multiple weights specifications for " + # f"{self.constructs.domain_axis_identity(axis)!r} axis" + # ) + # + # clm = clm.get_data(_fill_value=False).copy() + # if clm_axes != clm_axes0: + # iaxes = [clm_axes0.index(axis) for axis in clm_axes] + # clm.squeeze(iaxes, inplace=True) + # + # if methods: + # comp[tuple(clm_axes)] = measure + " cell measure" + # else: + # comp[tuple(clm_axes)] = clm + # + # weights_axes.update(clm_axes) + # + # return True + # + # def _weights_scale(self, w, scale): + # """Scale the weights so that they are <= scale. + # + # :Parameters: + # + # w: `Data` + # The weights to be scaled. + # + # scale: number or `None` + # The maximum value of the scaled weights. If `None` then no + # scaling is applied. + # + # :Returns: + # + # `Data` + # The scaled weights. + # + # """ + # if scale is None: + # return w + # + # if scale <= 0: + # raise ValueError( + # "Can't set 'scale' parameter to a negative number. " + # f"Got {scale!r}" + # ) + # + # w = w / w.max() + # if scale != 1: + # w = w * scale + # + # return w + # + # def _weights_yyy( + # self, domain_axis, geometry_type, methods=False, auto=False + # ): + # """Checks whether weights can be created for given coordinates. + # + # .. versionadded:: 3.2.0 + # + # :Parameters: + # + # domain_axis : `str` or `None` + # + # geometry_type: `str` + # Either ``'polygon'`` or ``'line'``. + # + # auto: `bool` + # + # :Returns: + # + # `tuple` + # + # """ + # aux_X = None + # aux_Y = None + # aux_Z = None + # x_axis = None + # y_axis = None + # z_axis = None + # + # auxiliary_coordinates_1d = self.auxiliary_coordinates( + # filter_by_naxes=(1,), todict=True + # ) + # + # for key, aux in auxiliary_coordinates_1d.items(): + # aux_axis = self.get_data_axes(key)[0] + # dt = self.domain_topology(default=None, filter_by_axis=(aux_axis,)) + # if dt is None and aux.get_geometry(None) != geometry_type: + # continue + # + # # Still here? Then this auxiliary coordinate has either + # # UGRID or geometry cells. + # if aux.X: + # aux_X = aux.copy() + # x_axis = aux_axis # self.get_data_axes(key)[0] + # if domain_axis is not None and x_axis != domain_axis: + # aux_X = None + # continue + # elif aux.Y: + # aux_Y = aux.copy() + # y_axis = aux_axis # self.get_data_axes(key)[0] + # if domain_axis is not None and y_axis != domain_axis: + # aux_Y = None + # continue + # elif aux.Z: + # aux_Z = aux.copy() + # z_axis = aux_axis # self.get_data_axes(key)[0] + # if domain_axis is not None and z_axis != domain_axis: + # aux_Z = None + # continue + # + # if aux_X is None or aux_Y is None: + # if auto: + # return (None,) * 4 + # + # raise ValueError( + # "Can't create weights: Need both X and Y nodes to " + # f"calculate {geometry_type} cell weights" + # ) + # + # if x_axis != y_axis: + # if auto: + # return (None,) * 4 + # + # raise ValueError( + # "Can't create weights: X and Y cells span different " + # "domain axes" + # ) + # + # axis = x_axis + # + # if aux_X.get_bounds(None) is None or aux_Y.get_bounds(None) is None: + # # Not both X and Y coordinates have bounds + # if auto: + # return (None,) * 4 + # + # raise ValueError( + # "Can't create weights: " + # "Not both X and Y coordinates have bounds" + # ) + # + # if aux_X.bounds.shape != aux_Y.bounds.shape: + # if auto: + # return (None,) * 4 + # + # raise ValueError( + # "Can't create weights: UGRID/geometry X and Y coordinate " + # "bounds must have the same shape. " + # f"Got {aux_X.bounds.shape} and {aux_Y.bounds.shape}" + # ) + # + # if aux_Z is None: + # for key, aux in auxiliary_coordinates_1d.items(): + # if aux.Z: + # aux_Z = aux.copy() + # z_axis = self.get_data_axes(key)[0] + # + # # Check Z coordinates + # if aux_Z is not None: + # if z_axis != x_axis: + # if auto: + # return (None,) * 4 + # + # raise ValueError( + # "Z coordinates span different domain axis to X and Y " + # "geometry coordinates" + # ) + # + # return axis, aux_X, aux_Y, aux_Z # ---------------------------------------------------------------- # End of worker functions for weights @@ -4854,6 +4852,19 @@ def weights( (2,): 'linear longitude'} """ + from .weights import ( + weights_area_XY, + weights_data, + weights_field, + weights_field_scalar, + weights_geometry_area, + weights_geometry_line, + weights_geometry_volume, + weights_linear, + weights_measure, + weights_scale, + ) + if isinstance(weights, str) and weights == "auto": _DEPRECATION_ERROR_KWARG_VALUE( self, @@ -4889,7 +4900,8 @@ def weights( return Data(1.0, "1") # Return a field containing a single weight of 1 - return self._weights_field_scalar() + # return self._weights_field_scalar() + return weights_field_scalar(self) # Still here? if methods: @@ -4915,18 +4927,22 @@ def weights( # Auto-detect all weights # -------------------------------------------------------- # Volume weights - if self._weights_measure( - "volume", comp, weights_axes, methods=methods, auto=True + # if self._weights_measure( + if weights_measure( + self, "volume", comp, weights_axes, methods=methods, auto=True ): # Found volume weights from cell measures pass - elif self._weights_measure( - "area", comp, weights_axes, methods=methods, auto=True + # elif self._weights_measure( + elif weights_measure( + self, "area", comp, weights_axes, methods=methods, auto=True ): # Found area weights from cell measures pass - elif self._weights_area_XY( + # elif self._weights_area_XY( + elif weights_area_XY( + self, comp, weights_axes, measure=measure, @@ -4941,7 +4957,9 @@ def weights( domain_axes = self.domain_axes(todict=True) for da_key in domain_axes: - if self._weights_geometry_area( + # if self._weights_geometry_area( + if weights_geometry_area( + self, da_key, comp, weights_axes, @@ -4953,7 +4971,9 @@ def weights( ): # Found area weights from polygon geometries pass - elif self._weights_geometry_line( + # elif self._weights_geometry_line( + elif weights_geometry_line( + self, da_key, comp, weights_axes, @@ -4965,7 +4985,9 @@ def weights( ): # Found linear weights from line geometries pass - elif self._weights_linear( + # elif self._weights_linear( + elif weights_linear( + self, da_key, comp, weights_axes, @@ -5028,13 +5050,16 @@ def weights( # -------------------------------------------------------- # Field # -------------------------------------------------------- - self._weights_field([weights], comp, weights_axes) + # self._weights_field([weights], comp, weights_axes) + weights_field(self, [weights], comp, weights_axes) elif isinstance(weights, Data): # -------------------------------------------------------- # Data # -------------------------------------------------------- - self._weights_data( + # self._weights_data( + weights_data( + self, weights, comp, weights_axes, @@ -5081,22 +5106,37 @@ def weights( axes.append(w) # Field weights - self._weights_field(fields, comp, weights_axes) + # self._weights_field(fields, comp, weights_axes) + weights_field(self, fields, comp, weights_axes) # Volume weights if "volume" in cell_measures: - self._weights_measure( - "volume", comp, weights_axes, methods=methods, auto=False + # self._weights_measure( + weights_measure( + self, + "volume", + comp, + weights_axes, + methods=methods, + auto=False, ) # Area weights if "area" in cell_measures: - if self._weights_measure( - "area", comp, weights_axes, methods=methods, auto=True + # if self._weights_measure( + if weights_measure( + self, + "area", + comp, + weights_axes, + methods=methods, + auto=True, ): # Found area weights from cell measures pass - elif self._weights_area_XY( + # elif self._weights_area_XY( + elif weights_area_XY( + self, comp, weights_axes, measure=measure, @@ -5109,7 +5149,9 @@ def weights( pass else: # Found area weights from 1-d UGRID/geometry cells - self._weights_geometry_area( + # self._weights_geometry_area( + weights_geometry_area( + self, None, comp, weights_axes, @@ -5128,7 +5170,9 @@ def weights( f"{axis!r}" ) - if self._weights_geometry_area( + # if self._weights_geometry_area( + if weights_geometry_area( + self, da_key, comp, weights_axes, @@ -5140,7 +5184,9 @@ def weights( ): # Found area weights from polygon geometries pass - elif self._weights_geometry_line( + # elif self._weights_geometry_line( + elif weights_geometry_line( + self, da_key, comp, weights_axes, @@ -5153,7 +5199,9 @@ def weights( # Found linear weights from line geometries pass else: - self._weights_linear( + # self._weights_linear( + weights_linear( + self, da_key, comp, weights_axes, @@ -5171,10 +5219,13 @@ def weights( del comp[(yaxis,)] weights_axes.discard(xaxis) weights_axes.discard(yaxis) - if not self._weights_measure( - "area", comp, weights_axes, methods=methods + # if not self._weights_measure( + if not weights_measure( + self, "area", comp, weights_axes, methods=methods ): - self._weights_area_XY( + # self._weights_area_XY( + weights_area_XY( + self, comp, weights_axes, measure=measure, @@ -5184,11 +5235,12 @@ def weights( if not methods: if scale is not None: - # -------------------------------------------------------- + # ---------------------------------------------------- # Scale the weights so that they are <= scale - # -------------------------------------------------------- + # ---------------------------------------------------- for key, w in comp.items(): - comp[key] = self._weights_scale(w, scale) + # comp[key] = self._weights_scale(w, scale) + comp[key] = weights_scale(w, scale) for w in comp.values(): mn = w.minimum() @@ -5219,7 +5271,8 @@ def weights( # No component weights have been defined so return an # equal weights field # -------------------------------------------------------- - f = self._weights_field_scalar() + # f = self._weights_field_scalar() + f = weights_field_scalar(self) if data: return f.data @@ -5240,7 +5293,8 @@ def weights( # -------------------------------------------------------- # Scale the weights so that they are <= scale # -------------------------------------------------------- - wdata = self._weights_scale(wdata, scale) + # wdata = self._weights_scale(wdata, scale) + wdata = weights_scale(wdata, scale) # ------------------------------------------------------------ # Reorder the data so that its dimensions are in the same diff --git a/cf/regrid/regrid.py b/cf/regrid/regrid.py index 3f01f0c6dd..c558863b9e 100644 --- a/cf/regrid/regrid.py +++ b/cf/regrid/regrid.py @@ -61,8 +61,12 @@ class Grid: # 'domainaxis2', 'Y': 'domainaxis1'} or ['domainaxis1', # 'domainaxis2'], or {'X': 'domainaxis2', 'Y': 'domainaxis2'} axes: Any = None - # + # The number of domain axis identifiers of the regrid axes. n_regrid_axes: int = 0 + # The dimensionality of the regridding on this grid. Generally the + # same as *n_regrid_axes*, but for a mesh axis *n_regrid_axes* is + # 1 and *regridding_dimensionality* is 2. + regridding_dimensionality: int = 0 # The shape of the regridding axes, in the same order as the # 'axis_keys' attribute. E.g. (96, 73) or (1243,) shape: tuple = None @@ -302,15 +306,7 @@ def regrid( f"using the {method!r} regridding method." ) - # Parse the axes - if spherical: - n_axes = 2 - elif not create_regrid_operator: - if regrid_operator.dst_mesh_location: - n_axes = 2 - else: - n_axes = len(dst_axes) - elif cartesian: + if cartesian: if isinstance(axes, str): axes = (axes,) @@ -332,18 +328,6 @@ def regrid( elif isinstance(dst_axes, str): dst_axes = (dst_axes,) - # if src_mesh and src_mesh_axes == tuple(src_axes): - # # The source regrid axis is a mesh topology axis - # two_d_regridding = True # n_axes = 2 - # else: - n_axes = len(src_axes) - - if n_axes != 2 and method in ("conservative_2nd", "patch"): - raise ValueError( - f"{method!r} regridding is only available for 2-d regridding, " - f"but {n_axes}-d regridding has been requested." - ) - if isinstance(dst, src._Domain): use_dst_mask = False dst = dst.copy() @@ -399,6 +383,20 @@ def regrid( conform_coordinate_units(src_grid, dst_grid) + if method in ("conservative_2nd", "patch") and not ( + src_grid.regridding_dimensionality == 2 + and dst_grid.regridding_dimensionality == 2 + ): + raise ValueError( + f"{method!r} regridding is only available for 2-d regridding" + ) + + # if dst_grid.mesh_location and not src_grid.mesh_location: + # raise ValueError( + # "Can't (yet!) regrid from a non-mesh source grid to a mesh " + # "destination grid" + # ) + if create_regrid_operator: # ------------------------------------------------------------ # Create a new regrid operator @@ -529,12 +527,24 @@ def regrid( # ---------------------------------------------------------------- # Still here? Then do the regridding # ---------------------------------------------------------------- - regridded_axis_sizes = { - src_iaxis: (dst_size,) - for src_iaxis, dst_size in zip(src_grid.axis_indices, dst_grid.shape) - } - if src_grid.n_regrid_axes == 1 and dst_grid.n_regrid_axes > 1: - regridded_axis_sizes[src_grid.axis_indices[0]] = dst_grid.shape + if src_grid.n_regrid_axes == dst_grid.n_regrid_axes: + regridded_axis_sizes = { + src_iaxis: (dst_size,) + for src_iaxis, dst_size in zip( + src_grid.axis_indices, dst_grid.shape + ) + } + elif src_grid.n_regrid_axes == 1: + # Fewer source grid axes than destination grid axes (e.g. mesh + # regridded to lat/lon). + regridded_axis_sizes = {src_grid.axis_indices[0]: dst_grid.shape} + elif dst_grid.n_regrid_axes == 1: + # More source grid axes than destination grid axes + # (e.g. lat/lon regridded to mesh). + src_axis_indices = sorted(src_grid.axis_indices) + regridded_axis_sizes = {src_axis_indices[0]: (dst_grid.shape[0],)} + for src_iaxis in src_axis_indices[1:]: + regridded_axis_sizes[src_iaxis] = () regridded_data = src.data._regrid( method=method, @@ -547,11 +557,9 @@ def regrid( # ---------------------------------------------------------------- # Update the regridded metadata # ---------------------------------------------------------------- - update_non_coordinates( - src, dst, regrid_operator, src_grid=src_grid, dst_grid=dst_grid - ) + update_non_coordinates(src, dst, regrid_operator, src_grid, dst_grid) - update_coordinates(src, dst, src_grid=src_grid, dst_grid=dst_grid) + update_coordinates(src, dst, src_grid, dst_grid) # ---------------------------------------------------------------- # Insert regridded data into the new field @@ -907,7 +915,7 @@ def spherical_grid(f, name=None, method=None, cyclic=None, axes=None): aux_coords_1d = True x_axis = mesh_axis y_axis = mesh_axis - elif mesh_location in ("edge", "volume"): + elif mesh_location: raise ValueError( f"Can't regrid {'from' if name == 'source' else 'to'} " f"a {name} grid comprising an unstructured mesh of " @@ -981,11 +989,17 @@ def spherical_grid(f, name=None, method=None, cyclic=None, axes=None): data_axes = f.get_data_axes() axis_indices = [data_axes.index(key) for key in axis_keys] + n_regrid_axes = len(axis_keys) + regridding_dimensionality = n_regrid_axes + if mesh_location: + regridding_dimensionality += 1 + return Grid( axis_keys=axis_keys, axis_indices=axis_indices, axes=axes, - n_regrid_axes=len(axis_keys), + n_regrid_axes=n_regrid_axes, + regridding_dimensionality=regridding_dimensionality, shape=shape, coords=coords, bounds=get_bounds(method, coords, mesh_location), @@ -1057,15 +1071,31 @@ def Cartesian_grid(f, name=None, method=None, axes=None): axis_keys.append(key) axis_sizes.append(domain_axis.size) - # mesh = False - # domain_topology, mesh_location, mesh_axis = get_mesh(f) - # if mesh_axis is not None: - # if [mesh_axis] == axis_keys: - # mesh = True - # elif mesh_axis in axis_keys: - # raise ValueError( - # "TODOUGRID Can't regrid a combo of mesh and non-mesh axes" - # ) + domain_topology, mesh_location, mesh_axis = get_mesh(f) + if mesh_location: + # There is a domain topology axis + if list(set(axis_keys)) == [mesh_axis]: + # There is a unique regridding axis, and its the domain + # topology axis. + if mesh_location not in ("face", "node"): + raise ValueError( + f"Can't regrid {'from' if name == 'source' else 'to'} " + f"a {name} grid comprising an unstructured mesh of " + f"{mesh_location!r} cells" + ) + + axis_keys = axis_keys[0:1] + axis_sizes = axis_sizes[0:1] + elif mesh_axis in axis_keys: + raise ValueError( + "Can't do Cartesian regridding for a combination of " + f"mesh and non-mesh axes: {axis_keys}" + ) + else: + # None of the regridding axes have a domain topology + domain_topology = None + mesh_location = None + mesh_axis = None if f.construct_type == "domain": axis_indices = list(range(len(axis_keys))) @@ -1082,25 +1112,40 @@ def Cartesian_grid(f, name=None, method=None, axes=None): data_axes = f.get_data_axes() axis_indices = [data_axes.index(key) for key in axis_keys] - domain_topology, mesh_location, mesh_axis = get_mesh(f) - cyclic = False coords = [] - if domain_topology is not None: - # TODOUGRID: check for the mesh spanning the two axes - # ... somehow. Also maybe set mesh=True and cyclic=None. - # - pass + if mesh_location: + if n_axes == 1: + coord_ids = ["X", "Y"] + elif n_axes == 2: + coord_ids = axes[::-1] + else: + raise ValueError("Can't provide 3 axes for mesh axis regridding") + + for coord_id in coord_ids: + aux = f.auxiliary_coordinate( + coord_id, + filter_by_axis=(mesh_axis,), + axis_mode="exact", + default=None, + ) + if aux is None: + raise ValueError( + f"Could not find {coord_id!r} 1-d mesh auxiliary " + "coordinates" + ) + + coords.append(aux) else: for key in axis_keys[::-1]: - coord = f.dimension_coordinate(filter_by_axis=(key,), default=None) - if coord is None: + dim = f.dimension_coordinate(filter_by_axis=(key,), default=None) + if dim is None: raise ValueError( f"No unique {name} dimension coordinate for domain axis " f"{key!r}." ) - coords.append(coord) + coords.append(dim) bounds = get_bounds(method, coords, mesh_location) @@ -1119,11 +1164,17 @@ def Cartesian_grid(f, name=None, method=None, axes=None): coords.append(data) dummy_size_2_dimension = True + n_regrid_axes = len(axis_keys) + regridding_dimensionality = n_regrid_axes + if mesh_location: + regridding_dimensionality += 1 + return Grid( axis_keys=axis_keys, axis_indices=axis_indices, axes=axis_keys, - n_regrid_axes=len(axis_keys), + n_regrid_axes=n_regrid_axes, + regridding_dimensionality=regridding_dimensionality, shape=tuple(axis_sizes), coords=coords, bounds=bounds, @@ -1584,7 +1635,7 @@ def create_esmpy_mesh(grid, mask=None): parametric_dim=2, spatial_dim=2, coord_sys=coord_sys ) - element_conn = grid.domain_topology.array + element_conn = grid.domain_topology.normalise().array element_count = element_conn.shape[0] element_types = np.ma.count(element_conn, axis=1) element_conn = np.ma.compressed(element_conn) @@ -2063,7 +2114,7 @@ def update_coordinates(src, dst, src_grid, dst_grid): axis; and if the destination grid is a mesh copy omain topology and cell connectivity constructs from the destination grid. - .. versionadded:: TODOUGRID + .. versionadded:: 3.14.0 :Parameters: @@ -2158,6 +2209,8 @@ def update_non_coordinates( ): """Update the coordinate references of the regridded field. + .. versionadded:: 3.14.0 + :Parameters: src: `Field` diff --git a/cf/test/test_general.py b/cf/test/test_general.py index a4fd53d035..c6aed83cdf 100644 --- a/cf/test/test_general.py +++ b/cf/test/test_general.py @@ -43,29 +43,6 @@ def test_GENERAL(self): f == c - # +, -, *, /, ** - h = g.copy() - h **= 2 - h **= 0.5 - h.standard_name = g.standard_name - self.assertTrue(g.data.allclose(h.data), repr(g.array - h.array)) - h *= 10 - h /= 10.0 - self.assertTrue(g.data.allclose(h.data), repr(g.array - h.array)) - h += 1 - h -= 1 - self.assertTrue(g.data.allclose(h.data), repr(g.array - h.array)) - h = h**2.0 - h = h**0.5 - h.standard_name = g.standard_name - self.assertTrue(g.data.allclose(h.data), repr(g.array - h.array)) - h = h * 10 - h = h / 10.0 - self.assertTrue(g.data.allclose(h.data), repr(g.array - h.array)) - h = h + 1 - h = h - 1 - self.assertTrue(g.data.allclose(h.data), repr(g.array - h.array)) - # flip, expand_dims, squeeze and remove_axes h = g.copy() h.flip((1, 0), inplace=True) diff --git a/cf/test/test_read_write.py b/cf/test/test_read_write.py index 0720ef4173..c6a46b46cc 100644 --- a/cf/test/test_read_write.py +++ b/cf/test/test_read_write.py @@ -663,6 +663,7 @@ def test_read_CDL(self): geometry_1_file = os.path.join( os.path.dirname(os.path.abspath(__file__)), "geometry_1.nc" ) + tmpfileh2 = "delme.cdl" # TODOUGRID subprocess.run( " ".join(["ncdump", "-h", geometry_1_file, ">", tmpfileh2]), shell=True, @@ -682,7 +683,9 @@ def test_read_CDL(self): c = cf.read(tmpfilec)[0] # Case (2) as above, so the right error should be raised on read - with self.assertRaises(ValueError): + with self.assertRaises( + ValueError + ): # TODOUGRID: why does this raise? should it? What's changed? cf.read(tmpfileh2)[0] with self.assertRaises(ValueError): diff --git a/cf/weights.py b/cf/weights.py new file mode 100644 index 0000000000..8651d37c8a --- /dev/null +++ b/cf/weights.py @@ -0,0 +1,1428 @@ +"""Worker functions for cf.Field.weights""" +import numpy as np + +from .units import Units + + +def weights_area_XY( + f, + comp, + weights_axes, + auto=False, + measure=False, + radius=None, + methods=False, +): + """Calculate area weights from X and Y dimension coordinate + constructs. + + .. versionadded:: UGRIDVER + + :Parameters: + + measure: `bool` + If True then make sure that the weights represent true + cell sizes. + + methods: `bool`, optional + If True then add a description of the method used to + create the weights to the *comp* dictionary, as opposed to + the actual weights. + + :Returns: + + `bool` or `None` + + """ + xkey, xcoord = f.dimension_coordinate("X", item=True, default=(None, None)) + ykey, ycoord = f.dimension_coordinate("Y", item=True, default=(None, None)) + + if xcoord is None or ycoord is None: + if auto: + return + + raise ValueError( + "No unique 'X' and 'Y' dimension coordinate constructs for " + "calculating area weights" + ) + + radians = Units("radians") + metres = Units("metres") + + if xcoord.Units.equivalent(radians) and ycoord.Units.equivalent(radians): + pass + elif xcoord.Units.equivalent(metres) and ycoord.Units.equivalent(metres): + radius = None + else: + if auto: + return + + raise ValueError( + "Insufficient coordinate constructs for calculating " + "area weights" + ) + + xaxis = f.get_data_axes(xkey)[0] + yaxis = f.get_data_axes(ykey)[0] + + for axis in (xaxis, yaxis): + if axis in weights_axes: + if auto: + return + + raise ValueError( + "Multiple weights specifications for " + f"{f.constructs.domain_axis_identity(axis)!r} axis" + ) + + if measure and radius is not None: + radius = f.radius(default=radius) + + if measure or xcoord.size > 1: + if not xcoord.has_bounds(): + if auto: + return + + raise ValueError( + "Can't create area weights: No bounds for " + f"{xcoord.identity()!r} axis" + ) + + if methods: + comp[(xaxis,)] = "linear " + xcoord.identity() + else: + cells = xcoord.cellsize + if xcoord.Units.equivalent(radians): + cells.Units = radians + if measure: + cells *= radius + cells.override_units(radius.Units, inplace=True) + else: + cells.Units = metres + + comp[(xaxis,)] = cells + + weights_axes.add(xaxis) + + if measure or ycoord.size > 1: + if not ycoord.has_bounds(): + if auto: + return + + raise ValueError( + "Can't create area weights: No bounds for " + f"{ycoord.identity()!r} axis" + ) + + if ycoord.Units.equivalent(radians): + ycoord = ycoord.clip(-90, 90, units=Units("degrees")) + ycoord.sin(inplace=True) + + if methods: + comp[(yaxis,)] = "linear sine " + ycoord.identity() + else: + cells = ycoord.cellsize + if measure: + cells *= radius + + comp[(yaxis,)] = cells + else: + if methods: + comp[(yaxis,)] = "linear " + ycoord.identity() + else: + cells = ycoord.cellsize + comp[(yaxis,)] = cells + + weights_axes.add(yaxis) + + return True + + +def weights_data( + f, + w, + comp, + weights_axes, + axes=None, + data=False, + components=False, + methods=False, +): + """Creates weights for the field construct's data array. + + .. versionadded:: UGRIDVER + + :Parameters: + + methods: `bool`, optional + If True then add a description of the method used to + create the weights to the *comp* dictionary, as opposed to + the actual weights. + + """ + # ---------------------------------------------------------------- + # Data weights + # ---------------------------------------------------------------- + field_data_axes = f.get_data_axes() + + if axes is not None: + if isinstance(axes, (str, int)): + axes = (axes,) + + if len(axes) != w.ndim: + raise ValueError( + "'axes' parameter must provide an axis identifier " + f"for each weights data dimension. Got {axes!r} for " + f"{w.ndim} dimension(s)." + ) + + iaxes = [ + field_data_axes.index(f.domain_axis(axis, key=True)) + for axis in axes + ] + + if data: + for i in range(f.ndim): + if i not in iaxes: + w = w.insert_dimension(position=i) + iaxes.insert(i, i) + + w = w.transpose(iaxes) + + if w.ndim > 0: + while w.shape[0] == 1: + w = w.squeeze(0) + + else: + iaxes = list(range(f.ndim - w.ndim, f.ndim)) + + if not (components or methods): + if not f._is_broadcastable(w.shape): + raise ValueError( + f"The 'Data' weights (shape {w.shape}) are not " + "broadcastable to the field construct's data " + f"(shape {f.shape})." + ) + + axes0 = field_data_axes[f.ndim - w.ndim :] + else: + axes0 = [field_data_axes[i] for i in iaxes] + + for axis0 in axes0: + if axis0 in weights_axes: + raise ValueError( + "Multiple weights specified for " + f"{f.constructs.domain_axis_identity(axis0)!r} axis" + ) + + if methods: + comp[tuple(axes0)] = "custom data" + else: + comp[tuple(axes0)] = w + + weights_axes.update(axes0) + + return True + + +def weights_field(f, fields, comp, weights_axes, methods=False): + """Creates a weights field. + + .. versionadded:: UGRIDVER + + """ + s = f.analyse_items() + + domain_axes = f.domain_axes(todict=True) + for w in fields: + t = w.analyse_items() + # TODO CHECK this with org + + if t["undefined_axes"]: + w_domain_axes_1 = w.domain_axes(filter_by_size=(1,), todict=True) + if set(w_domain_axes_1).intersection(t["undefined_axes"]): + raise ValueError( + f"Weights field {w} has at least one undefined " + f"domain axes: {w_domain_axes_1}." + ) + + w = w.squeeze() + + w_domain_axes = w.domain_axes(todict=True) + + axis1_to_axis0 = {} + + coordinate_references = f.coordinate_references(todict=True) + w_coordinate_references = w.coordinate_references(todict=True) + + for axis1 in w.get_data_axes(): + identity = t["axis_to_id"].get(axis1, None) + if identity is None: + raise ValueError( + "Weights field has unmatched, size > 1 " + f"{w.constructs.domain_axis_identity(axis1)!r} axis" + ) + + axis0 = s["id_to_axis"].get(identity, None) + if axis0 is None: + raise ValueError( + f"Weights field has unmatched, size > 1 {identity!r} " + "axis" + ) + + w_axis_size = w_domain_axes[axis1].get_size() + f_axis_size = domain_axes[axis0].get_size() + + if w_axis_size != f_axis_size: + raise ValueError( + f"Weights field has incorrectly sized {identity!r} " + f"axis ({w_axis_size} != {f_axis_size})" + ) + + axis1_to_axis0[axis1] = axis0 + + # Check that the defining coordinate data arrays are + # compatible + key0 = s["axis_to_coord"][axis0] + key1 = t["axis_to_coord"][axis1] + + if not f._equivalent_construct_data( + w, key0=key0, key1=key1, s=s, t=t + ): + raise ValueError( + f"Weights field has incompatible {identity!r} " + "coordinates" + ) + + # Still here? Then the defining coordinates have + # equivalent data arrays + + # If the defining coordinates are attached to coordinate + # references then check that those coordinate references + # are equivalent + refs0 = [ + key + for key, ref in coordinate_references.items() + if key0 in ref.coordinates() + ] + refs1 = [ + key + for key, ref in w_coordinate_references.items() + if key1 in ref.coordinates() + ] + + nrefs = len(refs0) + if nrefs > 1 or nrefs != len(refs1): + # The defining coordinate are associated with + # different numbers of coordinate references + equivalent_refs = False + elif not nrefs: + # Neither defining coordinate is associated with a + # coordinate reference + equivalent_refs = True + else: + # Each defining coordinate is associated with exactly + # one coordinate reference + equivalent_refs = f._equivalent_coordinate_references( + w, key0=refs0[0], key1=refs1[0], s=s, t=t + ) + + if not equivalent_refs: + raise ValueError( + "Input weights field has an incompatible " + "coordinate reference" + ) + + axes0 = tuple([axis1_to_axis0[axis1] for axis1 in w.get_data_axes()]) + + for axis0 in axes0: + if axis0 in weights_axes: + raise ValueError( + "Multiple weights specified for " + f"{f.constructs.domain_axis_identity(axis0)!r} " + "axis" + ) + + comp[tuple(axes0)] = w.data + + weights_axes.update(axes0) + + return True + + +def weights_field_scalar(f, methods=False): + """Return a scalar field of weights with long_name ``'weight'``. + + .. versionadded:: UGRIDVER + + :Returns: + + `Field` + + """ + data = f._Data(1.0, "1") + + f = type(f)() + f.set_data(data, copy=False) + f.long_name = "weight" + f.comment = f"Weights for {f!r}" + + return f + + +def weights_geometry_area( + f, + domain_axis, + comp, + weights_axes, + auto=False, + measure=False, + radius=None, + great_circle=False, + return_areas=False, + methods=False, +): + """Creates area weights for polygon geometry cells. + + .. versionadded:: UGRIDVER + + :Parameters: + + domain_axis : `str` or `None` + + measure: `bool` + If True then make sure that the weights represent true + cell sizes. + + methods: `bool`, optional + If True then add a description of the method used to + create the weights to the *comp* dictionary, as opposed to + the actual weights. + + :Returns: + + `bool` or `Data` + + """ + axis, aux_X, aux_Y, aux_Z = weights_yyy( + f, domain_axis, "polygon", methods=methods, auto=auto + ) + + if axis is None: + if auto: + return False + + if domain_axis is None: + raise ValueError("No polygon geometries") + + raise ValueError( + "No polygon geometries for " + f"{f.constructs.domain_axis_identity(domain_axis)!r} axis" + ) + + if axis in weights_axes: + if auto: + return False + + raise ValueError( + "Multiple weights specifications for " + f"{f.constructs.domain_axis_identity(axis)!r} axis" + ) + + ugrid = f.domain_topologies(filter_by_axis=(axis,), todict=True) + + # Check for interior rings + if ugrid: + interior_ring = None + else: + interior_ring_X = aux_X.get_interior_ring(None) + interior_ring_Y = aux_Y.get_interior_ring(None) + if interior_ring_X is None and interior_ring_Y is None: + interior_ring = None + elif interior_ring_X is None: + raise ValueError( + "Can't find weights: X coordinates have missing " + "interior ring variable" + ) + elif interior_ring_Y is None: + raise ValueError( + "Can't find weights: Y coordinates have missing " + "interior ring variable" + ) + elif not interior_ring_X.data.equals(interior_ring_Y.data): + raise ValueError( + "Can't find weights: Interior ring variables for X and Y " + "coordinates have different data values" + ) + else: + interior_ring = interior_ring_X.data + if interior_ring.shape != aux_X.bounds.shape[:-1]: + raise ValueError( + "Can't find weights: Interior ring variables have " + f"incorrect shape. Got {interior_ring.shape}, " + f"expected {aux_X.bounds.shape[:-1]}" + # TODOUGRID: is this not tested in set_interior_ring? + ) + + x = aux_X.bounds.data + y = aux_Y.bounds.data + spherical = False + + radians = Units("radians") + metres = Units("metres") + + if x.Units.equivalent(metres) and y.Units.equivalent(metres): + # -------------------------------------------------------- + # Plane polygons defined by straight lines. + # + # The area is given by the shoelace formula: + # https://en.wikipedia.org/wiki/Shoelace_formula + # + # Do this in preference to weights based on spherical + # polygons, which require the great circle assumption. + # -------------------------------------------------------- + if methods: + comp[(axis,)] = "area plane polygon geometry" + return True + + y.Units = x.Units + + all_areas = (x[..., :-1] * y[..., 1:]).sum(-1, squeeze=True) - ( + x[..., 1:] * y[..., :-1] + ).sum(-1, squeeze=True) + + # TODO: replace this loop! + for i, (parts_x, parts_y) in enumerate(zip(x, y)): + for j, (nodes_x, nodes_y) in enumerate(zip(parts_x, parts_y)): + nodes_x = nodes_x.compressed() + nodes_y = nodes_y.compressed() + + if (nodes_x.size and nodes_x[0] != nodes_x[-1]) or ( + nodes_y.size and nodes_y[0] != nodes_y[-1] + ): + # First and last nodes of this polygon part are + # different => need to account for the "last" edge + # of the polygon that joins the first and last + # points. + all_areas[i, j] += x[-1] * y[0] - x[0] * y[-1] + + all_areas = all_areas.abs() * 0.5 + + elif x.Units.equivalent(radians) and y.Units.equivalent(radians): + # -------------------------------------------------------- + # Spherical polygons defined by great circles + # + # The area is given by the sum of the interior angles minus + # (N-2)pi, where N is the number of sides: + # https://en.wikipedia.org/wiki/Spherical_trigonometry + # -------------------------------------------------------- + spherical = True + + if not great_circle: + raise ValueError( + "Must set great_circle=True to allow the derivation of " + "area weights from spherical polygons composed from " + "great circle segments." + ) + + if methods: + comp[(axis,)] = "area spherical polygon geometry" + return True + + x.Units = radians + y.Units = radians + x = x.array + y = y.array + + cols = np.ma.count(x, axis=1) + n_cells = cols.size + if (cols != np.ma.count(y, axis=1)).any(): + raise ValueError( + "Can't create area weights for " + f"{f.constructs.domain_axis_identity(axis)!r} axis: " + f"{aux_X!r} and {aux_Y!r} have inconsistent bounds " + "specifications" + ) + + col0 = cols[0] + if (cols == col0).all(): + # All cells have the same number of nodes, so we can use + # an integer and a slice in place of two integer arrays, + # which is much faster. + cols = col0 + rows = slice(None) + else: + rows = np.arange(x.shape[0]) + + if not ugrid: + # Ascertain whether or not the first and last boundary + # vertices are coincident in all cells. Note that this can + # never be the case for UGRID cells. + cols_m1 = cols - 1 + duplicate = f._Data(x[:, 0]).isclose(x[rows, cols_m1]) & f._Data( + y[:, 0] + ).isclose(y[rows, cols_m1]) + + empty_col_x = np.ma.masked_all((n_cells, 1), dtype=x.dtype) + empty_col_y = np.ma.masked_all((n_cells, 1), dtype=y.dtype) + + # Pad out the boundary vertiex array for each cell with first + # and last bounds neighbours + if ugrid or not duplicate.any(): + # The first and last boundary vertices are different in + # all cells, i.e. No. nodes = No. edges. + # + # E.g. for triangle cells: + # [[4, 5, 6]] -> [[6, 4, 5, 6, 4]] + # [[4, 5, 6, --]] -> [[6, 4, 5, 6, 4, --]] + cols_p1 = cols + 1 + x = np.ma.concatenate((empty_col_x, x, empty_col_x), axis=1) + y = np.ma.concatenate((empty_col_y, y, empty_col_y), axis=1) + x[:, 0] = x[rows, cols] + y[:, 0] = y[rows, cols] + x[rows, cols_p1] = x[rows, 1] + y[rows, cols_p1] = y[rows, 1] + x = f._Data(x, units=radians) + y = f._Data(y, units=radians) + + # The number of edges of each polygon + N = cols + elif duplicate.all(): + # The first and last boundary vertices coincide in all + # cells, i.e. No. nodes = No. edges + 1. + # + # E.g. for triangle cells: + # [[4, 5, 6, 4]] -> [[6, 4, 5, 6, 4]] + # [[4, 5, 6, 4, --]] -> [[6, 4, 5, 6, 4, --]] + x = np.ma.concatenate((empty_col_x, x), axis=1) + y = np.ma.concatenate((empty_col_y, y), axis=1) + x[:, 0] = x[rows, cols_m1] + y[:, 0] = y[rows, cols_m1] + x = f._Data(x, units=radians) + y = f._Data(y, units=radians) + + # The number of edges of each polygon + N = cols - 1 + else: + raise ValueError( + "Can't calculate spherical polygon cell areas when " + "the first and last boundary vertices coincide in " + "some cells but not all" + ) + + # Find the interior angles of each polygon + interior_angles = weights_interior_angles(f, x, y, interior_ring) + + # Find the polgon areas + all_areas = interior_angles.sum(-1, squeeze=True) - (N - 2) * np.pi + + # It's worth checking for negative areas, as this has been + # known to trap problems. TODOUGRID: is it? can we do this + # check on the weights later? Is it worth it even then? + if all_areas.min() < 0: + if ugrid: + raise ValueError( + "A UGRID spherical polygon cell has negative area" + ) + else: + raise ValueError( + "A spherical polygon geometry cell part has " + "negative area" + ) + else: + return False + + if ugrid: + areas = all_areas + else: + if interior_ring is not None: + # Change the sign of areas for geometry polygons that + # are interior rings + all_areas.where(interior_ring, -all_areas, inplace=True) + + # Sum the areas of each geometry polygon part to get the + # total area of each cell + areas = all_areas.sum(-1, squeeze=True) + + areas.override_units(Units("1"), inplace=True) + + if measure and spherical: + # Multiply by radius squared, accounting for any Z + # coordinates, to get the actual area + if aux_Z is None: + # No Z coordinates + r = radius + else: + z = aux_Z.get_data(None, _fill_value=False) + if z is None: + # No Z coordinates + r = radius + else: + if not z.Units.equivalent(metres): + raise ValueError( + "Z coordinates must have units equivalent to " + f"metres for area calculations. Got {z.Units!r}" + ) + + positive = aux_Z.get_property("positive", None) + if positive is None: + raise ValueError( + "Z coordinate 'positive' property is not defined" + ) + + if positive.lower() == "up": + r = radius + z + elif positive.lower() == "down": + r = radius - z + else: + raise ValueError( + "Bad value of Z coordinate 'positive' property: " + f"{positive!r}." + ) + + areas *= r**2 + + if return_areas: + return areas + + comp[(axis,)] = areas + + weights_axes.add(axis) + + return True + + +def weights_geometry_line( + f, + domain_axis, + comp, + weights_axes, + auto=False, + measure=False, + radius=None, + great_circle=False, + methods=False, +): + """Creates line-length weights for line geometries. + + .. versionadded:: UGRIDVER + + :Parameters: + + measure: `bool` + If True then make sure that the weights represent true + cell sizes. + + methods: `bool`, optional + If True then add a description of the method used to + create the weights to the *comp* dictionary, as opposed to + the actual weights. + + """ + axis, aux_X, aux_Y, aux_Z = weights_yyy( + f, domain_axis, "line", methods=methods, auto=auto + ) + + if axis is None: + if auto: + return False + + if domain_axis is None: + raise ValueError("No line geometries") + + raise ValueError( + "No line geometries for " + f"{f.constructs.domain_axis_identity(domain_axis)!r} axis" + ) + + if axis in weights_axes: + if auto: + return False + + raise ValueError( + "Multiple weights specifications for " + f"{f.constructs.domain_axis_identity(axis)!r} axis" + ) + + x = aux_X.bounds.data + y = aux_Y.bounds.data + + radians = Units("radians") + metres = Units("metres") + + if x.Units.equivalent(metres) and y.Units.equivalent(metres): + # ------------------------------------------------------------ + # Plane lines. + # + # Each line segment is the simple cartesian distance between + # two adjacent nodes. + # ------------------------------------------------------------ + if methods: + comp[(axis,)] = "linear plane line geometry" + return True + + y.Units = x.Units + + delta_x = x.diff(axis=-1) + delta_y = y.diff(axis=-1) + + all_lengths = (delta_x**2 + delta_y**2) ** 0.5 + all_lengths = all_lengths.sum(-1, squeeze=True) + + elif x.Units.equivalent(radians) and y.Units.equivalent(radians): + # ------------------------------------------------------------ + # Spherical lines. + # + # Each line segment is a great circle arc between two adjacent + # nodes. + # + # The length of the great circle arc is the the central angle + # multiplied by the radius. The central angle is calculated + # with a special case of the Vincenty formula: + # https://en.wikipedia.org/wiki/Great-circle_distance + # ------------------------------------------------------------ + if not great_circle: + raise ValueError( + "Must set great_circle=True to allow the derivation " + "of line-length weights from great circle segments." + ) + + if methods: + comp[(axis,)] = "linear spherical line geometry" + return True + + x.Units = radians + y.Units = radians + + central_angles = weights_central_angles(f, x, y) + + # It's worth checking for negative areas, as this has been + # known to trap problems. TODOUGRID: is it? can we do this + # check on the weights later? Is it worth it even then? + if central_angles.min() < 0: + raise ValueError( + "A spherical line geometry segment has " + f"negative length: {central_angles.min() * radius!r}" + ) + + all_lengths = central_angles.sum(-1, squeeze=True) + + if measure: + all_lengths *= radius + else: + return False + + # Sum the lengths of each part to get the total length of each + # cell + lengths = all_lengths.sum(-1, squeeze=True) + + comp[(axis,)] = lengths + + weights_axes.add(axis) + + return True + + +def weights_geometry_volume( + f, + comp, + weights_axes, + auto=False, + measure=False, + radius=None, + great_circle=False, + methods=False, +): + """Creates volume weights for polygon geometry cells. + + .. versionadded:: UGRIDVER + + :Parameters: + + measure: `bool` + If True then make sure that the weights represent true + cell sizes. + + methods: `bool`, optional + If True then add a description of the method used to + create the weights to the *comp* dictionary, as opposed to + the actual weights. + + """ + axis, aux_X, aux_Y, aux_Z = weights_yyy( + f, "polygon", methods=methods, auto=auto + ) + + if axis is None and auto: + return False + + if axis in weights_axes: + if auto: + return False + + raise ValueError( + "Multiple weights specifications for " + f"{f.constructs.domain_axis_identity(axis)!r} axis" + ) + + x = aux_X.bounds.data + y = aux_Y.bounds.data + z = aux_Z.bounds.data + + radians = Units("radians") + metres = Units("metres") + + if not z.Units.equivalent(metres): + if auto: + return False + + raise ValueError( + "Z coordinate bounds must have units equivalent to " + f"metres for volume calculations. Got {z.Units!r}." + ) + + if not methods: + # Initialise cell volumes as the cell areas + volumes = weights_geometry_area( + f, + comp, + weights_axes, + auto=auto, + measure=measure, + radius=radius, + great_circle=great_circle, + methods=False, + return_areas=True, + ) + + if measure: + delta_z = abs(z[..., 1] - z[..., 0]) + delta_z.squeeze(axis=-1, inplace=True) + + if x.Units.equivalent(metres) and y.Units.equivalent(metres): + # ------------------------------------------------------------ + # Plane polygons defined by straight lines. + # + # Do this in preference to weights based on spherical + # polygons, which require the great circle assumption. + # ------------------------------------------------------------ + if methods: + comp[(axis,)] = "volume plane polygon geometry" + return True + + if measure: + volumes *= delta_z + + elif x.Units.equivalent(radians) and y.Units.equivalent(radians): + # ------------------------------------------------------------ + # Spherical polygons defined by great circles + # + # The area of such a spherical polygon is given by the sum of + # the interior angles minus (N-2)pi, where N is the number of + # sides (Todhunter): + # + # https://en.wikipedia.org/wiki/Spherical_trigonometry#Area_and_spherical_excess + # + # The interior angle of a side is calculated with a special + # case of the Vincenty formula: + # https://en.wikipedia.org/wiki/Great-circle_distance + # ------------------------------------------------------------ + if not great_circle: + raise ValueError( + "Must set great_circle=True to allow the derivation " + "of volume weights from spherical polygons composed " + "from great circle segments." + ) + + if methods: + comp[(axis,)] = "volume spherical polygon geometry" + return True + + if measure: + r = radius + + # actual_volume = + # [actual_area/(4*pi*r**2)] + # * (4/3)*pi*[(r+delta_z)**3 - r**3)] + volumes *= delta_z**3 / (3 * r**2) + delta_z**2 / r + delta_z + else: + raise ValueError( + "X and Y coordinate bounds must both have units " + "equivalent to either metres, for plane polygon, or radians, " + "for spherical polygon, volume calculations. Got " + f"{x.Units!r} and {y.Units!r}." + ) + + comp[(axis,)] = volumes + + weights_axes.add(axis) + + return True + + +def weights_interior_angles(f, lon, lat, interior_ring=None): + """Find the interior angles at spherical polygon vertices. + + The interior angle at a vertex is that formed by two adjacent + sides. Each interior angle is therefore a function of three + vertices - the vertex at which the interior angle is required and + the two adjacent vertices either side of it. + + **Method** + + Bevis, M., Cambareri, G. Computing the area of a spherical polygon + of arbitrary shape. Math Geol 19, 335–346 (1987). + + http://bemlar.ism.ac.jp/zhuang/Refs/Refs/bevis1987mathgeol.pdf + + Abstract: An algorithm for determining the area of a spherical + polygon of arbitrary shape is presented. The kernel of the + problem is to compute the interior angle at each vertex of the + spherical polygon; a well-known relationship between the area of + a spherical polygon and the sum of its interior angles then may + be exploited. The algorithm has been implemented as a FORTRAN + subroutine and a listing is provided. + + .. versionadded:: UGRIDVER + + .. seealso:: `_weights_central_angles` + + :Parameters: + + lon: `Data` + Longitudes. Must have units of radians, which is not + checked. + + lat: `Data` + Latitudes. Must have units of radians, which is not + checked. + + interior_ring: `Data`, optional + The interior ring flags for geometry cells. + + :Returns: + + `Data` + The interior angles of each spherical polygon, in + units of radians. + + """ + from .query import lt + + # P denotes a vertex at which the interior angle is required, A + # denotes the adjacent point clockwise from P, and B denotes the + # adjacent point anticlockwise from P. + lon_A = lon[..., :-2] + lon_P = lon[..., 1:-1] + lon_B = lon[..., 2:] + + lon_A_minus_lon_P = lon_A - lon_P + lon_B_minus_lon_P = lon_B - lon_P + + cos_lat = lat.cos() + cos_lat_A = cos_lat[..., :-2] + cos_lat_P = cos_lat[..., 1:-1] + cos_lat_B = cos_lat[..., 2:] + + sin_lat = lat.sin() + sin_lat_A = sin_lat[..., :-2] + sin_lat_P = sin_lat[..., 1:-1] + sin_lat_B = sin_lat[..., 2:] + + y = lon_A_minus_lon_P.sin() * cos_lat_A + x = sin_lat_A * cos_lat_P - cos_lat_A * sin_lat_P * lon_A_minus_lon_P.cos() + lat_A_primed = f._Data.arctan2(y, x) + + y = lon_B_minus_lon_P.sin() * cos_lat_B + x = sin_lat_B * cos_lat_P - cos_lat_B * sin_lat_P * lon_B_minus_lon_P.cos() + lat_B_primed = f._Data.arctan2(y, x) + + # The CF vertices here are, in general, given in anticlockwise + # order, so we do alpha_P = lat_B_primed - lat_A_primed, rather + # than the "alpha_P = lat_A_primed - lat_B_primed" given in Bevis + # and Cambareri, which assumes clockwise order. + alpha_P = lat_B_primed - lat_A_primed + + if interior_ring is not None: + # However, interior rings *are* given in clockwise order in + # CF, so we need to negate alpha_P in these cases. + alpha_P.where(interior_ring, -alpha_P, inplace=True) + + alpha_P = alpha_P.where(lt(0), alpha_P + 2 * np.pi, alpha_P) + return alpha_P + + +def weights_central_angles(f, lon, lat): + r"""Find the central angles for spherical great circle line segments. + + The central angle of two points on the sphere is the angle + subtended from the centre of the sphere by the two points on its + surface points. It is calculated with a special case of the + Vincenty formula that is accurate for all spherical distances: + + \Delta \sigma =\arctan {\frac {\sqrt {\left(\cos \phi + _{2}\sin(\Delta \lambda )\right)^{2}+\left(\cos \phi _{1}\sin \phi + _{2}-\sin \phi _{1}\cos \phi _{2}\cos(\Delta \lambda + )\right)^{2}}}{\sin \phi _{1}\sin \phi _{2}+\cos \phi _{1}\cos + \phi _{2}\cos(\Delta \lambda )}} + + The quadrant for \Delta \sigma should be governed by the signs of + the numerator and denominator of the right hand side using the + atan2 function. + + (https://en.wikipedia.org/wiki/Great-circle_distance) + + .. versionadded:: UGRIDVER + + :Parameters: + + lon: `Data` + Longitudes. Must have units of radians, which is not + checked. + + lat: `Data` + Latitudes. Must have units of radians, which is not + checked. + + :Returns: + + `Data` + The central angles in units of radians. + + """ + # A and B denote adjacent vertices that define each line segment + delta_lon = lon.diff(axis=-1) + + cos_lat = lat.cos() + sin_lat = lat.sin() + + cos_lat_A = cos_lat[..., :-1] + cos_lat_B = cos_lat[..., 1:] + + sin_lat_A = sin_lat[..., :-1] + sin_lat_B = sin_lat[..., 1:] + + cos_delta_lon = delta_lon.cos() + sin_delta_lon = delta_lon.sin() + + y = ( + (cos_lat_B * sin_delta_lon) ** 2 + + (cos_lat_A * sin_lat_B - sin_lat_A * cos_lat_B * cos_delta_lon) ** 2 + ) ** 0.5 + x = sin_lat_A * sin_lat_B + cos_lat_A * cos_lat_B * cos_delta_lon + + delta_sigma = f._Data.arctan2(y, x) + delta_sigma.override_units(Units("1"), inplace=True) + return delta_sigma + + +def weights_linear( + f, + axis, + comp, + weights_axes, + auto=False, + measure=False, + methods=False, +): + """1-d linear weights from dimension coordinate constructs. + + .. versionadded:: UGRIDVER + + :Parameters: + + axis: `str` + + comp: `dict` + + weights_axes: `set` + + auto: `bool` + + measure: `bool` + If True then make sure that the weights represent true + cell sizes. + + methods: `bool`, optional + If True then add a description of the method used to + create the weights to the *comp* dictionary, as opposed to + the actual weights. + + :Returns: + + `bool` + + """ + da_key = f.domain_axis(axis, key=True, default=None) + if da_key is None: + if auto: + return False + + raise ValueError( + "Can't create weights: Can't find domain axis matching " + f"{axis!r}" + ) + + dim = f.dimension_coordinate(filter_by_axis=(da_key,), default=None) + if dim is None: + if auto: + return False + + raise ValueError( + f"Can't create linear weights for {axis!r} axis: Can't find " + "dimension coordinate construct." + ) + + if not measure and dim.size == 1: + return False + + if da_key in weights_axes: + if auto: + return False + + raise ValueError( + f"Can't create linear weights for {axis!r} axis: Multiple " + "axis specifications" + ) + + if not dim.has_bounds(): + # Dimension coordinate has no bounds + if auto: + return False + + raise ValueError( + f"Can't create linear weights for {axis!r} axis: No bounds" + ) + else: + # Bounds exist + if methods: + comp[(da_key,)] = "linear " + f.constructs.domain_axis_identity( + da_key + ) + else: + comp[(da_key,)] = dim.cellsize + + weights_axes.add(da_key) + + return True + + +def weights_measure(f, measure, comp, weights_axes, methods=False, auto=False): + """Cell measure weights. + + .. versionadded:: UGRIDVER + + :Parameters: + + methods: `bool`, optional + If True then add a description of the method used to + create the weights to the *comp* dictionary, as opposed to + the actual weights. + + :Returns: + + `bool` + + """ + m = f.cell_measures(filter_by_measure=(measure,), todict=True) + len_m = len(m) + + if not len_m: + if measure == "area": + return False + + if auto: + return + + raise ValueError(f"Can't find weights: No {measure!r} cell measure") + + elif len_m > 1: + if auto: + return False + + raise ValueError( + f"Can't find weights: Multiple {measure!r} cell measures" + ) + + key, clm = m.popitem() + + clm_axes0 = f.get_data_axes(key) + + clm_axes = tuple( + [axis for axis, n in zip(clm_axes0, clm.data.shape) if n > 1] + ) + + for axis in clm_axes: + if axis in weights_axes: + if auto: + return False + + raise ValueError( + "Multiple weights specifications for " + f"{f.constructs.domain_axis_identity(axis)!r} axis" + ) + + clm = clm.get_data(_fill_value=False).copy() + if clm_axes != clm_axes0: + iaxes = [clm_axes0.index(axis) for axis in clm_axes] + clm.squeeze(iaxes, inplace=True) + + if methods: + comp[tuple(clm_axes)] = measure + " cell measure" + else: + comp[tuple(clm_axes)] = clm + + weights_axes.update(clm_axes) + + return True + + +def weights_scale(w, scale): + """Scale the weights so that they are <= scale. + + .. versionadded:: UGRIDVER + + :Parameters: + + w: `Data` + The weights to be scaled. + + scale: number or `None` + The maximum value of the scaled weights. If `None` then no + scaling is applied. + + :Returns: + + `Data` + The scaled weights. + + """ + if scale is None: + return w + + if scale <= 0: + raise ValueError( + "Can't set 'scale' parameter to a negative number. Got " + f"{scale!r}" + ) + + w = w / w.max() + if scale != 1: + w = w * scale + + return w + + +def weights_yyy(f, domain_axis, geometry_type, methods=False, auto=False): + """Checks whether weights can be created for given coordinates. + + .. versionadded:: UGRIDVER + + :Parameters: + + domain_axis : `str` or `None` + + geometry_type: `str` + Either ``'polygon'`` or ``'line'``. + + auto: `bool` + + :Returns: + + `tuple` + + """ + aux_X = None + aux_Y = None + aux_Z = None + x_axis = None + y_axis = None + z_axis = None + + auxiliary_coordinates_1d = f.auxiliary_coordinates( + filter_by_naxes=(1,), todict=True + ) + + for key, aux in auxiliary_coordinates_1d.items(): + aux_axis = f.get_data_axes(key)[0] + dt = f.domain_topology(default=None, filter_by_axis=(aux_axis,)) + if dt is None and aux.get_geometry(None) != geometry_type: + continue + + # Still here? Then this auxiliary coordinate has either UGRID + # or geometry cells. + if aux.X: + aux_X = aux.copy() + x_axis = aux_axis + if domain_axis is not None and x_axis != domain_axis: + aux_X = None + continue + elif aux.Y: + aux_Y = aux.copy() + y_axis = aux_axis + if domain_axis is not None and y_axis != domain_axis: + aux_Y = None + continue + elif aux.Z: + aux_Z = aux.copy() + z_axis = aux_axis + if domain_axis is not None and z_axis != domain_axis: + aux_Z = None + continue + + if aux_X is None or aux_Y is None: + if auto: + return (None,) * 4 + + raise ValueError( + "Can't create weights: Need both X and Y nodes to " + f"calculate {geometry_type} cell weights" + ) + + if x_axis != y_axis: + if auto: + return (None,) * 4 + + raise ValueError( + "Can't create weights: X and Y cells span different " "domain axes" + ) + + axis = x_axis + + if aux_X.get_bounds(None) is None or aux_Y.get_bounds(None) is None: + # Not both X and Y coordinates have bounds + if auto: + return (None,) * 4 + + raise ValueError( + "Can't create weights: " "Not both X and Y coordinates have bounds" + ) + + if aux_X.bounds.shape != aux_Y.bounds.shape: + if auto: + return (None,) * 4 + + raise ValueError( + "Can't create weights: UGRID/geometry X and Y coordinate " + "bounds must have the same shape. " + f"Got {aux_X.bounds.shape} and {aux_Y.bounds.shape}" + ) + + if aux_Z is None: + for key, aux in auxiliary_coordinates_1d.items(): + if aux.Z: + aux_Z = aux.copy() + z_axis = f.get_data_axes(key)[0] + + # Check Z coordinates + if aux_Z is not None: + if z_axis != x_axis: + if auto: + return (None,) * 4 + + raise ValueError( + "Z coordinates span different domain axis to X and Y " + "geometry coordinates" + ) + + return axis, aux_X, aux_Y, aux_Z