diff --git a/src/metpy/calc/indices.py b/src/metpy/calc/indices.py index f0af5e13836..60bdfc47f6d 100644 --- a/src/metpy/calc/indices.py +++ b/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 @@ -45,11 +45,10 @@ def precipitable_water(dewpt, pressure, bottom=None, top=None): """ # Sort pressure and dewpoint to be in decreasing pressure order (increasing height) sort_inds = np.argsort(pressure)[::-1] - pressure_inter = pressure[sort_inds] - dewpt_inter = dewpt[sort_inds] + pressure = pressure[sort_inds] + dewpt = dewpt[sort_inds] - dewpt = dewpt_inter[~np.isnan(dewpt_inter) & ~np.isnan(pressure_inter)] - pressure = pressure[~np.isnan(dewpt_inter) & ~np.isnan(pressure_inter)] + pressure, dewpt = _remove_nans(pressure, dewpt) if top is None: top = np.nanmin(pressure.magnitude) * pressure.units diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index b14de43002d..bf4e1fdc78e 100644 --- a/src/metpy/calc/thermo.py +++ b/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 @@ -408,6 +408,8 @@ def lfc(pressure, temperature, dewpt, parcel_temperature_profile=None, dewpt_sta parcel_profile """ + _, temperature = _remove_nans(pressure, temperature) + pressure, dewpt = _remove_nans(pressure, 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) @@ -515,6 +517,8 @@ def el(pressure, temperature, dewpt, parcel_temperature_profile=None, which='top parcel_profile """ + _, temperature = _remove_nans(pressure, temperature) + pressure, dewpt = _remove_nans(pressure, 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) @@ -1431,6 +1435,10 @@ def cape_cin(pressure, temperature, dewpt, parcel_profile): lfc, el """ + _, parcel_profile = _remove_nans(temperature, parcel_profile) + _, temperature = _remove_nans(pressure, temperature) + pressure, dewpt = _remove_nans(pressure, dewpt) + # Calculate LFC limit of integration lfc_pressure, _ = lfc(pressure, temperature, dewpt, parcel_temperature_profile=parcel_profile) @@ -1758,6 +1766,8 @@ def surface_based_cape_cin(pressure, temperature, dewpoint): cape_cin, parcel_profile """ + _, temperature = _remove_nans(pressure, temperature) + pressure, dewpoint = _remove_nans(pressure, dewpoint) p, t, td, profile = parcel_profile_with_lcl(pressure, temperature, dewpoint) return cape_cin(p, t, td, profile) @@ -1795,6 +1805,8 @@ def most_unstable_cape_cin(pressure, temperature, dewpoint, **kwargs): cape_cin, most_unstable_parcel, parcel_profile """ + _, temperature = _remove_nans(pressure, temperature) + pressure, dewpoint = _remove_nans(pressure, 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:], diff --git a/src/metpy/calc/tools.py b/src/metpy/calc/tools.py index 00fd0502869..999cfb6d016 100644 --- a/src/metpy/calc/tools.py +++ b/src/metpy/calc/tools.py @@ -1480,3 +1480,14 @@ def _unabbrieviate_direction(abb_dir_str): .replace('W', 'West ') .replace(' ,', ',') ).strip() + + +def _remove_nans(pressure, nan_variable): + """Remove NaNs from arrays that cause issues with calculations. + + Takes a variable number of arguments + Returns a dictionary with arrays in the same order as provided + """ + no_nan_variable = nan_variable[~np.isnan(nan_variable) & ~np.isnan(pressure)] + pressure = pressure[~np.isnan(nan_variable) & ~np.isnan(pressure)] + return pressure, no_nan_variable diff --git a/tests/calc/test_calc_tools.py b/tests/calc/test_calc_tools.py index 80a62844f2c..82580cf128d 100644 --- a/tests/calc/test_calc_tools.py +++ b/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) diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index 6383efa6413..1a548a682d8 100644 --- a/tests/calc/test_thermo.py +++ b/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