diff --git a/Changelog.rst b/Changelog.rst index 13098a9eb3..918ba2a6d2 100644 --- a/Changelog.rst +++ b/Changelog.rst @@ -3,6 +3,13 @@ version 3.17.0 **2024-??-??** +* Added the ``cell_measures`` and ``coordinates`` keyword arguments to + `cf.Field.weights` + (https://github.com/NCAS-CMS/cf-python/issues/709) +* Added the ``cell_measures``, ``coordinates``, + ``return_cell_measure``, and ``methods`` keyword arguments to + `cf.Field.cell_area` + (https://github.com/NCAS-CMS/cf-python/issues/709) * Allow `cf.Data` to be initialised with `xarray.DataAarray` (https://github.com/NCAS-CMS/cf-python/issues/706) * Fix bug that caused `cf.Field.del_file_location` to fail when diff --git a/cf/field.py b/cf/field.py index c970460636..19e40b2535 100644 --- a/cf/field.py +++ b/cf/field.py @@ -11,6 +11,7 @@ from . import ( AuxiliaryCoordinate, Bounds, + CellMeasure, CellMethod, Constructs, Count, @@ -2338,11 +2339,14 @@ def cell_area( self, radius="earth", great_circle=False, - set=False, + cell_measures=True, + coordinates=True, + methods=False, + return_cell_measure=False, insert=False, force=False, ): - """Return a field containing horizontal cell areas. + """Return the horizontal cell areas. .. versionadded:: 1.0 @@ -2374,21 +2378,92 @@ def cell_area( .. versionadded:: 3.2.0 + cell_measures: `bool`, optional + If True, the default, then area cell measure + constructs are considered for cell area + creation. Otherwise they are ignored. + + .. versionadded:: 3.17.0 + + coordinates: `bool`, optional + If True, the default, then coordinate constructs are + considered for cell area creation. Otherwise they are + ignored. + + .. versionadded:: 3.17.0 + + methods: `bool`, optional + If True, then return a dictionary describing the method + used to create the cell areas instead of the default, + a field construct. + + .. versionadded:: 3.17.0 + + return_cell_measure: `bool`, optional + If True, then return a cell measure construct instead + of the default, a field construct. + + .. versionadded:: 3.17.0 + insert: deprecated at version 3.0.0 force: deprecated at version 3.0.0 :Returns: - `Field` - A field construct containing the horizontal cell - areas. + `Field`, `CellMeasure`, or `dict` + A field construct, or cell measure construct containing + the horizontal cell areas if *return_cell_measure* is True, + or a dictionary describing the method used to create the + cell areas if *methods* is True. **Examples** + >>> f = cf.example_field(0) + >>> print(f) + Field: specific_humidity (ncvar%q) + ---------------------------------- + Data : specific_humidity(latitude(5), longitude(8)) 1 + Cell methods : area: mean + Dimension coords: latitude(5) = [-75.0, ..., 75.0] degrees_north + : longitude(8) = [22.5, ..., 337.5] degrees_east + : time(1) = [2019-01-01 00:00:00] + >>> a = f.cell_area() + >>> a + + >>> print(a.array) + [[4.27128714e+12 4.27128714e+12 4.27128714e+12 4.27128714e+12 + 4.27128714e+12 4.27128714e+12 4.27128714e+12 4.27128714e+12] + [1.16693735e+13 1.16693735e+13 1.16693735e+13 1.16693735e+13 + 1.16693735e+13 1.16693735e+13 1.16693735e+13 1.16693735e+13] + [3.18813213e+13 3.18813213e+13 3.18813213e+13 3.18813213e+13 + 3.18813213e+13 3.18813213e+13 3.18813213e+13 3.18813213e+13] + [1.16693735e+13 1.16693735e+13 1.16693735e+13 1.16693735e+13 + 1.16693735e+13 1.16693735e+13 1.16693735e+13 1.16693735e+13] + [4.27128714e+12 4.27128714e+12 4.27128714e+12 4.27128714e+12 + 4.27128714e+12 4.27128714e+12 4.27128714e+12 4.27128714e+12]] + >>> f.cell_area(methods=True) + {(1,): 'linear longitude', (0,): 'linear sine latitude'} + >>> a = f.cell_area(radius=cf.Data(3389.5, 'km')) - >>> a = f.cell_area(insert=True) + + >>> c = f.cell_area(return_cell_measure=True) + >>> c + + >>> f.set_construct(c) + 'cellmeasure0' + >>> print(f) + Field: specific_humidity (ncvar%q) + ---------------------------------- + Data : specific_humidity(latitude(5), longitude(8)) 1 + Cell methods : area: mean + Dimension coords: latitude(5) = [-75.0, ..., 75.0] degrees_north + : longitude(8) = [22.5, ..., 337.5] degrees_east + : time(1) = [2019-01-01 00:00:00] + Cell measures : measure:area(latitude(5), longitude(8)) = [[4271287143027.272, ..., 4271287143027.272]] m2 + >>> f.cell_area(methods=True) + {(0, 1): 'area cell measure'} """ if insert: @@ -2415,9 +2490,24 @@ def cell_area( measure=True, scale=None, great_circle=great_circle, + cell_measures=cell_measures, + coordinates=coordinates, + methods=methods, ) - w.set_property("standard_name", "cell_area", copy=False) + if methods: + if return_cell_measure: + raise ValueError( + "Can't set both the 'methods' and 'return_cell_measure'" + "parameters." + ) + return w + + if return_cell_measure: + w = CellMeasure(source=w, copy=False) + w.set_measure("area") + else: + w.set_property("standard_name", "cell_area", copy=False) return w @@ -3128,6 +3218,8 @@ def weights( data=False, great_circle=False, axes=None, + cell_measures=True, + coordinates=True, **kwargs, ): """Return weights for the data array values. @@ -3194,19 +3286,33 @@ def weights( following methods, in order of preference, - 1. Volume cell measures + If the *cell_measures* parameter is + True: + + 1. Volume cell measures (see the note + on the *measure* parameter) + 2. Area cell measures + + If the *coordinates* parameter is True: + 3. Area calculated from X and Y dimension coordinate constructs with bounds + 4. Area calculated from 1-d auxiliary - coordinate constructs for geometries - or a UGRID mesh topology. - 5. Length calculated from 1-d auxiliary - coordinate constructs for geometries - or a UGRID mesh topology. + coordinate constructs for + geometries or a UGRID mesh + topology. + + 5. Length calculated from 1-d + auxiliary coordinate constructs for + geometries or a UGRID mesh + topology. + 6. Cell sizes of dimension coordinate constructs with bounds + 7. Equal weights and the outer product of these weights @@ -3245,20 +3351,31 @@ def weights( field by the following methods, in order of preference, - 1. Area cell measures. - 2. X and Y dimension coordinate - constructs with bounds. - 3. X and Y 1-d auxiliary coordinate - constructs for polygon cells - defined by geometries or a UGRID - mesh topology. + If the *cell_measures* parameter is + True: + + 1. Area cell measures. + + If the *coordinates* parameter is + True: + + 2. X and Y dimension coordinate + constructs with bounds. + + 3. X and Y 1-d auxiliary coordinate + constructs for polygon cells + defined by geometries or a UGRID + mesh topology. Set the *methods* parameter to find out how the weights were actually created. ``'volume'`` Cell volume weights from the field - construct's volume cell measure construct. + construct's volume cell measure + construct (see the note on the + *measure* parameter). Requires the + *cell_measures* parameter to be True. `str` Weights from the cell sizes of the dimension coordinate construct with this @@ -3381,6 +3498,23 @@ def weights( .. versionadded:: 3.3.0 + cell_measures: `bool`, optional + If True, the default, then area and volume cell + measure constructs are considered for weights creation + when *weights* is `True`, ``'area'``, or + ``'volume'``. If False then cell measure constructs + are ignored for these *weights*. + + .. versionadded:: 3.17.0 + + coordinates: `bool`, optional + If True, the default, then coordinate constructs are + considered for weights creation for *weights* of + `True` or ``'area'``. If False then coordinate + constructs are ignored for these *weights*. + + .. versionadded:: 3.17.0 + kwargs: deprecated at version 3.0.0. :Returns: @@ -3476,18 +3610,23 @@ def weights( # Auto-detect all weights # -------------------------------------------------------- # Volume weights - if Weights.cell_measure( + if cell_measures and Weights.cell_measure( self, "volume", comp, weights_axes, methods=methods, auto=True ): # Found volume weights from cell measures pass - elif Weights.cell_measure( - self, "area", comp, weights_axes, methods=methods, auto=True + elif cell_measures and Weights.cell_measure( + self, + "area", + comp, + weights_axes, + methods=methods, + auto=True, ): # Found area weights from cell measures pass - elif Weights.area_XY( + elif coordinates and Weights.area_XY( self, comp, weights_axes, @@ -3502,44 +3641,45 @@ def weights( domain_axes = self.domain_axes(todict=True) - for da_key in domain_axes: - if Weights.polygon_area( - self, - da_key, - comp, - weights_axes, - measure=measure, - radius=radius, - great_circle=great_circle, - methods=methods, - auto=True, - ): - # Found area weights from polygon geometries - pass - elif Weights.line_length( - self, - da_key, - comp, - weights_axes, - measure=measure, - radius=radius, - great_circle=great_circle, - methods=methods, - auto=True, - ): - # Found linear weights from line geometries - pass - elif Weights.linear( - self, - da_key, - comp, - weights_axes, - measure=measure, - methods=methods, - auto=True, - ): - # Found linear weights from dimension coordinates - pass + if coordinates: + for da_key in domain_axes: + if Weights.polygon_area( + self, + da_key, + comp, + weights_axes, + measure=measure, + radius=radius, + great_circle=great_circle, + methods=methods, + auto=True, + ): + # Found area weights from polygon geometries + pass + elif Weights.line_length( + self, + da_key, + comp, + weights_axes, + measure=measure, + radius=radius, + great_circle=great_circle, + methods=methods, + auto=True, + ): + # Found linear weights from line geometries + pass + elif Weights.linear( + self, + da_key, + comp, + weights_axes, + measure=measure, + methods=methods, + auto=True, + ): + # Found linear weights from dimension coordinates + pass weights_axes = [] for key in comp: @@ -3615,11 +3755,11 @@ def weights( # -------------------------------------------------------- fields = [] axes = [] - cell_measures = [] + measures = [] if isinstance(weights, str): if weights in ("area", "volume"): - cell_measures = (weights,) + measures = (weights,) else: axes.append(weights) else: @@ -3631,7 +3771,7 @@ def weights( weights = iter(weights) except TypeError: raise TypeError( - f"Invalid type of 'weights' parameter: {weights}" + f"Invalid type of 'weights' parameter: {weights!r}" ) for w in tuple(weights): @@ -3639,10 +3779,10 @@ def weights( fields.append(w) elif isinstance(w, Data): raise ValueError( - f"Invalid weight {w} in sequence of weights." + f"Invalid weight {w!r} in sequence of weights." ) elif w in ("area", "volume"): - cell_measures.append(w) + measures.append(w) else: axes.append(w) @@ -3650,7 +3790,13 @@ def weights( Weights.field(self, fields, comp, weights_axes) # Volume weights - if "volume" in cell_measures: + if "volume" in measures: + if not cell_measures: + raise ValueError( + "Can't create weights: Unable to use " + "volume cell measures when cell_meaures=False" + ) + Weights.cell_measure( self, "volume", @@ -3661,41 +3807,51 @@ def weights( ) # Area weights - if "area" in cell_measures: - if Weights.cell_measure( - self, - "area", - comp, - weights_axes, - methods=methods, - auto=True, - ): - # Found area weights from cell measures - pass - elif Weights.area_XY( - self, - comp, - weights_axes, - measure=measure, - radius=radius, - methods=methods, - auto=True, - ): - # Found area weights from X and Y dimension - # coordinates - pass - else: - # Found area weights from UGRID/geometry cells - Weights.polygon_area( + if "area" in measures: + area_weights = False + if cell_measures: + if Weights.cell_measure( + self, + "area", + comp, + weights_axes, + methods=methods, + auto=True, + ): + # Found area weights from cell measures + area_weights = True + + if not area_weights and coordinates: + if Weights.area_XY( self, - None, comp, weights_axes, measure=measure, radius=radius, - great_circle=great_circle, methods=methods, - auto=False, + auto=True, + ): + # Found area weights from X and Y dimension + # coordinates + area_weights = True + else: + Weights.polygon_area( + self, + None, + comp, + weights_axes, + measure=measure, + radius=radius, + great_circle=great_circle, + methods=methods, + auto=False, + ) + # Found area weights from UGRID/geometry cells + area_weights = True + + if not area_weights: + raise ValueError( + "Can't create weights: Unable to find cell areas" ) for axis in axes: @@ -3743,26 +3899,26 @@ def weights( auto=False, ) - # Check for area weights specified by X and Y axes - # separately and replace them with area weights - xaxis = self.domain_axis("X", key=True, default=None) - yaxis = self.domain_axis("Y", key=True, default=None) - if xaxis != yaxis and (xaxis,) in comp and (yaxis,) in comp: - del comp[(xaxis,)] - del comp[(yaxis,)] - weights_axes.discard(xaxis) - weights_axes.discard(yaxis) - if not Weights.cell_measure( - self, "area", comp, weights_axes, methods=methods - ): - Weights.area_XY( - self, - comp, - weights_axes, - measure=measure, - radius=radius, - methods=methods, - ) + # Check for area weights specified by X and Y axes + # separately and replace them with area weights + xaxis = self.domain_axis("X", key=True, default=None) + yaxis = self.domain_axis("Y", key=True, default=None) + if xaxis != yaxis and (xaxis,) in comp and (yaxis,) in comp: + del comp[(xaxis,)] + del comp[(yaxis,)] + weights_axes.discard(xaxis) + weights_axes.discard(yaxis) + if not Weights.cell_measure( + self, "area", comp, weights_axes, methods=methods + ): + Weights.area_XY( + self, + comp, + weights_axes, + measure=measure, + radius=radius, + methods=methods, + ) if not methods: if scale is not None: diff --git a/cf/test/test_Field.py b/cf/test/test_Field.py index 6927f9813b..b253484d72 100644 --- a/cf/test/test_Field.py +++ b/cf/test/test_Field.py @@ -900,6 +900,13 @@ def test_Field_cell_area(self): y = f.dimension_coordinate("Y") self.assertTrue(y.has_bounds()) + ca = f.cell_area(return_cell_measure=True) + self.assertIsInstance(ca, cf.CellMeasure) + self.assertEqual(ca.get_measure(), "area") + + m = f.cell_area(methods=True) + self.assertIsInstance(m, dict) + def test_Field_radius(self): f = self.f.copy() diff --git a/cf/test/test_weights.py b/cf/test/test_weights.py index 7f3bbb5cdf..e72c3399da 100644 --- a/cf/test/test_weights.py +++ b/cf/test/test_weights.py @@ -291,6 +291,41 @@ def test_weights_line_area_ugrid(self): ) self.assertEqual(w.Units, cf.Units("m")) + def test_weights_cell_measures_coordinates(self): + f = cf.example_field(0) + + areas1 = f.cell_area() + areas2 = areas1.copy() + areas2[...] = -1 + + c = cf.CellMeasure(source=areas2) + c.set_measure("area") + f.set_construct(c) + + w = f.weights(True, measure=True) + self.assertTrue(w.data.allclose(areas2.data)) + + w = f.weights(True, measure=True, cell_measures=False) + self.assertTrue(w.data.allclose(areas1.data)) + + w = f.weights(True, measure=True, coordinates=False) + self.assertTrue(w.data.allclose(areas2.data)) + + with self.assertRaises(ValueError): + w = f.weights(True, cell_measures=False, coordinates=False) + + w = f.weights("area", measure=True) + self.assertTrue(w.data.allclose(areas2.data)) + + w = f.weights("area", measure=True, cell_measures=False) + self.assertTrue(w.data.allclose(areas1.data)) + + w = f.weights("area", measure=True, coordinates=False) + self.assertTrue(w.data.allclose(areas2.data)) + + with self.assertRaises(ValueError): + w = f.weights("area", cell_measures=False, coordinates=False) + if __name__ == "__main__": print("Run date:", datetime.datetime.now()) diff --git a/cf/weights.py b/cf/weights.py index 4a25c749a1..023854865a 100644 --- a/cf/weights.py +++ b/cf/weights.py @@ -138,12 +138,11 @@ def area_XY( if ycoord.Units.equivalent(radians): ycoord = ycoord.clip(-90, 90, units=Units("degrees")) - ycoord.sin(inplace=True) - + ysin = ycoord.sin() if methods: weights[(yaxis,)] = f"linear sine {ycoord.identity()}" else: - cells = ycoord.cellsize + cells = ysin.cellsize if measure: cells = cells * radius @@ -1522,6 +1521,16 @@ def cell_measure( key, clm = m.popitem() + if not clm.has_data(): + if auto: + return False + + raise ValueError( + f"Can't find weights: Cell measure {m!r} has no data, " + "possibly because it is external. " + "Consider setting cell_measures=False" + ) + clm_axes0 = f.get_data_axes(key) clm_axes = tuple(