Sorting an array of a vertical temperature profile to determine location/magnitude of temperature inversions on a sounding.

[kylejgillett]

Hello GitHub!

I have been working with MetPy's sounding tools recently and have come pretty far in the development of a neat observed sounding plotter. I am working with a professor of mine to develop this code to be used for teaching observed soundings for a class. One feature we would like to add to our plot would be the location and depth of temperature inversions to aid students in identifying such inversions.

Given an array of temperatures pulled from the WyomingUpperAir.request_data siphon service, I would like to determine segments of the array where temperatures are increasing and compare them to the array of pressure values from WyomingUpperAir.request_data to find the pressure levels corresponding to the inversion (to find depth and location AGL). My issue right now is how to go about determining the index(s) where the temperature array is increasing (and then no longer increasing) and doing this for the entire array to catch multiple inversions.

Does there happen to be any MetPy support for this already? I have looked around and haven't found anything. Based on my research, finding increasing intervals of arrays seems to be a rather complex issue.

Does anyone have some insight as to how this could be accomplished? Thanks!

[DanielAdriaansen]

Hi thanks for posting! This was an interesting little problem, and I came up with one solution that appears to work but I only tested it on a single profile, so you may consider this as a "starting point" that may need some modifications for edge cases. It could also probably be improved/streamlined in several areas, but I didn't focus too much on elegance but rather functionality.

Maybe I'll try and work on contributing this, since I think this would be a handy capability. Maybe @dopplershift has some guidance, but I can get a issue/PR opened at the least. Needs a lot more testing though for sure.

Here's the code:

#/usr/bin/env python

# Sounding example:
# https://unidata.github.io/siphon/latest/examples/upperair/Wyoming_Request.html

# Imports
from datetime import datetime
from metpy.units import units
from siphon.simplewebservice.wyoming import WyomingUpperAir
import matplotlib.pyplot as plt
import numpy as np
from metpy.plots import SkewT

# Set the date and station
date = datetime(2022, 4, 15, 12)
station = 'LBF'

# Get the sounding data
df = WyomingUpperAir.request_data(date, station)

# Get temperature profile as a numpy ndarray
tprof = df['temperature']
print("\nSOUNDING TEMP:")
print(tprof)

# Take a difference of the temperature profile to find where the temperature is
# increasing/decreasing based on the sign, and use '.values' to get a numpy ndarray
# object rather than a Pandas object
tdiff = df['temperature'].diff().values
print("\nTEMPERATURE DIFF:")
print(tdiff)

# **** ASSUMPTION WARNING!!! *****
# NOTE: This section makes some assumptions that could be done differently, but this is how I thought of it
# We need to locate zeros. We will replace zero differences with the value before
# the zero difference, to assume that the behavior of the profile remains either
# decreasing or increasing, but that it cannot be constant. 
zerodiff = np.where(tdiff==0)[0]
print("\nLOCATIONS OF TDIFF==0:")
print(zerodiff)
tdiff[zerodiff] = tdiff[zerodiff-1]
print("\nTEMPERATURE DIFF WITHOUT ZEROS:")
print(tdiff)

# Now we take the sign of the difference, which just returns -1, 0, or 1
# for values <0 ,==0, and >0 respectively
tsign = np.sign(tdiff)
print("\nSIGN OF TEMPERATURE DIFFERENCE:")
print(tsign)

# Now we take a difference of the result of np.sign()
signdiff = np.diff(tsign)
print("DIFFERENCE OF TEMPERATURE DIFFERENCE SIGN:")
print(signdiff)

# The "signdiff" array should now tell us where inversion tops and bases are.
# The tops are where signdiff = -2, which is simply (-1 minus +1), Tinc to Tdec
# The bases are where signdiff = +2, which is simply (+1 minus -1), Tdec to Tinc
baseind = np.where(signdiff==2)[0]
print("\nINVERSION BASE INDEXES:")
print(baseind.tolist())

topind = np.where(signdiff==-2)[0]
print("\nINVERSION TOP INDEXES:")
print(topind.tolist())

# Now we can print characteristics about each inversion layer (low to high)
print("\nTOTAL NUMBER of INVERSIONS: %d\n" % (int(len(baseind))))
icnt = 0
while icnt < len(baseind):
  print("INVERSION %02d" % (int(icnt)))
  print("INVERSION LEVELS %02d" % (int(topind[icnt]-baseind[icnt])))
  print("INVERSION DELTA-P %3.2f hPa" % (df['pressure'].iloc[baseind[icnt]]-df['pressure'].iloc[topind[icnt]]))
  print("INVERSTION MAGNITUDE %3.2f C\n" % (df['temperature'].iloc[topind[icnt]]-df['temperature'].iloc[baseind[icnt]]))
  icnt += 1

# TOP data
ptop = df['pressure'].iloc[topind.tolist()]
ttop = df['temperature'].iloc[topind.tolist()]

# BASE data
pbase = df['pressure'].iloc[baseind.tolist()]
tbase = df['temperature'].iloc[baseind.tolist()]

