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

@@ -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)
 

tests/calc/test_kinematics.py

@@ -307,6 +307,53 @@ 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')