diff --git a/docs/_templates/overrides/metpy.calc.rst b/docs/_templates/overrides/metpy.calc.rst index 57a438210f9..e5e7a4f787f 100644 --- a/docs/_templates/overrides/metpy.calc.rst +++ b/docs/_templates/overrides/metpy.calc.rst @@ -79,6 +79,7 @@ Soundings cross_totals downdraft_cape el + galvez_davison_index k_index lcl lfc diff --git a/docs/api/references.rst b/docs/api/references.rst index 787cc6581da..dacd4189546 100644 --- a/docs/api/references.rst +++ b/docs/api/references.rst @@ -86,6 +86,10 @@ References Services and Supporting Research, 2003. `FCM-R19-2003 <../_static/FCM-R19-2003-WindchillReport.pdf>`_, 75 pp. +.. [Galvez2015] Galvez, J. M. and Michel Davison, 2015. “The Gálvez-Davison Index for Tropical + Convection.” `GDI_Manuscript_V20161021 + `_. + .. [Galway1956] Galway, J. G., 1956: The Lifted Index as a Predictor of Latent Instability. *American Meteorology Society*, doi:`10.1175/1520-0477-37.10.528 diff --git a/examples/calculations/Sounding_Calculations.py b/examples/calculations/Sounding_Calculations.py index fad52818007..0791a382ab0 100644 --- a/examples/calculations/Sounding_Calculations.py +++ b/examples/calculations/Sounding_Calculations.py @@ -88,6 +88,11 @@ def effective_layer(p, t, td, h, height_layer=False): sped = df['speed'].values * units.knot height = df['height'].values * units.meter +########################################### +# Compute needed variables from our data file and attach units +relhum = mpcalc.relative_humidity_from_dewpoint(T, Td) +mixrat = mpcalc.mixing_ratio_from_relative_humidity(p, T, relhum) + ########################################### # Compute the wind components u, v = mpcalc.wind_components(sped, wdir) @@ -96,6 +101,7 @@ def effective_layer(p, t, td, h, height_layer=False): # Compute common sounding index parameters ctotals = mpcalc.cross_totals(p, T, Td) kindex = mpcalc.k_index(p, T, Td) +gdi = mpcalc.galvez_davison_index(p, T, mixrat, p[0]) showalter = mpcalc.showalter_index(p, T, Td) total_totals = mpcalc.total_totals_index(p, T, Td) vert_totals = mpcalc.vertical_totals(p, T) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 2a166f874ea..4e7382daabf 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -4497,6 +4497,181 @@ def k_index(pressure, temperature, dewpoint, vertical_dim=0): return ((t850 - t500) + td850 - (t700 - td700)).to(units.degC) +@exporter.export +@add_vertical_dim_from_xarray +@preprocess_and_wrap(broadcast=('pressure', 'temperature', 'mixing_ratio')) +@check_units('[pressure]', '[temperature]', '[dimensionless]', '[pressure]') +def galvez_davison_index(pressure, temperature, mixing_ratio, surface_pressure, + vertical_dim=0): + """ + Calculate GDI from the pressure, temperature, mixing ratio, and surface pressure. + + Calculation of the GDI relies on temperatures and mixing ratios at 950, + 850, 700, and 500 hPa. These four levels define three layers: A) Boundary, + B) Trade Wind Inversion (TWI), C) Mid-Troposphere. + + GDI formula derived from [Galvez2015]_: + + .. math:: GDI = CBI + MWI + II + TC + + where: + + * :math:`CBI` is the Column Buoyancy Index + * :math:`MWI` is the Mid-tropospheric Warming Index + * :math:`II` is the Inversion Index + * :math:`TC` is the Terrain Correction [optional] + + .. list-table:: GDI Values & Corresponding Convective Regimes + :widths: 15 75 + :header-rows: 1 + + * - GDI Value + - Expected Convective Regime + * - >=45 + - Scattered to widespread thunderstorms likely. + * - 35 to 45 + - Scattered thunderstorms and/or scattered to widespread rain showers. + * - 25 to 35 + - Isolated to scattered thunderstorms and/or scattered showers. + * - 15 to 25 + - Isolated thunderstorms and/or isolated to scattered showers. + * - 5 to 10 + - Isolated to scattered showers. + * - <5 + - Strong TWI likely, light rain possible. + + Parameters + ---------- + pressure : `pint.Quantity` + Pressure level(s) + + temperature : `pint.Quantity` + Temperature corresponding to pressure + + mixing_ratio : `pint.Quantity` + Mixing ratio values corresponding to pressure + + surface_pressure : `pint.Quantity` + Pressure of the surface. + + vertical_dim : int, optional + The axis corresponding to vertical, defaults to 0. Automatically determined from + xarray DataArray arguments. + + Returns + ------- + `pint.Quantity` + GDI Index + + Examples + -------- + >>> from metpy.calc import mixing_ratio_from_relative_humidity + >>> from metpy.units import units + >>> # pressure + >>> p = [1008., 1000., 950., 900., 850., 800., 750., 700., 650., 600., + ... 550., 500., 450., 400., 350., 300., 250., 200., + ... 175., 150., 125., 100., 80., 70., 60., 50., + ... 40., 30., 25., 20.] * units.hPa + >>> # temperature + >>> T = [29.3, 28.1, 23.5, 20.9, 18.4, 15.9, 13.1, 10.1, 6.7, 3.1, + ... -0.5, -4.5, -9.0, -14.8, -21.5, -29.7, -40.0, -52.4, + ... -59.2, -66.5, -74.1, -78.5, -76.0, -71.6, -66.7, -61.3, + ... -56.3, -51.7, -50.7, -47.5] * units.degC + >>> # relative humidity + >>> rh = [.85, .65, .36, .39, .82, .72, .75, .86, .65, .22, .52, + ... .66, .64, .20, .05, .75, .76, .45, .25, .48, .76, .88, + ... .56, .88, .39, .67, .15, .04, .94, .35] * units.dimensionless + >>> # calculate mixing ratio + >>> mixrat = mixing_ratio_from_relative_humidity(p, T, rh) + >>> galvez_davison_index(p, T, mixrat, p[0]) + + """ + if np.any(np.max(pressure, axis=vertical_dim) < 950 * units.hectopascal): + indices_without_950 = np.where( + np.max(pressure, axis=vertical_dim) < 950 * units.hectopascal + ) + raise ValueError( + f'Data not provided for 950hPa or higher pressure. ' + f'GDI requires 950hPa temperature and dewpoint data, ' + f'see referenced paper section 3.d. in docstring for discussion of' + f' extrapolating sounding data below terrain surface in high-' + f'elevation regions.\nIndices without a 950hPa or higher datapoint' + f':\n{indices_without_950}' + f'\nMax provided pressures:' + f'\n{np.max(pressure, axis=0)[indices_without_950]}' + ) + + potential_temp = potential_temperature(pressure, temperature) + + # Interpolate to appropriate level with appropriate units + ( + (t950, t850, t700, t500), + (r950, r850, r700, r500), + (th950, th850, th700, th500) + ) = interpolate_1d( + units.Quantity([950, 850, 700, 500], 'hPa'), + pressure, temperature.to('K'), mixing_ratio, potential_temp.to('K'), + axis=vertical_dim, + ) + + # L_v definition preserved from referenced paper + # Value differs heavily from metpy.constants.Lv in tropical context + # and using MetPy value affects resulting GDI + l_0 = units.Quantity(2.69e6, 'J/kg') + + # Calculate adjusted equivalent potential temperatures + alpha = units.Quantity(-10, 'K') + eptp_a = th950 * np.exp(l_0 * r950 / (mpconsts.Cp_d * t850)) + eptp_b = ((th850 + th700) / 2 + * np.exp(l_0 * (r850 + r700) / 2 / (mpconsts.Cp_d * t850)) + alpha) + eptp_c = th500 * np.exp(l_0 * r500 / (mpconsts.Cp_d * t850)) + alpha + + # Calculate Column Buoyanci Index (CBI) + # Apply threshold to low and mid levels + beta = units.Quantity(303, 'K') + l_e = eptp_a - beta + m_e = eptp_c - beta + + # Gamma unit - likely a typo from the paper, should be units of K^(-2) to + # result in dimensionless CBI + gamma = units.Quantity(6.5e-2, '1/K^2') + + column_buoyancy_index = np.atleast_1d(gamma * l_e * m_e) + column_buoyancy_index[l_e <= 0] = 0 + + # Calculate Mid-tropospheric Warming Index (MWI) + # Apply threshold to 500-hPa temperature + tau = units.Quantity(263.15, 'K') + t_diff = t500 - tau + + mu = units.Quantity(-7, '1/K') # Empirical adjustment + mid_tropospheric_warming_index = np.atleast_1d(mu * t_diff) + mid_tropospheric_warming_index[t_diff <= 0] = 0 + + # Calculate Inversion Index (II) + s = t950 - t700 + d = eptp_b - eptp_a + inv_sum = s + d + + sigma = units.Quantity(1.5, '1/K') # Empirical scaling constant + inversion_index = np.atleast_1d(sigma * inv_sum) + inversion_index[inv_sum >= 0] = 0 + + # Calculate Terrain Correction + terrain_correction = 18 - 9000 / (surface_pressure.m_as('hPa') - 500) + + # Calculate G.D.I. + gdi = (column_buoyancy_index + + mid_tropospheric_warming_index + + inversion_index + + terrain_correction) + + if gdi.size == 1: + return gdi[0] + else: + return gdi + + @exporter.export @add_vertical_dim_from_xarray @preprocess_and_wrap( diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index ea97646d49d..4cd9c60749b 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -37,7 +37,7 @@ vertical_velocity_pressure, virtual_potential_temperature, virtual_temperature, virtual_temperature_from_dewpoint, wet_bulb_potential_temperature, wet_bulb_temperature) -from metpy.calc.thermo import _find_append_zero_crossings +from metpy.calc.thermo import _find_append_zero_crossings, galvez_davison_index from metpy.testing import (assert_almost_equal, assert_array_almost_equal, assert_nan, version_check) from metpy.units import is_quantity, masked_array, units @@ -2309,6 +2309,50 @@ def index_xarray_data(): coords={'isobaric': pressure, 'time': ['2020-01-01T00:00Z']}) +@pytest.fixture() +def index_xarray_data_expanded(): + """Create expanded data for testing that index calculations work with xarray data. + + Specifically for Galvez Davison Index calculation, which requires 950hPa pressure + """ + pressure = xr.DataArray( + [950., 850., 700., 500.], dims=('isobaric',), attrs={'units': 'hPa'} + ) + temp = xr.DataArray([[[[306., 305., 304.], [303., 302., 301.]], + [[296., 295., 294.], [293., 292., 291.]], + [[286., 285., 284.], [283., 282., 281.]], + [[276., 275., 274.], [273., 272., 271.]]]] * units.K, + dims=('time', 'isobaric', 'y', 'x')) + + profile = xr.DataArray([[[[299., 298., 297.], [296., 295., 294.]], + [[289., 288., 287.], [286., 285., 284.]], + [[279., 278., 277.], [276., 275., 274.]], + [[269., 268., 267.], [266., 265., 264.]]]] * units.K, + dims=('time', 'isobaric', 'y', 'x')) + + dewp = xr.DataArray([[[[304., 303., 302.], [301., 300., 299.]], + [[294., 293., 292.], [291., 290., 289.]], + [[284., 283., 282.], [281., 280., 279.]], + [[274., 273., 272.], [271., 270., 269.]]]] * units.K, + dims=('time', 'isobaric', 'y', 'x')) + + dirw = xr.DataArray([[[[135., 135., 135.], [135., 135., 135.]], + [[180., 180., 180.], [180., 180., 180.]], + [[225., 225., 225.], [225., 225., 225.]], + [[270., 270., 270.], [270., 270., 270.]]]] * units.degree, + dims=('time', 'isobaric', 'y', 'x')) + + speed = xr.DataArray([[[[15., 15., 15.], [15., 15., 15.]], + [[20., 20., 20.], [20., 20., 20.]], + [[25., 25., 25.], [25., 25., 25.]], + [[50., 50., 50.], [50., 50., 50.]]]] * units.knots, + dims=('time', 'isobaric', 'y', 'x')) + + return xr.Dataset({'temperature': temp, 'profile': profile, 'dewpoint': dewp, + 'wind_direction': dirw, 'wind_speed': speed}, + coords={'isobaric': pressure, 'time': ['2023-01-01T00:00Z']}) + + def test_lifted_index(): """Test the Lifted Index calculation.""" pressure = np.array([1014., 1000., 997., 981.2, 947.4, 925., 914.9, 911., @@ -2398,6 +2442,109 @@ def test_k_index_xarray(index_xarray_data): np.array([[[312., 311., 310.], [309., 308., 307.]]]) * units.K) +def test_gdi(): + """Test the Galvez Davison Index calculation.""" + pressure = np.array([1014., 1000., 997., 981.2, 947.4, 925., 914.9, 911., + 902., 883., 850., 822.3, 816., 807., 793.2, 770., + 765.1, 753., 737.5, 737., 713., 700., 688., 685., + 680., 666., 659.8, 653., 643., 634., 615., 611.8, + 566.2, 516., 500., 487., 484.2, 481., 475., 460., + 400.]) * units.hPa + temperature = np.array([24.2, 24.2, 24., 23.1, 21., 19.6, 18.7, 18.4, + 19.2, 19.4, 17.2, 15.3, 14.8, 14.4, 13.4, 11.6, + 11.1, 10., 8.8, 8.8, 8.2, 7., 5.6, 5.6, + 5.6, 4.4, 3.8, 3.2, 3., 3.2, 1.8, 1.5, + -3.4, -9.3, -11.3, -13.1, -13.1, -13.1, -13.7, -15.1, + -23.5]) * units.degC + dewpoint = np.array([23.2, 23.1, 22.8, 22., 20.2, 19., 17.6, 17., + 16.8, 15.5, 14., 11.7, 11.2, 8.4, 7., 4.6, + 5., 6., 4.2, 4.1, -1.8, -2., -1.4, -0.4, + -3.4, -5.6, -4.3, -2.8, -7., -25.8, -31.2, -31.4, + -34.1, -37.3, -32.3, -34.1, -37.3, -41.1, -37.7, -58.1, + -57.5]) * units.degC + + relative_humidity = relative_humidity_from_dewpoint(temperature, dewpoint) + mixrat = mixing_ratio_from_relative_humidity(pressure, temperature, relative_humidity) + gdi = galvez_davison_index(pressure, temperature, mixrat, pressure[0]) + + # Compare with value from hand calculation + assert_almost_equal(gdi, 6.635, decimal=1) + + +def test_gdi_xarray(index_xarray_data_expanded): + """Test the GDI calculation with a grid of xarray data.""" + pressure = index_xarray_data_expanded.isobaric + temperature = index_xarray_data_expanded.temperature + dewpoint = index_xarray_data_expanded.dewpoint + mixing_ratio = mixing_ratio_from_relative_humidity( + pressure, temperature, relative_humidity_from_dewpoint(temperature, dewpoint)) + + result = galvez_davison_index( + pressure, + temperature, + mixing_ratio, + pressure[0] + ) + + assert_array_almost_equal( + result, + np.array([[[189.5890429, 157.4307982, 129.9739099], + [106.6763526, 87.0637477, 70.7202505]]]) + ) + + +def test_gdi_arrays(index_xarray_data_expanded): + """Test GDI on 3-D Quantity arrays with an array of surface pressure.""" + ds = index_xarray_data_expanded.isel(time=0).squeeze() + pressure = ds.isobaric.metpy.unit_array[:, None, None] + temperature = ds.temperature.metpy.unit_array + dewpoint = ds.dewpoint.metpy.unit_array + mixing_ratio = mixing_ratio_from_relative_humidity( + pressure, temperature, relative_humidity_from_dewpoint(temperature, dewpoint)) + surface_pressure = units.Quantity( + np.broadcast_to(pressure.m, temperature.shape), pressure.units)[0] + + result = galvez_davison_index(pressure, temperature, mixing_ratio, surface_pressure) + + assert_array_almost_equal( + result, + np.array([[189.5890429, 157.4307982, 129.9739099], + [106.6763526, 87.0637477, 70.7202505]]) + ) + + +def test_gdi_profile(index_xarray_data_expanded): + """Test GDI calculation on an individual profile.""" + ds = index_xarray_data_expanded.isel(time=0, y=0, x=0) + pressure = ds.isobaric.metpy.unit_array + temperature = ds.temperature.metpy.unit_array + dewpoint = ds.dewpoint.metpy.unit_array + mixing_ratio = mixing_ratio_from_relative_humidity( + pressure, temperature, relative_humidity_from_dewpoint(temperature, dewpoint)) + + assert_almost_equal(galvez_davison_index(pressure, temperature, mixing_ratio, pressure[0]), + 189.5890429, 4) + + +def test_gdi_no_950_raises_valueerror(index_xarray_data): + """GDI requires a 950hPa or higher measurement. + + Ensure error is raised if this data is not provided. + """ + with pytest.raises(ValueError): + pressure = index_xarray_data.isobaric + temperature = index_xarray_data.temperature + dewpoint = index_xarray_data.dewpoint + relative_humidity = relative_humidity_from_dewpoint(temperature, dewpoint) + mixrat = mixing_ratio_from_relative_humidity(pressure, temperature, relative_humidity) + galvez_davison_index( + pressure, + temperature, + mixrat, + pressure[0] + ) + + def test_gradient_richardson_number(): """Test gradient Richardson number calculation.""" theta = units('K') * np.asarray([254.5, 258.3, 262.2])