# Plot the bases and tops as dots on the temperature profile
fig = plt.figure(1,figsize=(22,15))
skew = SkewT(fig=fig)

# Plot the Skew-T T/Td lines
skew.plot(df['pressure'],tprof,'r',linestyle='-',linewidth=1.5)
skew.plot(df['pressure'],df['dewpoint'],'g',linestyle='-',linewidth=1.5)

# Typical Skew-T lines (disabled to make it easier to show inversion lines)
#skew.plot_mixing_lines(pressure=np.arange(1000, 99, -20) * units.hPa,linestyle='dotted', color='tab:blue')
#skew.plot_moist_adiabats(t0=np.arange(233, 400, 5) * units.K,alpha=0.25, color='tab:green')
#skew.plot_dry_adiabats(t0=np.arange(233, 533, 10) * units.K,alpha=0.25, color='orangered')

# Plot inversion bases and tops as colored markers on the temperature profile
skew.plot(pbase,tbase,'k',marker='.',markersize=10,linestyle='')
skew.plot(ptop,ttop,'b',marker='.',markersize=10,linestyle='')

# Plot horizontal lines for inversion bases and tops instead
#[skew.plot([df['pressure'].iloc[x]]*40,np.arange(-100,100,5),'k',linestyle='--',linewidth=1) for x in baseind]
#[skew.plot([df['pressure'].iloc[x]]*40,np.arange(-100,100,5),'b',linestyle='--',linewidth=1) for x in topind]

#plt.show()
plt.savefig('inversions.png')

TOTAL NUMBER INVERSIONS: 8

INVERSION 00
INVERSION LEVELS 01
INVERSION DELTA-P 20.00 hPa
INVERSTION MAGNITUDE 2.40 C

INVERSION 01
INVERSION LEVELS 05
INVERSION DELTA-P 69.00 hPa
INVERSTION MAGNITUDE 9.90 C

INVERSION 02
INVERSION LEVELS 02
INVERSION DELTA-P 15.00 hPa
INVERSTION MAGNITUDE 1.00 C

INVERSION 03
INVERSION LEVELS 01
INVERSION DELTA-P 7.00 hPa
INVERSTION MAGNITUDE 0.80 C

INVERSION 04
INVERSION LEVELS 01
INVERSION DELTA-P 18.00 hPa
INVERSTION MAGNITUDE 1.40 C

INVERSION 05
INVERSION LEVELS 01
INVERSION DELTA-P 6.00 hPa
INVERSTION MAGNITUDE 0.80 C

INVERSION 06
INVERSION LEVELS 01
INVERSION DELTA-P 9.00 hPa
INVERSTION MAGNITUDE 1.00 C

INVERSION 07
INVERSION LEVELS 03
INVERSION DELTA-P 19.00 hPa
INVERSTION MAGNITUDE 3.80 C

[kylejgillett]

Here is the top and base output, it did miss a base. Interesting.

[DanielAdriaansen]

My replacement of zero differences assumed only two levels could have the same temperature, but certainly a layer could be isothermal for many levels. So I need to re-do that section. I will paste a fix here soon.

[DanielAdriaansen]

Here's a new code section:

# **** ASSUMPTION WARNING!!! *****
# NOTE: This section makes some assumptions that could be done differently, but this is how I thought of it
# We need to locate zeros. We will replace zero differences with the value before
# the zero difference, to assume that the behavior of the profile remains either
# decreasing or increasing, but that it cannot be constant.
zerodiff = np.where(tdiff==0)[0]
print("\nLOCATIONS OF TDIFF==0:")
print(zerodiff)
zcnt = 0
while zcnt < len(tdiff):
  # If this level is zero, overwrite it with the previous level
  if tdiff[zcnt]==0:
    tdiff[zcnt] = tdiff[zcnt-1]
  zcnt+=1
print("\nTEMPERATURE DIFF WITHOUT ZEROS:")
print(tdiff)

Handling the isothermal layers is still a bit tricky. For example, in this DTX profile, the temperature curve takes a bend from ~dry-adiabatic around 650 hPa to isothermal (temperature line cranks to the right), but the code doesn't indicate the beginning of the inversion until the top of the isothermal layer, rather than the bottom. But technically an inversion is not isothermal by definition (I don't think), so I think it's OK.

[kylejgillett]

Ah, interesting. Clever fix! Something I'd like to do with this is label each inversion (something simple like 'INV01') next to the inversion. For my other annotations, I use plt.annotate. I should be able to put something like this:

inversionLabel = "IVN%02d" % (int(icnt)) plt.annotate(inversionLabel, # this is the label text (tbase[icnt], pbase[icnt]), # these are the coordinates to position the label textcoords="offset points", xytext=(7,-10), fontsize=8, color='red', ha='center')

inside of a while loop like you used above (while icnt < len(baseind):) so that for every inversion, there is a label for it, right?

[kylejgillett]

Ope, I actually just got that to work!