Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix incorrect result units for some differential operators #808

Merged
merged 3 commits into from
Sep 9, 2024
Merged

Fix incorrect result units for some differential operators #808

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
186acb3
dev
davidhassell Aug 30, 2024
d250a64
ignore_coordinate_units
davidhassell Sep 2, 2024
bdfcd88
ignore_coordinate_units
davidhassell Sep 2, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions Changelog.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
62 changes: 53 additions & 9 deletions cf/field.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
[[-- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- --]
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.]
Expand Down Expand Up @@ -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
)
Expand Down Expand Up @@ -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.

Expand Down Expand Up @@ -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
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Very clearly documented, +1.

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}}
Expand Down Expand Up @@ -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

Expand Down
36 changes: 24 additions & 12 deletions cf/maths.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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)
Expand Down
32 changes: 24 additions & 8 deletions cf/test/test_Field.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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,
)
Expand All @@ -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),
Expand Down
8 changes: 4 additions & 4 deletions cf/test/test_Maths.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading