Merge pull request #808 from davidhassell/curl

davidhassell · web-flow · commit 807f7409b00f · 2024-09-09T14:24:08.000+01:00
Fix incorrect result units for some differential operators
diff --git 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)
diff --git 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
-        <CF Field: long_name=X-Y Laplacian of specific_humidity(latitude(5), longitude(8)) m-2.rad-2>
+        <CF Field: long_name=X-Y Laplacian of specific_humidity(latitude(5), longitude(8)) m-2>
         >>> 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
-        (<CF Field: long_name=X gradient of specific_humidity(latitude(5), longitude(8)) m-1.rad-1>,
-         <CF Field: long_name=Y gradient of specific_humidity(latitude(5), longitude(8)) m-1.rad-1>)
+        (<CF Field: long_name=X gradient of specific_humidity(latitude(5), longitude(8)) m-1>,
+         <CF Field: long_name=Y gradient of specific_humidity(latitude(5), longitude(8)) m-1>)
         >>> print(fx.array)
         [[0. 0. 0. 0. 0. 0. 0. 0.]
          [0. 0. 0. 0. 0. 0. 0. 0.]
@@ -11679,14 +11691,18 @@ 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 = (
                 f.derivative(
                     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,28 @@ 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 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``
+                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 +14432,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 from the coordinate
+            # differences, 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
@@ -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
-    (<CF Field: long_name=X gradient of specific_humidity(latitude(5), longitude(8)) m-1.rad-1>,
-     <CF Field: long_name=Y gradient of specific_humidity(latitude(5), longitude(8)) m-1.rad-1>)
+    (<CF Field: long_name=X gradient of specific_humidity(latitude(5), longitude(8)) m-1>,
+     <CF Field: long_name=Y gradient of specific_humidity(latitude(5), longitude(8)) m-1>)
     >>> c = cf.curl_xy(fx, fy, radius='earth')
     >>> c
-    <CF Field: long_name=Divergence of (long_name=X gradient of specific_humidity, long_name=Y gradient of specific_humidity)(latitude(5), longitude(8)) m-2.rad-2>
+    <CF Field: long_name=Horizontal curl of (long_name=X gradient of specific_humidity, long_name=Y gradient of specific_humidity)(latitude(5), longitude(8)) m-2>
     >>> 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
-    (<CF Field: long_name=X gradient of specific_humidity(latitude(5), longitude(8)) m-1.rad-1>,
-     <CF Field: long_name=Y gradient of specific_humidity(latitude(5), longitude(8)) m-1.rad-1>)
+    (<CF Field: long_name=X gradient of specific_humidity(latitude(5), longitude(8)) m-1>,
+     <CF Field: long_name=Y gradient of specific_humidity(latitude(5), longitude(8)) m-1>)
     >>> d = cf.div_xy(fx, fy, radius='earth')
     >>> d
-    <CF Field: long_name=Divergence of (long_name=X gradient of specific_humidity, long_name=Y gradient of specific_humidity)(latitude(5), longitude(8)) m-2.rad-2>
+    <CF Field: long_name=Horizontal divergence of (long_name=X gradient of specific_humidity, long_name=Y gradient of specific_humidity)(latitude(5), longitude(8)) m-2>
+
     >>> 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
@@ -2018,17 +2018,23 @@ 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))
 
         # The reversed field should contain the same gradients in this
         # 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,15 +2588,15 @@ 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(
                         radius=radius,
                         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
@@ -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
