Skip to content

Commit

Permalink
Major updates to calc_ctp and calc_humidity_index. First, functions r…
Browse files Browse the repository at this point in the history 
…e-worked to take a full Xarray object for NWP data, which mainly involved making `station_index` a positional rather than a kwarg. Added new `station_name` kwarg if needed for calc_ctp optional plotting. Avoids unnecessary interpolation of pressure when a specific pressure is requested because the user provided it so we don't need to interpolate. Similarly, the bottom/top pressures are also known in this case and don't need to be interpolated, only the data used at those pressures. Changes to debugging statements, bounds checking to see if the starting and top/bottom pressures are more than a certain pressure delta away from the values requested. Buildout of calc_humidity_index to be more like calc_ctp in many regards.
  • Loading branch information
@DanielAdriaansen
DanielAdriaansen committed Feb 6, 2025
1 parent c90b980 commit 7115d0d
Showing 1 changed file with 141 additions and 39 deletions.
180 changes: 141 additions & 39 deletions metcalcpy/diagnostics/land_surface.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,21 +9,22 @@
from pandas.core.series import Series
from xarray.core.dataarray import DataArray

def calc_ctp(pressure,temperature,start_pressure_hpa=-1,bot_pressure_hpa=100.0,top_pressure_hpa=300.0,interp=False,db=False,plotskewt=False,plotdir="",station=""):
def calc_ctp(pressure,temperature,station_index,start_pressure_hpa=-1,bot_pressure_hpa=100.0,top_pressure_hpa=300.0,interp=False,db=False,plotskewt=False,plotdir="",station_name=""):

""" Function for computing the Convective Triggering Potential
Args:
pressure (pint.Quantity): the vertical pressure profile
temperature (pint.Quantity): the vertical temperature profile
station_index (int): the integer index of the station currently being processed. Use -1 if a single station is being passed.
start_pressure_hpa (float, optional): the starting pressure to use. Default: -1 (bottom level in profile).
bot_pressure_hpa (float, optional): bottom pressure value of the layer, added to start_pressure_hpa. Default: 100 hPa.
top_pressure_hpa (float, optional): top pressure value of the layer, added to start_pressure_hpa. Default: 300 hPa.
interp (bool): Whether to interpolate data to exact pressures or use the closest. Default: False.
db (bool): Print debugging statements. Default: False
plotskewt (bool): Plot a Skew-T Log-P graphic of the CTP calculation. Default: False.
plotdir (string, optional): Directory where Skew-T plots should be written. Default: "".
station (string, optional): Location ID string used for labeling the Skew-T plot and image file name. Default: "".
station_name (string, optional): Location ID string used for labeling the Skew-T plot and image file name. Default: "".
Returns:
float32
Expand All @@ -35,6 +36,21 @@ def calc_ctp(pressure,temperature,start_pressure_hpa=-1,bot_pressure_hpa=100.0,t
Lorem Ipsum
"""

# If the station index is a non-negative value, then extract the profile at the station index
if station_index>=0:
temperature=temperature.isel(sid=station_index).values*units('degK')
pressure=pressure.isel(sid=station_index).values*units('hPa')

# Ensure there's no missing data values in the profile
#print(temperature)
#temperature = temperature[~np.isnan(temperature)]
#print(temperature)

# Subset the profile to only levels where pressure is > min_prs_profile
min_prs_profile = 100.0*units('hPa')
temperature = temperature[pressure>=min_prs_profile]
pressure = pressure[pressure>=min_prs_profile]

