Merge pull request #3121 from alexander-lakocy/feat/gdi

dopplershift · web-flow · commit 2f8dc029cc60 · 2023-12-27T20:11:12.000-07:00
Add Galvez-Davison Index calculation
diff --git 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
@@ -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
+            <https://www.wpc.ncep.noaa.gov/international/gdi/GDI_Manuscript_V20161021.pdf>`_.
+
 .. [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
@@ -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
@@ -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])
+    <Quantity(-8.78797532, 'dimensionless')>
+    """
+    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
@@ -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])
