Merge pull request #3696 from sgdecker/front_invest

Prevent frontogenesis from returning nans

src/metpy/calc/kinematics.py

@@ -567,9 +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)
-
-    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

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