# Set a pressure difference for non-interpolation case
# If the closest pressure to the bot_pressure_hpa is more
Expand All @@ -50,27 +66,27 @@ def calc_ctp(pressure,temperature,start_pressure_hpa=-1,bot_pressure_hpa=100.0,t
print("USING LOWEST STARTING PRESSURE: %f\n" % (start_prs.m))
else:
if interp:
start_prs = log_interp_1d(start_pressure_hpa,pressure.m,pressure.m)
start_prs = start_prs*units('hPa')
# If the starting pressure is NaN, most likely the starting pressure was lower than
# all the pressures in the sounding.
if np.isnan(start_prs):
# If the requested starting pressure is greater than all the pressures in the
# profile, then we won't be able to interpolate to the requested starting pressure.
if start_pressure_hpa > np.max(pressure.m):
if db:
print("")
print("ERROR! REQUESTED STARTING PRESSURE INVALID")
print("UNABLE TO COMPUTE CTP.")
return(-9999.*units('J/kg'))
else:
start_prs = start_pressure_hpa*units('hPa')
if db:
print("")
print("USING INTERPOLATED STARTING PRESSURE: %f\n" % (start_prs.m))
print("USING ACTUAL REQUESTED STARTING PRESSURE: %f\n" % (start_prs.m))
else:
# Find the closest value. We'll just take the difference between the start pressure and pressure
# and find the index of the minimum
prs_diff = pressure-(start_pressure_hpa*units('hPa'))
start_prs = pressure[np.argmin(np.abs(prs_diff))]
if np.abs(start_pressure_hpa*units('hPa')-start_prs)>=max_prs_diff:
print("INFO: ACTUAL STARTING PRESSURE IS AT LEAST %3.2f hPa FROM REQUESTED START PRESSURE." % (max_prs_diff.m))
print("requested: start_pressure_hpa = %4.2f hPa" % (start_pressure_hpa.m))
print("requested: start_pressure_hpa = %4.2f hPa" % (start_pressure_hpa))
print("actual: start_pressure_hpa = %4.2f hPa" % (start_prs.m))
if db:
print("")
Expand All @@ -87,8 +103,14 @@ def calc_ctp(pressure,temperature,start_pressure_hpa=-1,bot_pressure_hpa=100.0,t

# Obtain information at the top and bottom of the layer
if interp:
prsBot, tmpBot = log_interp_1d(layer_bot_prs.m,pressure.m,pressure.m,temperature.m)
prsTop, tmpTop = log_interp_1d(layer_top_prs.m,pressure.m,pressure.m,temperature.m)

#prsBot, tmpBot = log_interp_1d(layer_bot_prs.m,pressure.m,pressure.m,temperature.m)
#prsTop, tmpTop = log_interp_1d(layer_top_prs.m,pressure.m,pressure.m,temperature.m)
tmpBot = log_interp_1d(layer_bot_prs.m,pressure.m,temperature.m)
tmpTop = log_interp_1d(layer_top_prs.m,pressure.m,temperature.m)
prsBot = np.array([layer_bot_prs.m])
prsTop = np.array([layer_top_prs.m])

if db:
print("")
print("USING INTERPOLATED LAYER BOTTOM PRESSURE: %f\n" % (prsBot))
Expand Down Expand Up @@ -137,10 +159,19 @@ def calc_ctp(pressure,temperature,start_pressure_hpa=-1,bot_pressure_hpa=100.0,t

# Find the index of the closest value.
# We'll just take the difference between the top/bottom pressure and find the index of the minimum
bot_diff = pressure-layer_bot_prs
top_diff = pressure-layer_top_prs
layer_bot_idx = np.argmin(np.abs(bot_diff))
layer_top_idx = np.argmin(np.abs(top_diff))
layer_bot_idx = np.argmin(np.abs(pressure-layer_bot_prs))
layer_top_idx = np.argmin(np.abs(pressure-layer_top_prs))

# Warn if the distance between the closest value to the bottom and top values exceeds max_prs_diff
if np.abs(pressure[layer_bot_idx]-layer_bot_prs)>=max_prs_diff:
print("INFO: ACTUAL BOTTOM PRESSURE IS AT LEAST %3.2f hPa FROM REQUESTED BOTTOM PRESSURE." % (max_prs_diff.m))
print("requested: layer_bot_prs = %4.2f hPa" % (layer_bot_prs.m))
print("actual: layer_bot_prs = %4.2f hPa" % (pressure[layer_bot_idx].m))
if np.abs(pressure[layer_top_idx]-layer_top_prs)>=max_prs_diff:
print("INFO: ACTUAL TOP PRESSURE IS AT LEAST %3.2f hPa FROM REQUESTED TOP PRESSURE." % (max_prs_diff.m))
print("requested: layer_top_prs = %4.2f hPa" % (layer_top_prs.m))
print("actual: layer_top_prs = %4.2f hPa" % (pressure[layer_top_idx].m))

if db:
print("")
print("INDEX OF LAYER BOT: %02d" % (int(layer_bot_idx)))
Expand Down Expand Up @@ -190,7 +221,7 @@ def calc_ctp(pressure,temperature,start_pressure_hpa=-1,bot_pressure_hpa=100.0,t
skew.ax.set_ylabel('Pressure (hPa)')
skew.ax.set_xlabel('Temperature (C)')
plt.title('CTP = %5.5f J/kg' % (float(CTP.m)),loc='left')
plt.title('STATION = %s' % (station))
plt.title('STATION = %s' % (station_name))
fig.savefig(os.path.join(plotdir,'CTP_%s.png' % (station)))
plt.close()

Expand Down Expand Up @@ -243,15 +274,19 @@ def calc_tci(soil_data,sfc_flux_data,skipna=True):
# Return the Terrestrial Coupling Index (TCI)
return covarTerm/soil_std

def calc_humidity_index(pressure,temperature,dewpoint,bot_pressure_hpa=950.0,top_pressure_hpa=850.0,interp=False):
def calc_humidity_index(pressure,temperature,dewpoint,station_index,start_pressure_hpa=-1,bot_pressure_hpa=50.0,top_pressure_hpa=150.0,interp=False,db=False):
""" Function for computing the Humidity Index
Args:
pressure (pint.Quantity): the vertical pressure profile
temperature (pint.Quantity): the vertical temperature profile
bot_pressure_hpa (float, optional): bottom pressure value of the layer. Default: 950 hPa.
top_pressure_hpa (float, optional): top pressure value of the layer. Default: 850 hPa.
dewpoint (pint.Quantity): the vertical dewpoint temperature profile
station_index (int): the integer index of the station currently being processed. Use -1 if a single station is being passed.
start_pressure_hpa (float, optional): the starting pressure to use. Default: -1 (bottom level in profile).
bot_pressure_hpa (float, optional): bottom pressure value of the layer, added to start_pressure_hpa. Default: 50 hPa.
top_pressure_hpa (float, optional): top pressure value of the layer, added to start_pressure_hpa. Default: 150 hPa.
interp (bool): perform vertical interpolation to bot_pressure_hpa and top_pressure_hpa or use closest. Default: False.
db (bool): Print debugging statements. Default: False
Returns:
float32
Expand All @@ -264,28 +299,90 @@ def calc_humidity_index(pressure,temperature,dewpoint,bot_pressure_hpa=950.0,top
"""

# If the station index is a non-negative value, then extract the profile at the station index
if station_index>=0:
temperature=temperature.isel(sid=station_index).values*units('degK')
pressure=pressure.isel(sid=station_index).values*units('hPa')
dewpoint=dewpoint.isel(sid=station_index).values*units('degK')

# Ensure there's no missing data values in the profile
#print(temperature)
#temperature = temperature[~np.isnan(temperature)]
#print(temperature)

# Subset the profile to only levels where pressure is > min_prs_profile
min_prs_profile = 100.0*units('hPa')
temperature = temperature[pressure>=min_prs_profile]
dewpoint = dewpoint[pressure>=min_prs_profile]
pressure = pressure[pressure>=min_prs_profile]

# Set a pressure difference for non-interpolation case
# If the closest pressure to the bot_pressure_hpa is more
# than this amount from the start_pressure_hpa, a warning will be
# printed.
max_prs_diff = 100.0*units('hPa')

# Set the bottom and top pressures for the layer
bot_pressure_hpa=bot_pressure_hpa*units('hPa')
top_pressure_hpa=top_pressure_hpa*units('hPa')
# Find the starting pressure in the profile
if start_pressure_hpa < 0:
start_prs = pressure[0]
if db:
print("")
print("USING LOWEST STARTING PRESSURE: %f\n" % (start_prs.m))
else:
if interp:
# If the requested starting pressure is greater than all the pressures in the
# profile, then we won't be able to interpolate to the requested starting pressure.
if start_pressure_hpa > np.max(pressure.m):
if db:
print("")
print("ERROR! REQUESTED STARTING PRESSURE INVALID")
print("UNABLE TO COMPUTE CTP.")
return(-9999.*units('J/kg'))
else:
start_prs = start_pressure_hpa*units('hPa')
if db:
print("")
print("USING ACTUAL REQUESTED STARTING PRESSURE: %f\n" % (start_prs.m))
else:
# Find the closest value. We'll just take the difference between the start pressure and pressure
# and find the index of the minimum
prs_diff = pressure-(start_pressure_hpa*units('hPa'))
start_prs = pressure[np.argmin(np.abs(prs_diff))]
if np.abs(start_pressure_hpa*units('hPa')-start_prs)>=max_prs_diff:
print("INFO: ACTUAL STARTING PRESSURE IS AT LEAST %3.2f hPa FROM REQUESTED START PRESSURE." % (max_prs_diff.m))
print("requested: start_pressure_hpa = %4.2f hPa" % (start_pressure_hpa))
print("actual: start_pressure_hpa = %4.2f hPa" % (start_prs.m))
if db:
print("")
print("USING NEAREST STARTING PRESSURE: %f\n" % (start_prs.m))

# Based on the starting pressure, set the initial layer bottom and top pressures
layer_bot_prs = start_prs-(bot_pressure_hpa*units('hPa'))
layer_top_prs = start_prs-(top_pressure_hpa*units('hPa'))

if db:
print("")
print("TARGET LAYER BOTTOM PRESSURE: %f\n" % (layer_bot_prs.m))
print("TARGET LAYER TOP PRESSURE: %f\n" % (layer_top_prs.m))

# Obtain information at the top and bottom of the layer
if interp:

# If the highest pressure in the sounding is < the bot_pressure_hpa, then skip this site
if bot_pressure_hpa.m > np.max(pressure.m):
# If the highest pressure in the sounding is < the layer_bot_prs, then skip this site
if layer_bot_prs.m > np.max(pressure.m):
print("ERROR! HIGHEST PRESSURE IN SOUNDING IS LOWER THAN REQUESTED BOTTOM PRESSURE.")
print("max pressure: %4.2f" % (np.max(pressure.m)))
print("requested bottom pressure: %4.2f" % (bot_pressure_hpa.m))
print("requested bottom pressure: %4.2f" % (layer_bot_prs.m))
print("UNABLE TO COMPUTE HI.")
return(-9999.*units('degK'))

if db:
print("")
print("INTERPOLATING TO BOTTOM PRESSURE: %f\n" % (layer_bot_prs.m))
print("INTERPOLATING TO TOP PRESSURE: %f\n" % (layer_top_prs.m))

tmpBot, dewBot = log_interp_1d(bot_pressure_hpa,pressure,temperature,dewpoint)
tmpTop, dewTop = log_interp_1d(top_pressure_hpa,pressure,temperature,dewpoint)
tmpBot, dewBot = log_interp_1d(layer_bot_prs,pressure,temperature,dewpoint)
tmpTop, dewTop = log_interp_1d(layer_top_prs,pressure,temperature,dewpoint)

# log_interp_1d returns array-like, so get the float values
tmpBot = tmpBot[0]
Expand All @@ -296,22 +393,27 @@ def calc_humidity_index(pressure,temperature,dewpoint,bot_pressure_hpa=950.0,top
else:

# Find the index of the closest pressure value to the bottom and top values requested
bot_idx = np.argmin(np.abs(pressure-bot_pressure_hpa))
top_idx = np.argmin(np.abs(pressure-top_pressure_hpa))
layer_bot_idx = np.argmin(np.abs(pressure-layer_bot_prs))
layer_top_idx = np.argmin(np.abs(pressure-layer_top_prs))

# Warn if the distance between the closest value to the bottom and top values exceeds max_prs_diff
if np.abs(pressure[bot_idx]-bot_pressure_hpa)>=max_prs_diff:
if np.abs(pressure[layer_bot_idx]-layer_bot_prs)>=max_prs_diff:
print("INFO: ACTUAL BOTTOM PRESSURE IS AT LEAST %3.2f hPa FROM REQUESTED BOTTOM PRESSURE." % (max_prs_diff.m))
print("requested: bot_pressure_hpa = %4.2f hPa" % (bot_pressure_hpa.m))
print("actual: bot_pressure_hpa = %4.2f hPa" % (pressure[bot_idx].m))
if np.abs(pressure[top_idx]-top_pressure_hpa)>=max_prs_diff:
print("requested: layer_bot_prs = %4.2f hPa" % (layer_bot_prs.m))
print("actual: layer_bot_prs = %4.2f hPa" % (pressure[layer_bot_idx].m))
if np.abs(pressure[layer_top_idx]-layer_top_prs)>=max_prs_diff:
print("INFO: ACTUAL TOP PRESSURE IS AT LEAST %3.2f hPa FROM REQUESTED TOP PRESSURE." % (max_prs_diff.m))
print("requested: top_pressure_hpa = %4.2f hPa" % (top_pressure_hpa.m))
print("actual: top_pressure_hps = %4.2f hPa" % (pressure[top_idx].m))
print("requested: layer_top_prs = %4.2f hPa" % (layer_top_prs.m))
print("actual: layer_top_prs = %4.2f hPa" % (pressure[layer_top_idx].m))

if db:
print("")
print("USING DATA AT NEAREST BOTTOM PRESSURE: %f\n" % (pressure[layer_bot_idx].m))
print("USING DATA AT NEAREST TOP PRESSURE: %f\n" % (pressure[layer_top_idx].m))

tmpBot = temperature[bot_idx]
dewBot = dewpoint[bot_idx]
tmpTop = temperature[top_idx]
dewTop = dewpoint[top_idx]
tmpBot = temperature[layer_bot_idx]
dewBot = dewpoint[layer_bot_idx]
tmpTop = temperature[layer_top_idx]
dewTop = dewpoint[layer_top_idx]

return (tmpBot-dewBot) + (tmpTop-dewTop)

0 comments on commit 7115d0d

Please sign in to comment.