Allow PW calculation when NaNs present

src/metpy/calc/indices.py

@@ -1,11 +1,11 @@
-# Copyright (c) 2017 MetPy Developers.
+# Copyright (c) 2017,2019 MetPy Developers.
 # Distributed under the terms of the BSD 3-Clause License.
 # SPDX-License-Identifier: BSD-3-Clause
 """Contains calculation of various derived indices."""
 import numpy as np
 
 from .thermo import mixing_ratio, saturation_vapor_pressure
-from .tools import get_layer
+from .tools import _remove_nans, get_layer
 from .. import constants as mpconsts
 from ..package_tools import Exporter
 from ..units import atleast_1d, check_units, concatenate, units
@@ -48,6 +48,8 @@ def precipitable_water(dewpt, pressure, bottom=None, top=None):
     pressure = pressure[sort_inds]
     dewpt = dewpt[sort_inds]
 
+    pressure, dewpt = _remove_nans(pressure, dewpt)
+
     if top is None:
         top = np.nanmin(pressure.magnitude) * pressure.units
 

src/metpy/calc/thermo.py

@@ -11,7 +11,7 @@
 import scipy.integrate as si
 import scipy.optimize as so
 
-from .tools import (_greater_or_close, _less_or_close, find_bounding_indices,
+from .tools import (_greater_or_close, _less_or_close, _remove_nans, find_bounding_indices,
                     find_intersections, first_derivative, get_layer)
 from .. import constants as mpconsts
 from ..cbook import broadcast_indices
@@ -407,6 +407,7 @@ def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None, dewpt_sta
     parcel_profile
 
     """
+    pressure, temperature, dewpt = _remove_nans(pressure, temperature, dewpt)
     # Default to surface parcel if no profile or starting pressure level is given
     if parcel_temperature_profile is None:
         new_stuff = parcel_profile_with_lcl(pressure, temperature, dewpt)
@@ -514,6 +515,7 @@ def el(pressure, temperature, dewpt, parcel_temperature_profile=None, which='top
     parcel_profile
 
     """
+    pressure, temperature, dewpt = _remove_nans(pressure, temperature, dewpt)
     # Default to surface parcel if no profile or starting pressure level is given
     if parcel_temperature_profile is None:
         new_stuff = parcel_profile_with_lcl(pressure, temperature, dewpt)
@@ -1430,6 +1432,8 @@ def cape_cin(pressure, temperature, dewpt, parcel_profile):
     lfc, el
 
     """
+    pressure, temperature, dewpt, parcel_profile = _remove_nans(pressure, temperature, dewpt,
+                                                                parcel_profile)
     # Calculate LFC limit of integration
     lfc_pressure, _ = lfc(pressure, temperature, dewpt,
                           parcel_temperature_profile=parcel_profile)
@@ -1757,6 +1761,7 @@ def surface_based_cape_cin(pressure, temperature, dewpoint):
     cape_cin, parcel_profile
 
     """
+    pressure, temperature, dewpoint = _remove_nans(pressure, temperature, dewpoint)
     p, t, td, profile = parcel_profile_with_lcl(pressure, temperature, dewpoint)
     return cape_cin(p, t, td, profile)
 
@@ -1794,6 +1799,7 @@ def most_unstable_cape_cin(pressure, temperature, dewpoint, **kwargs):
     cape_cin, most_unstable_parcel, parcel_profile
 
     """
+    pressure, temperature, dewpoint = _remove_nans(pressure, temperature, dewpoint)
     _, _, _, parcel_idx = most_unstable_parcel(pressure, temperature, dewpoint, **kwargs)
     p, t, td, mu_profile = parcel_profile_with_lcl(pressure[parcel_idx:],
                                                    temperature[parcel_idx:],

src/metpy/calc/tools.py

@@ -1480,3 +1480,23 @@ def _unabbrieviate_direction(abb_dir_str):
             .replace('W', 'West ')
             .replace(' ,', ',')
             ).strip()
+
+
+def _remove_nans(*variables):
+    """Remove NaNs from arrays that cause issues with calculations.
+
+    Takes a variable number of arguments
+    Returns masked arrays in the same order as provided
+    """
+    mask = None
+    for v in variables:
+        if mask is None:
+            mask = np.isnan(v)
+        else:
+            mask |= np.isnan(v)
+
+    # Mask everyone with that joint mask
+    ret = []
+    for v in variables:
+        ret.append(v[~mask])
+    return ret

tests/calc/test_calc_tools.py

@@ -1,4 +1,4 @@
-# Copyright (c) 2016,2017,2018 MetPy Developers.
+# Copyright (c) 2016,2017,2018,2019 MetPy Developers.
 # Distributed under the terms of the BSD 3-Clause License.
 # SPDX-License-Identifier: BSD-3-Clause
 """Test the `tools` module."""
@@ -19,7 +19,7 @@
                         reduce_point_density, resample_nn_1d, second_derivative)
 from metpy.calc.tools import (_delete_masked_points, _get_bound_pressure_height,
                               _greater_or_close, _less_or_close, _next_non_masked_element,
-                              BASE_DEGREE_MULTIPLIER, DIR_STRS, UND)
+                              _remove_nans, BASE_DEGREE_MULTIPLIER, DIR_STRS, UND)
 from metpy.deprecation import MetpyDeprecationWarning
 from metpy.testing import (assert_almost_equal, assert_array_almost_equal, assert_array_equal,
                            check_and_silence_deprecation)
@@ -1190,3 +1190,14 @@ def test_gradient_xarray_pint_conversion(test_da_xy):
 
     assert_array_almost_equal(deriv_x, truth_x, 12)
     assert_array_almost_equal(deriv_y, truth_y, 12)
+
+
+def test_remove_nans():
+    """Test removal of NaNs."""
+    x = np.array([3, 2, np.nan, 5, 6, np.nan])
+    y = np.arange(0, len(x))
+    y_test, x_test = _remove_nans(y, x)
+    x_expected = np.array([3, 2, 5, 6])
+    y_expected = np.array([0, 1, 3, 4])
+    assert_array_almost_equal(x_expected, x_test, 0)
+    assert_almost_equal(y_expected, y_test, 0)

tests/calc/test_indices.py

@@ -1,4 +1,4 @@
-# Copyright (c) 2017 MetPy Developers.
+# Copyright (c) 2017,2019 MetPy Developers.
 # Distributed under the terms of the BSD 3-Clause License.
 # SPDX-License-Identifier: BSD-3-Clause
 """Test the `indices` module."""
@@ -45,6 +45,22 @@ def test_precipitable_water_bound_error():
     assert_almost_equal(pw, truth, 8)
 
 
+def test_precipitable_water_nans():
+    """Test that PW returns appropriate number if NaNs are present."""
+    pressure = np.array([1001, 1000, 997, 977.9, 977, 957, 937.8, 925, 906, 899.3, 887, 862.5,
+                         854, 850, 800, 793.9, 785, 777, 771, 762, 731.8, 726, 703, 700, 655,
+                         630, 621.2, 602, 570.7, 548, 546.8, 539, 513, 511, 485, 481, 468,
+                         448, 439, 424, 420, 412]) * units.hPa
+    dewpoint = np.array([-25.1, -26.1, -26.8, np.nan, -27.3, -28.2, np.nan, -27.2, -26.6,
+                        np.nan, -27.4, np.nan, -23.5, -23.5, -25.1, np.nan, -22.9, -17.8,
+                        -16.6, np.nan, np.nan, -16.4, np.nan, -18.5, -21., -23.7, np.nan,
+                        -28.3, np.nan, -32.6, np.nan, -33.8, -35., -35.1, -38.1, -40.,
+                        -43.3, -44.6, -46.4, -47., -49.2, -50.7]) * units.degC
+    pw = precipitable_water(dewpoint, pressure)
+    truth = 4.003709214463873 * units.mm
+    assert_almost_equal(pw, truth, 8)
+
+
 def test_mean_pressure_weighted():
     """Test pressure-weighted mean wind function with vertical interpolation."""
     data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC')

tests/calc/test_thermo.py

@@ -1032,6 +1032,37 @@ def test_surface_based_cape_cin():
     assert_almost_equal(cin, -136.607809 * units('joule / kilogram'), 2)
 
 
+def test_profile_with_nans():
+    """Test a profile with nans to make sure it calculates functions appropriately (#1187)."""
+    pressure = np.array([1001, 1000, 997, 977.9, 977, 957, 937.8, 925, 906, 899.3, 887, 862.5,
+                         854, 850, 800, 793.9, 785, 777, 771, 762, 731.8, 726, 703, 700, 655,
+                         630, 621.2, 602, 570.7, 548, 546.8, 539, 513, 511, 485, 481, 468,
+                         448, 439, 424, 420, 412]) * units.hPa
+    temperature = np.array([-22.5, -22.7, -23.1, np.nan, -24.5, -25.1, np.nan, -24.5, -23.9,
+                            np.nan, -24.7, np.nan, -21.3, -21.3, -22.7, np.nan, -20.7, -16.3,
+                            -15.5, np.nan, np.nan, -15.3, np.nan, -17.3, -20.9, -22.5,
+                            np.nan, -25.5, np.nan, -31.5, np.nan, -31.5, -34.1, -34.3,
+                            -37.3, -37.7, -39.5, -42.1, -43.1, -45.1, -45.7, -46.7]
+                           ) * units.degC
+    dewpoint = np.array([-25.1, -26.1, -26.8, np.nan, -27.3, -28.2, np.nan, -27.2, -26.6,
+                         np.nan, -27.4, np.nan, -23.5, -23.5, -25.1, np.nan, -22.9, -17.8,
+                         -16.6, np.nan, np.nan, -16.4, np.nan, -18.5, -21, -23.7, np.nan,
+                         -28.3, np.nan, -32.6, np.nan, -33.8, -35, -35.1, -38.1, -40,
+                         -43.3, -44.6, -46.4, -47, -49.2, -50.7]) * units.degC
+    lfc_p, _ = lfc(pressure, temperature, dewpoint)
+    profile = parcel_profile(pressure, temperature[0], dewpoint[0])
+    cape, cin = cape_cin(pressure, temperature, dewpoint, profile)
+    sbcape, sbcin = surface_based_cape_cin(pressure, temperature, dewpoint)
+    mucape, mucin = most_unstable_cape_cin(pressure, temperature, dewpoint)
+    assert_nan(lfc_p, units.hPa)
+    assert_almost_equal(cape, 0 * units('J/kg'), 0)
+    assert_almost_equal(cin, 0 * units('J/kg'), 0)
+    assert_almost_equal(sbcape, 0 * units('J/kg'), 0)
+    assert_almost_equal(sbcin, 0 * units('J/kg'), 0)
+    assert_almost_equal(mucape, 0 * units('J/kg'), 0)
+    assert_almost_equal(mucin, 0 * units('J/kg'), 0)
+
+
 def test_most_unstable_cape_cin_surface():
     """Test the most unstable CAPE/CIN calculation when surface is most unstable."""
     pressure = np.array([959., 779.2, 751.3, 724.3, 700., 269.]) * units.mbar