From 186acb3187771512c5db768891eaf9eb5fcdaaeb Mon Sep 17 00:00:00 2001 From: David Hassell Date: Fri, 30 Aug 2024 16:39:11 +0100 Subject: [PATCH 1/3] dev --- cf/field.py | 65 +++++++++++++++++++++++++++++++++++++------ cf/maths.py | 36 ++++++++++++++++-------- cf/test/test_Field.py | 32 +++++++++++++++------ cf/test/test_Maths.py | 8 +++--- 4 files changed, 108 insertions(+), 33 deletions(-) diff --git a/cf/field.py b/cf/field.py index 71bc209dc4..209319d83f 100644 --- a/cf/field.py +++ b/cf/field.py @@ -2982,7 +2982,7 @@ def laplacian_xy( [0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1]] >>> lp = f.laplacian_xy(radius='earth') >>> lp - + >>> print(lp.array) [[-- -- -- -- -- -- -- --] [-- -- -- -- -- -- -- --] @@ -3050,19 +3050,31 @@ def laplacian_xy( r2_sin_theta = sin_theta * r**2 d2f_dphi2 = f.derivative( - x_key, wrap=x_wrap, one_sided_at_boundary=one_sided_at_boundary + x_key, + wrap=x_wrap, + one_sided_at_boundary=one_sided_at_boundary, + ignore_coordinate_units=True, ).derivative( - x_key, wrap=x_wrap, one_sided_at_boundary=one_sided_at_boundary + x_key, + wrap=x_wrap, + one_sided_at_boundary=one_sided_at_boundary, + ignore_coordinate_units=True, ) term1 = d2f_dphi2 / (r2_sin_theta * sin_theta) df_dtheta = f.derivative( - y_key, wrap=None, one_sided_at_boundary=one_sided_at_boundary + y_key, + wrap=None, + one_sided_at_boundary=one_sided_at_boundary, + ignore_coordinate_units=True, ) term2 = (df_dtheta * sin_theta).derivative( - y_key, wrap=None, one_sided_at_boundary=one_sided_at_boundary + y_key, + wrap=None, + one_sided_at_boundary=one_sided_at_boundary, + ignore_coordinate_units=True, ) / r2_sin_theta f = term1 + term2 @@ -11607,8 +11619,8 @@ def grad_xy(self, x_wrap=None, one_sided_at_boundary=False, radius=None): [0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1]] >>> fx, fy = f.grad_xy(radius='earth') >>> fx, fy - (, - ) + (, + ) >>> print(fx.array) [[0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0.] @@ -11679,7 +11691,10 @@ def grad_xy(self, x_wrap=None, one_sided_at_boundary=False, radius=None): r = f.radius(default=radius) X = f.derivative( - x_key, wrap=x_wrap, one_sided_at_boundary=one_sided_at_boundary + x_key, + wrap=x_wrap, + one_sided_at_boundary=one_sided_at_boundary, + ignore_coordinate_units=True, ) / (theta.sin() * r) Y = ( @@ -11687,6 +11702,7 @@ def grad_xy(self, x_wrap=None, one_sided_at_boundary=False, radius=None): y_key, wrap=None, one_sided_at_boundary=one_sided_at_boundary, + ignore_coordinate_units=True, ) / r ) @@ -14226,9 +14242,10 @@ def derivative( axis, wrap=None, one_sided_at_boundary=False, + cyclic=None, + ignore_coordinate_units=False, inplace=False, i=False, - cyclic=None, ): """Calculate the derivative along the specified axis. @@ -14262,6 +14279,31 @@ def derivative( at the non-cyclic boundaries. By default missing values are set at non-cyclic boundaries. + ignore_coordinate_units: `bool`, optional + If True then the coordinates providing the cell + spacings along the axis are assumed to dimensionless, + even if they do in fact have units. This does not + change the returned numerical values, but will alter + the units of the returned field. If False (the + default) then the coordinate units will propagate + through to the result. + + If *ignore_coordinate_units* is False then the units + of the returned field construct will be the original + units divided by the coordinate units. + + If *ignore_coordinate_units* is True then the units of + the returned field construct will be identical to the + original units. + + For example, for a field construct with units of + ``m.s-1`` and X coordinate units of ``radians``, the + units of the X derivative will be ``m.s-1.radians-1`` + by default, or ``m.s-1`` if *ignore_coordinate_units* + is True. + + .. versionadded:: NEXTVERSION + {{inplace: `bool`, optional}} {{i: deprecated at version 3.0.0}} @@ -14393,6 +14435,11 @@ def derivative( for _ in range(self.ndim - 1 - axis_index): d.insert_dimension(position=1, inplace=True) + if ignore_coordinate_units: + # Remove the coordinate units before we calculate the + # derivative + d.override_units(None, inplace=True) + # Find the derivative f.data /= d diff --git a/cf/maths.py b/cf/maths.py index 52b4a4f1b0..b1cfb33646 100644 --- a/cf/maths.py +++ b/cf/maths.py @@ -372,11 +372,11 @@ def curl_xy(fx, fy, x_wrap=None, one_sided_at_boundary=False, radius=None): [0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1]] >>> fx, fy = f.grad_xy(radius='earth', one_sided_at_boundary=True) >>> fx, fy - (, - ) + (, + ) >>> c = cf.curl_xy(fx, fy, radius='earth') >>> c - + >>> print(c.array) [[-- -- -- -- -- -- -- --] [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0] @@ -437,8 +437,7 @@ def curl_xy(fx, fy, x_wrap=None, one_sided_at_boundary=False, radius=None): # -------------------------------------------------------- # Spherical polar coordinates # -------------------------------------------------------- - # Convert latitude and longitude units to radians, so that the - # units of the result are nice. + # Convert latitude and longitude units to radians radians = Units("radians") fx_x_coord.Units = radians fx_y_coord.Units = radians @@ -461,10 +460,16 @@ def curl_xy(fx, fy, x_wrap=None, one_sided_at_boundary=False, radius=None): r = fx.radius(default=radius) term1 = (fx * sin_theta).derivative( - fx_y_key, wrap=None, one_sided_at_boundary=one_sided_at_boundary + fx_y_key, + wrap=None, + one_sided_at_boundary=one_sided_at_boundary, + ignore_coordinate_units=True, ) term2 = fy.derivative( - fy_x_key, wrap=x_wrap, one_sided_at_boundary=one_sided_at_boundary + fy_x_key, + wrap=x_wrap, + one_sided_at_boundary=one_sided_at_boundary, + ignore_coordinate_units=True, ) c = (term1 - term2) / (sin_theta * r) @@ -597,11 +602,12 @@ def div_xy(fx, fy, x_wrap=None, one_sided_at_boundary=False, radius=None): [0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1]] >>> fx, fy = f.grad_xy(radius='earth', one_sided_at_boundary=True) >>> fx, fy - (, - ) + (, + ) >>> d = cf.div_xy(fx, fy, radius='earth') >>> d - + + >>> print(d.array) [[-- -- -- -- -- -- -- --] [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0] @@ -686,11 +692,17 @@ def div_xy(fx, fy, x_wrap=None, one_sided_at_boundary=False, radius=None): r = fx.radius(default=radius) term1 = fx.derivative( - fx_x_key, wrap=x_wrap, one_sided_at_boundary=one_sided_at_boundary + fx_x_key, + wrap=x_wrap, + one_sided_at_boundary=one_sided_at_boundary, + ignore_coordinate_units=True, ) term2 = (fy * sin_theta).derivative( - fy_y_key, wrap=None, one_sided_at_boundary=one_sided_at_boundary + fy_y_key, + wrap=None, + one_sided_at_boundary=one_sided_at_boundary, + ignore_coordinate_units=True, ) d = (term1 + term2) / (sin_theta * r) diff --git a/cf/test/test_Field.py b/cf/test/test_Field.py index 13366c52f2..064780852f 100644 --- a/cf/test/test_Field.py +++ b/cf/test/test_Field.py @@ -2018,9 +2018,15 @@ def test_Field_moving_window(self): def test_Field_derivative(self): f = cf.example_field(0) f[...] = np.arange(9)[1:] * 45 + x = f.dimension_coordinate("X") + + # Ignore coordinate units + d = f.derivative("X", ignore_coordinate_units=True) + self.assertEqual(d.Units, f.Units) # Check a cyclic periodic axis d = f.derivative("X") + self.assertEqual(d.Units, f.Units / x.Units) self.assertTrue(np.allclose(d[:, 1:-1].array, 1)) self.assertTrue(np.allclose(d[:, [0, -1]].array, -3)) @@ -2028,7 +2034,7 @@ def test_Field_derivative(self): # case f1 = f[:, ::-1] d1 = f1.derivative("X") - self.assertTrue(d1.data.equals(d.data)) + self.assertTrue(d1.data.equals(d.data, verbose=-1)) # Check non-cyclic d = f.derivative("X", wrap=False) @@ -2490,12 +2496,21 @@ def test_Field_grad_xy(self): radius=radius, x_wrap=wrap, one_sided_at_boundary=one_sided ) - self.assertTrue(x.Units == y.Units == cf.Units("m-1 rad-1")) + self.assertEqual(x.Units, y.Units) + self.assertEqual(y.Units, cf.Units("m-1")) x0 = f.derivative( - "X", wrap=wrap, one_sided_at_boundary=one_sided + "X", + wrap=wrap, + one_sided_at_boundary=one_sided, ) / (sin_theta * r) - y0 = f.derivative("Y", one_sided_at_boundary=one_sided) / r + y0 = ( + f.derivative( + "Y", + one_sided_at_boundary=one_sided, + ) + / r + ) # Check the data with cf.rtol(1e-10): @@ -2528,7 +2543,8 @@ def test_Field_grad_xy(self): for one_sided in (True, False): x, y = f.grad_xy(x_wrap=wrap, one_sided_at_boundary=one_sided) - self.assertTrue(x.Units == y.Units == cf.Units("m-1")) + self.assertEqual(x.Units, y.Units) + self.assertEqual(y.Units, cf.Units("m-1")) x0 = f.derivative( "X", wrap=wrap, one_sided_at_boundary=one_sided @@ -2572,7 +2588,7 @@ def test_Field_laplacian_xy(self): radius=radius, x_wrap=wrap, one_sided_at_boundary=one_sided ) - self.assertTrue(lp.Units == cf.Units("m-2 rad-2")) + self.assertEqual(lp.Units, cf.Units("m-2")) lp0 = cf.div_xy( *f.grad_xy( @@ -2580,7 +2596,7 @@ def test_Field_laplacian_xy(self): x_wrap=wrap, one_sided_at_boundary=one_sided, ), - radius=2, + radius=radius, x_wrap=wrap, one_sided_at_boundary=one_sided, ) @@ -2604,7 +2620,7 @@ def test_Field_laplacian_xy(self): x_wrap=wrap, one_sided_at_boundary=one_sided ) - self.assertTrue(lp.Units == cf.Units("m-2")) + self.assertEqual(lp.Units, cf.Units("m-2")) lp0 = cf.div_xy( *f.grad_xy(x_wrap=wrap, one_sided_at_boundary=one_sided), diff --git a/cf/test/test_Maths.py b/cf/test/test_Maths.py index 0900692aeb..349bd495dd 100644 --- a/cf/test/test_Maths.py +++ b/cf/test/test_Maths.py @@ -32,7 +32,7 @@ def test_curl_xy(self): one_sided_at_boundary=one_sided, ) - self.assertTrue(c.Units == cf.Units("m-2 rad-2")) + self.assertEqual(c.Units, cf.Units("m-2")) term1 = (x * sin_theta).derivative( "Y", one_sided_at_boundary=one_sided @@ -68,7 +68,7 @@ def test_curl_xy(self): x, y, x_wrap=wrap, one_sided_at_boundary=one_sided ) - self.assertTrue(d.Units == cf.Units("m-2")) + self.assertEqual(d.Units, cf.Units("m-2")) term1 = x.derivative( "X", wrap=wrap, one_sided_at_boundary=one_sided @@ -121,7 +121,7 @@ def test_div_xy(self): one_sided_at_boundary=one_sided, ) - self.assertTrue(d.Units == cf.Units("m-2 rad-2"), d.Units) + self.assertEqual(d.Units, cf.Units("m-2")) term1 = x.derivative( "X", wrap=wrap, one_sided_at_boundary=one_sided @@ -157,7 +157,7 @@ def test_div_xy(self): x, y, x_wrap=wrap, one_sided_at_boundary=one_sided ) - self.assertTrue(d.Units == cf.Units("m-2")) + self.assertEqual(d.Units, cf.Units("m-2")) term1 = x.derivative( "X", wrap=wrap, one_sided_at_boundary=one_sided From d250a648e9250608859b27f002dc4bd7ba7599a0 Mon Sep 17 00:00:00 2001 From: David Hassell Date: Mon, 2 Sep 2024 14:25:55 +0100 Subject: [PATCH 2/3] ignore_coordinate_units --- cf/field.py | 27 ++++++++++++--------------- 1 file changed, 12 insertions(+), 15 deletions(-) diff --git a/cf/field.py b/cf/field.py index 209319d83f..f10362b2f8 100644 --- a/cf/field.py +++ b/cf/field.py @@ -14281,21 +14281,18 @@ def derivative( ignore_coordinate_units: `bool`, optional If True then the coordinates providing the cell - spacings along the axis are assumed to dimensionless, - even if they do in fact have units. This does not - change the returned numerical values, but will alter - the units of the returned field. If False (the - default) then the coordinate units will propagate - through to the result. - - If *ignore_coordinate_units* is False then the units - of the returned field construct will be the original - units divided by the coordinate units. - - If *ignore_coordinate_units* is True then the units of - the returned field construct will be identical to the + spacings along the specified axis are assumed to be + dimensionless, even if they do in fact have + units. This does not change the magnitude of the + returned numerical values, but the units of the + returned field construct will be identical to the original units. + If False (the default) then the coordinate units will + propagate through to the result. i.e. the units of the + returned field construct will be the original units + divided by the coordinate units. + For example, for a field construct with units of ``m.s-1`` and X coordinate units of ``radians``, the units of the X derivative will be ``m.s-1.radians-1`` @@ -14436,8 +14433,8 @@ def derivative( d.insert_dimension(position=1, inplace=True) if ignore_coordinate_units: - # Remove the coordinate units before we calculate the - # derivative + # Remove the coordinate units from the coordinate + # differences, before we calculate the derivative. d.override_units(None, inplace=True) # Find the derivative From bdfcd88047971571365fb2d4d1ae7943b14024ec Mon Sep 17 00:00:00 2001 From: David Hassell Date: Mon, 2 Sep 2024 14:41:21 +0100 Subject: [PATCH 3/3] ignore_coordinate_units --- Changelog.rst | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/Changelog.rst b/Changelog.rst index dc36568b30..df22d3d73a 100644 --- a/Changelog.rst +++ b/Changelog.rst @@ -7,6 +7,13 @@ version NEXTVERSION (https://github.com/NCAS-CMS/cf-python/issues/784) * Include the UM version as a field property when reading UM files (https://github.com/NCAS-CMS/cf-python/issues/777) +* New keyword parameter to `cf.Field.derivative`: + ``ignore_coordinate_units`` + (https://github.com/NCAS-CMS/cf-python/issues/807) +* Fix bug that sometimes puts an incorrect ``radian-1`` or + ``radian-2`` in the returned units of the differential operator + methods and functions + (https://github.com/NCAS-CMS/cf-python/issues/807) * Fix bug where `cf.example_fields` returned a `list` of Fields rather than a `Fieldlist` (https://github.com/NCAS-CMS/cf-python/issues/725)