From f5d36003f993e3565cd1179e6e210ce34f89f06c Mon Sep 17 00:00:00 2001 From: kgoebber Date: Thu, 21 Sep 2023 11:50:26 -0500 Subject: [PATCH] make changes to GDI calc --- .../calculations/Sounding_Calculations.py | 7 ++- src/metpy/calc/thermo.py | 62 +++++++------------ tests/calc/test_thermo.py | 28 ++++++--- 3 files changed, 51 insertions(+), 46 deletions(-) diff --git a/examples/calculations/Sounding_Calculations.py b/examples/calculations/Sounding_Calculations.py index 45dad18fcee..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,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) diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index fbd7fb2b092..c77a73906e5 100644 --- a/src/metpy/calc/thermo.py +++ b/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]) """ - # 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,12 +4495,10 @@ 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: @@ -4512,10 +4506,9 @@ def galvez_davison_index(pressure, temperature, dewpoint, vertical_dim=0): 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 diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py index aeab59a0e38..4cc55ec6d7d 100644 --- a/tests/calc/test_thermo.py +++ b/tests/calc/test_thermo.py @@ -2297,7 +2297,9 @@ 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) @@ -2305,10 +2307,16 @@ def test_gdi(): 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] )