From a2f713ceeaabf1dafdbb56e68bbf59397d040bbb Mon Sep 17 00:00:00 2001 From: Steve Decker Date: Tue, 19 Nov 2024 10:42:01 -0500 Subject: [PATCH] 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..6b1a2386484 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..e54b3b0a654 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')