Skip to content

Commit 0562abc

Browse files
davidhasselldavidhassell
committed
update Field.lapalcian_xy and Field.grad_xy
1 parent 14ebabd commit 0562abc

File tree

3 files changed

+57
-15
lines changed

3 files changed

+57
-15
lines changed

‎cf/field.py‎

Lines changed: 20 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -3951,7 +3951,8 @@ def laplacian_xy(
39513951

39523952
.. versionadded:: 3.12.0
39533953

3954-
.. seealso:: `derivative`, `grad_xy`, `iscyclic`, `cf.div_xy`
3954+
.. seealso:: `derivative`, `grad_xy`, `iscyclic`,
3955+
`cf.curl_xy`, `cf.div_xy`
39553956

39563957
:Parameters:
39573958

@@ -4050,6 +4051,12 @@ def laplacian_xy(
40504051
x_coord.Units = _units_radians
40514052
y_coord.Units = _units_radians
40524053

4054+
# Ensure that the lat and lon dimension coordinates have
4055+
# standard names, so that metadata-aware broadcasting
4056+
# works as expected when all of their units are radians.
4057+
x_coord.standard_name = "longitude"
4058+
y_coord.standard_name = "latitude"
4059+
40534060
# Get theta as a field that will broadcast to f, and
40544061
# adjust its values so that theta=0 is at the north pole.
40554062
theta = np.pi / 2 - f.convert(y_key, full_domain=True)
@@ -4078,8 +4085,8 @@ def laplacian_xy(
40784085
f = term1 + term2
40794086

40804087
# Reset latitude and longitude coordinate units
4081-
f.dimension_coordinate("X").Units = x_units
4082-
f.dimension_coordinate("Y").Units = y_units
4088+
f.dimension_coordinate("longitude").Units = x_units
4089+
f.dimension_coordinate("latitude").Units = y_units
40834090
else:
40844091
# --------------------------------------------------------
40854092
# Cartesian coordinates
@@ -12392,6 +12399,12 @@ def grad_xy(self, x_wrap=None, one_sided_at_boundary=False, radius=None):
1239212399
x_coord.Units = _units_radians
1239312400
y_coord.Units = _units_radians
1239412401

12402+
# Ensure that the lat and lon dimension coordinates have
12403+
# standard names, so that metadata-aware broadcasting
12404+
# works as expected when all of their units are radians.
12405+
x_coord.standard_name = "longitude"
12406+
y_coord.standard_name = "latitude"
12407+
1239512408
# Get theta as a field that will broadcast to f, and
1239612409
# adjust its values so that theta=0 is at the north pole.
1239712410
theta = np.pi / 2 - f.convert(y_key, full_domain=True)
@@ -12412,11 +12425,11 @@ def grad_xy(self, x_wrap=None, one_sided_at_boundary=False, radius=None):
1241212425
)
1241312426

1241412427
# Reset latitude and longitude coordinate units
12415-
X.dimension_coordinate("X").Units = x_units
12416-
X.dimension_coordinate("Y").Units = y_units
12428+
X.dimension_coordinate("longitude").Units = x_units
12429+
X.dimension_coordinate("latitude").Units = y_units
1241712430

12418-
Y.dimension_coordinate("X").Units = x_units
12419-
Y.dimension_coordinate("Y").Units = y_units
12431+
Y.dimension_coordinate("longitude").Units = x_units
12432+
Y.dimension_coordinate("latitude").Units = y_units
1242012433
else:
1242112434
# --------------------------------------------------------
1242212435
# Cartesian coordinates

‎cf/test/test_Field.py‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2382,6 +2382,23 @@ def test_Field_grad_xy(self):
23822382
self.assertTrue(x.equals(x0, rtol=1e-10))
23832383
self.assertTrue(y.equals(y0, rtol=1e-10))
23842384

2385+
# Test case when spherical dimension coordinates have units
2386+
# but no standard names
2387+
f = cf.example_field(0)
2388+
del f.dimension_coordinate("X").standard_name
2389+
del f.dimension_coordinate("Y").standard_name
2390+
x, y = f.grad_xy(radius="earth")
2391+
self.assertEqual(x.shape, f.shape)
2392+
self.assertEqual(y.shape, f.shape)
2393+
self.assertEqual(x.dimension_coordinate("Y").standard_name, "latitude")
2394+
self.assertEqual(
2395+
x.dimension_coordinate("X").standard_name, "longitude"
2396+
)
2397+
self.assertEqual(y.dimension_coordinate("Y").standard_name, "latitude")
2398+
self.assertEqual(
2399+
y.dimension_coordinate("X").standard_name, "longitude"
2400+
)
2401+
23852402
def test_Field_laplacian_xy(self):
23862403
f = cf.example_field(0)
23872404

@@ -2439,6 +2456,18 @@ def test_Field_laplacian_xy(self):
24392456
del lp0.long_name
24402457
self.assertTrue(lp.equals(lp0, rtol=1e-10))
24412458

2459+
# Test case when spherical dimension coordinates have units
2460+
# but no standard names
2461+
f = cf.example_field(0)
2462+
del f.dimension_coordinate("X").standard_name
2463+
del f.dimension_coordinate("Y").standard_name
2464+
g = f.laplacian_xy(radius="earth")
2465+
self.assertEqual(g.shape, f.shape)
2466+
self.assertEqual(g.dimension_coordinate("Y").standard_name, "latitude")
2467+
self.assertEqual(
2468+
g.dimension_coordinate("X").standard_name, "longitude"
2469+
)
2470+
24422471
def test_Field_to_dask_array(self):
24432472
f = self.f0.copy()
24442473
self.assertIs(f.to_dask_array(), f.data.to_dask_array())

‎cf/test/test_Maths.py‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -86,10 +86,10 @@ def test_curl_xy(self):
8686
# but no standard names
8787
f = cf.example_field(0)
8888
x, y = f.grad_xy(radius="earth")
89-
del x.dim("X").standard_name
90-
del x.dim("Y").standard_name
91-
del y.dim("X").standard_name
92-
del y.dim("Y").standard_name
89+
del x.dimension_coordinate("X").standard_name
90+
del x.dimension_coordinate("Y").standard_name
91+
del y.dimension_coordinate("X").standard_name
92+
del y.dimension_coordinate("Y").standard_name
9393
c = cf.curl_xy(x, y, radius="earth")
9494
self.assertEqual(c.shape, f.shape)
9595
self.assertEqual(c.dimension_coordinate("Y").standard_name, "latitude")
@@ -174,10 +174,10 @@ def test_div_xy(self):
174174
# but no standard names
175175
f = cf.example_field(0)
176176
x, y = f.grad_xy(radius="earth")
177-
del x.dim("X").standard_name
178-
del x.dim("Y").standard_name
179-
del y.dim("X").standard_name
180-
del y.dim("Y").standard_name
177+
del x.dimension_coordinate("X").standard_name
178+
del x.dimension_coordinate("Y").standard_name
179+
del y.dimension_coordinate("X").standard_name
180+
del y.dimension_coordinate("Y").standard_name
181181
d = cf.div_xy(x, y, radius="earth")
182182
self.assertEqual(d.shape, f.shape)
183183
self.assertEqual(d.dimension_coordinate("Y").standard_name, "latitude")

0 commit comments

Comments
 (0)