From 12f42851827cc7a9e7ae49789bcb4bcf03151c5d Mon Sep 17 00:00:00 2001 From: Steve Decker Date: Tue, 19 Nov 2024 10:42:01 -0500 Subject: [PATCH 1/2] Prevent frontogenesis from returning nans Fixes #3768 by making sure the argument to the arcsin function is valid. Previously, frontogenesis could return nans when there was a constant theta field (division by zero would occur) or if round-off error resulted in the argument to arcsin being slightly outside the valid domain of the function (-1 to 1). In this commit, edits are made to set points to zero where nans occur due to division by zero (the frontogenesis is zero when the magnitude of the theta gradient is zero anyway) and to use np.clip to ensure the argument to arcsin is valid. I could not come up with a simplified test case that triggers the round-off error issue with arcsin, but I do include a test case for the constant theta situation. Because the test case results in a division by zero by design, it is currently failing since that triggers a RuntimeWarning. --- src/metpy/calc/kinematics.py | 13 ++++++++- tests/calc/test_kinematics.py | 53 +++++++++++++++++++++++++++++++++++ 2 files changed, 65 insertions(+), 1 deletion(-) diff --git a/src/metpy/calc/kinematics.py b/src/metpy/calc/kinematics.py index 73f8d7a64c4..218e4ba9538 100644 --- a/src/metpy/calc/kinematics.py +++ b/src/metpy/calc/kinematics.py @@ -567,7 +567,18 @@ def frontogenesis(potential_temperature, u, v, dx=None, dy=None, x_dim=-1, y_dim # Compute the angle (beta) between the wind field and the gradient of potential temperature psi = 0.5 * np.arctan2(shrd, strd) - beta = np.arcsin((-ddx_theta * np.cos(psi) - ddy_theta * np.sin(psi)) / mag_theta) + arg = (-ddx_theta * np.cos(psi) - ddy_theta * np.sin(psi)) / mag_theta + + # A few problems may occur when calculating the argument to the arcsin function. + # First, we may have divided by zero, since a constant theta field would mean + # mag_theta is zero. To counter this, we set the argument to zero in this case. + # Second, due to round-off error, the argument may be slightly outside the domain + # of arcsin. To counter this, we use np.clip to force the argument to be an + # acceptable value. With these adjustments, we can make sure beta doesn't end up + # with nans somewhere. + arg[mag_theta == 0] = 0 + arg = np.clip(arg, -1, 1) + beta = np.arcsin(arg) return 0.5 * mag_theta * (tdef * np.cos(2 * beta) - div) diff --git a/tests/calc/test_kinematics.py b/tests/calc/test_kinematics.py index f71a75bc530..71a62069c2a 100644 --- a/tests/calc/test_kinematics.py +++ b/tests/calc/test_kinematics.py @@ -307,6 +307,59 @@ def test_frontogenesis_asym(): assert_almost_equal(fronto, true_fronto) +def test_frontogenesis_nan(): + """Test that frontogenesis calculation does not result in nan.""" + x = np.array([-4142340.8, -4061069.8, -3979798.8, -3898527.8, -3817256.8], + dtype=np.float32) + y = np.array([-832207.56, -750936.56, -669665.56, -588394.5, -507123.56], + dtype=np.float32) + lat = np.array([[12.38805122, 12.58268367, 12.77387305, 12.96159447, 13.14582354], + [13.07545967, 13.27159197, 13.46425116, 13.65341235, 13.83905111], + [13.76423766, 13.96186003, 14.15597929, 14.34657051, 14.53360928], + [14.45429168, 14.65339375, 14.84896275, 15.04097373, 15.22940228], + [15.14552468, 15.3460955, 15.54310332, 15.73652321, 15.92633074]]) + lon = np.array([[-132.75696788, -132.05286812, -131.34671228, -130.63852983, + -129.92835084], + [-132.9590417, -132.25156385, -131.54199837, -130.83037505, + -130.11672431], + [-133.16323241, -132.45234731, -131.73934239, -131.02424779, + -130.30709426], + [-133.36957268, -132.65525085, -131.93877637, -131.22017972, + -130.49949199], + [-133.57809517, -132.86030681, -132.14033233, -131.41820252, + -130.69394884]]) + uvals = np.array([[0.165024, 0.055023, -0.454977, -1.534977, -2.744976], + [0.155024, -0.434977, -2.114977, -3.474977, -4.034977], + [-1.554977, -2.714977, -2.084976, -5.274977, -3.334976], + [-3.424976, -7.644977, -7.654977, -5.384976, -3.224977], + [-9.564977, -9.934977, -7.454977, -6.004977, -4.144977]]) * units('m/s') + vvals = (np.array([[2.6594982, 1.9194984, 2.979498, 2.149498, 2.6394978], + [3.4994984, 4.0794983, 4.8494987, 5.2594986, 3.1694984], + [5.159498, 6.4594975, 6.559498, 5.9694977, 3.189499], + [6.5994987, 9.799498, 7.4594975, 4.2894993, 3.729498], + [11.3394985, 6.779499, 4.0994987, 4.819498, 4.9994984]]) + * units('m/s')) + tvals = np.array([[290.2, 290.1, 290.2, 290.30002, 290.30002], + [290.5, 290.30002, 290.30002, 290.30002, 290.2], + [290.80002, 290.40002, 290.2, 290.40002, 289.90002], + [290.7, 290.90002, 290.7, 290.1, 289.7], + [290.90002, 290.40002, 289.7, 289.7, 289.30002]]) * units('degK') + + x = xr.DataArray(data=x, coords={'x': x}, dims='x', attrs={'units': 'meters'}) + y = xr.DataArray(data=y, coords={'y': y}, dims='y', attrs={'units': 'meters'}) + + dims = ['y', 'x'] + coords = {'x': x, 'y': y, 'latitude': (dims, lat), 'longitude': (dims, lon)} + + u = xr.DataArray(data=uvals, coords=coords, dims=dims) + v = xr.DataArray(data=vvals, coords=coords, dims=dims) + t = xr.DataArray(data=tvals, coords=coords, dims=dims) + + th = potential_temperature(850 * units('hPa'), t) + f = frontogenesis(th, u, v) + assert not np.any(np.isnan(f)) + + def test_advection_uniform(): """Test advection calculation for a uniform 1D field.""" u = np.ones((3,)) * units('m/s') From 7a820ef4dfc48ac4c55621dfcb1bbbcd08ec7aa8 Mon Sep 17 00:00:00 2001 From: Steve Decker Date: Thu, 21 Nov 2024 10:08:23 -0500 Subject: [PATCH 2/2] Avoid nans more accurately Cleaning up numerical issues like division by zero and round-off error after the fact are kludges, so this commit replaces the previous approach with more careful numerics. --- src/metpy/calc/kinematics.py | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/src/metpy/calc/kinematics.py b/src/metpy/calc/kinematics.py index 218e4ba9538..dc2d8275bb7 100644 --- a/src/metpy/calc/kinematics.py +++ b/src/metpy/calc/kinematics.py @@ -567,20 +567,18 @@ def frontogenesis(potential_temperature, u, v, dx=None, dy=None, x_dim=-1, y_dim # Compute the angle (beta) between the wind field and the gradient of potential temperature psi = 0.5 * np.arctan2(shrd, strd) - arg = (-ddx_theta * np.cos(psi) - ddy_theta * np.sin(psi)) / mag_theta - - # A few problems may occur when calculating the argument to the arcsin function. - # First, we may have divided by zero, since a constant theta field would mean - # mag_theta is zero. To counter this, we set the argument to zero in this case. - # Second, due to round-off error, the argument may be slightly outside the domain - # of arcsin. To counter this, we use np.clip to force the argument to be an - # acceptable value. With these adjustments, we can make sure beta doesn't end up - # with nans somewhere. - arg[mag_theta == 0] = 0 - arg = np.clip(arg, -1, 1) - beta = np.arcsin(arg) - - return 0.5 * mag_theta * (tdef * np.cos(2 * beta) - div) + # We need to be careful to avoid division by zero. When mag_theta + # is zero, the frontogenesis will also be zero. The minus signs + # are omitted from the numerator since this expression is squared + # to compute the frontogenesis. + sin_beta = np.divide(ddx_theta * np.cos(psi) + ddy_theta * np.sin(psi), mag_theta, + out=np.zeros_like(mag_theta), where=mag_theta != 0) + + # The textbook definition of frontogenesis includes the term + # cos(2*beta). However, using trig identities, one can show that + # cos(2*beta) = 1 - 2 * sin(beta)**2, and the expression involving + # sin(beta) is more numerically stable. + return 0.5 * mag_theta * (tdef * (1 - 2 * sin_beta**2) - div) @exporter.export