make changes to GDI calc

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,7 +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, 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)

src/metpy/calc/thermo.py

@@ -4333,11 +4333,12 @@ def k_index(pressure, temperature, dewpoint, vertical_dim=0):
 
 @exporter.export
 @add_vertical_dim_from_xarray
-@preprocess_and_wrap(broadcast=('pressure', 'temperature', 'dewpoint'))
-@check_units('[pressure]', '[temperature]', '[temperature]')
-def galvez_davison_index(pressure, temperature, dewpoint, vertical_dim=0):
+@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 and dewpoint.
+    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,
@@ -4376,13 +4377,16 @@ def galvez_davison_index(pressure, temperature, dewpoint, vertical_dim=0):
     Parameters
     ----------
     pressure : `pint.Quantity`
-        Pressure level(s), in order from highest to lowest pressure
+        Pressure level(s)
 
     temperature : `pint.Quantity`
         Temperature corresponding to pressure
 
-    dewpoint : `pint.Quantity`
-        Dewpoint 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
@@ -4411,19 +4415,11 @@ def galvez_davison_index(pressure, temperature, dewpoint, vertical_dim=0):
     >>> 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 dewpoint
-    >>> Td = dewpoint_from_relative_humidity(T, rh)
-    >>> galvez_davison_index(p, T, Td)
+    >>> # calculate mixing ratio
+    >>> mixrat = mixing_ratio_from_relative_humidity(p, T, rh)
+    >>> galvez_davison_index(p, T, mixrat, p[0])
     <Quantity(-8.06886508, 'dimensionless')>
     """
-    # Calculate mixing ratio from dewpoint in two steps
-    relative_humidity = relative_humidity_from_dewpoint(
-        temperature, dewpoint
-    )
-    mixing_ratio = mixing_ratio_from_relative_humidity(
-        pressure, temperature, relative_humidity
-    )
-
     # Calculate potential temperature
     potential_temp = potential_temperature(pressure, temperature)
 
@@ -4465,12 +4461,12 @@ def galvez_davison_index(pressure, temperature, dewpoint, vertical_dim=0):
     th_c = th500
     r_c = r500
 
-    alpha = -10 * units.kelvin  # Empirical adjustment
+    alpha = units.Quantity(-10, 'K')  # Empirical adjustment
 
     # Latent heat of vaporization of water - different from MetPy constant
     #  Using different value, from paper, since GDI unlikely to be used to
     #  derive other metrics
-    l_0 = 2.69E6 * (units.joule / units.kilogram)
+    l_0 = units.Quantity(2.69E6, 'J/kg')
 
     # Temperature math from here on requires kelvin units
     eptp_a = th_a * np.exp((l_0 * r_a) / (Cp_d * t850))
@@ -4483,14 +4479,14 @@ def galvez_davison_index(pressure, temperature, dewpoint, vertical_dim=0):
         is_array = True
 
     # Calculate C.B.I.
-    beta = 303 * units.kelvin  # Empirical adjustment
+    beta = units.Quantity(303, 'K')  # Empirical adjustment
     l_e = eptp_a - beta  # Low-troposphere EPT
     m_e = eptp_c - beta  # Mid-troposphere EPT
 
     # Gamma unit - likely a typo from the paper, should be units of K^(-2) to
-    #  result in dimensionless CBI
-    gamma = 6.5e-2 * (1 / units.kelvin)
-    zero_kelvin = 0 * units.kelvin
+    # result in dimensionless CBI
+    gamma = units.Quantity(6.5e-2, '1/K^2')
+    zero_kelvin = units.Quantity(0, 'K')
 
     if is_array:
         # Replace conditional in paper for array compatibility.
@@ -4499,23 +4495,20 @@ def galvez_davison_index(pressure, temperature, dewpoint, vertical_dim=0):
     else:
         l_e = max(l_e, zero_kelvin)
     column_buoyancy_index = gamma * (l_e * m_e)
-    # Convert to magnitude and dimensionless to avoid unit issue from typo
-    column_buoyancy_index = column_buoyancy_index.magnitude
 
     # Calculate Mid-tropospheric Warming Index
-    tau = 263.15 * units.kelvin  # Threshold
-    mu = -7 * (1 / units.kelvin)  # Empirical adjustment
+    tau = units.Quantity(263.15, 'K')  # Threshold
+    mu = units.Quantity(-7, '1/K')  # Empirical adjustment
 
     t_diff = t500 - tau
     if is_array:
         t_diff[t_diff <= zero_kelvin] = zero_kelvin
     else:
         t_diff = max(t_diff, zero_kelvin)
     mid_tropospheric_warming_index = mu * t_diff
-    mid_tropospheric_warming_index = mid_tropospheric_warming_index.magnitude
 
     # Calculate Inversion Index
-    sigma = 1.5 * (1 / units.kelvin)  # Empirical scaling constant
+    sigma = units.Quantity(1.5, '1/K')  # Empirical scaling constant
     s = t950 - t700
     d = eptp_b - eptp_a
 
@@ -4526,21 +4519,14 @@ def galvez_davison_index(pressure, temperature, dewpoint, vertical_dim=0):
         inv_sum = min(inv_sum, zero_kelvin)
 
     inversion_index = sigma * inv_sum
-    inversion_index = inversion_index.magnitude
 
     # Calculate Terrain Correction
     p_3 = 18
     p_2 = 9000 * units.hectopascal
     p_1 = 500 * units.hectopascal
-    p_sfc = pressure[0]
+    p_sfc = surface_pressure
     terrain_correction = p_3 - (p_2 / (p_sfc - p_1))
 
-    # Convert all to 'dimensionless'
-    column_buoyancy_index *= units.dimensionless
-    mid_tropospheric_warming_index *= units.dimensionless
-    inversion_index *= units.dimensionless
-    terrain_correction *= units.dimensionless
-
     # Calculate G.D.I.
     return (column_buoyancy_index
             + mid_tropospheric_warming_index

tests/calc/test_thermo.py

@@ -2297,18 +2297,26 @@ def test_gdi():
                          -34.1, -37.3, -32.3, -34.1, -37.3, -41.1, -37.7, -58.1,
                          -57.5]) * units.degC
 
-    gdi = galvez_davison_index(pressure, temperature, dewpoint)
+    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
+    relative_humidity = relative_humidity_from_dewpoint(temperature, dewpoint)
+    mixrat = mixing_ratio_from_relative_humidity(pressure, temperature, relative_humidity)
     result = galvez_davison_index(
-        index_xarray_data_expanded.isobaric,
-        index_xarray_data_expanded.temperature,
-        index_xarray_data_expanded.dewpoint
+        pressure,
+        temperature,
+        mixrat,
+        pressure[0]
     )
 
     assert_array_almost_equal(
@@ -2328,10 +2336,16 @@ def test_gdi_no_950_raises_valueerror(index_xarray_data):
     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(
-            index_xarray_data.isobaric,
-            index_xarray_data.temperature,
-            index_xarray_data.dewpoint
+            pressure,
+            temperature,
+            mixrat,
+            pressure[0]
         )