Skew-T’s
Contents
Skew-T’s¶
Requesting Upper Air Data from Siphon and Plotting a Skew-T with MetPy¶
Based on Unidata’s MetPy Monday # 15 -17 YouTube links:
https://www.youtube.com/watch?v=OUTBiXEuDIU;
https://www.youtube.com/watch?v=oog6_b-844Q;
https://www.youtube.com/watch?v=b0RsN9mCY5k¶
Import the libraries we need¶
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import metpy.calc as mpcalc
from metpy.plots import Hodograph, SkewT
from metpy.units import units, pandas_dataframe_to_unit_arrays
from datetime import datetime, timedelta
import pandas as pd
from siphon.simplewebservice.wyoming import WyomingUpperAir
Next we use methods from the datetime library to determine the latest time, and then use it to select the most recent date and time that we’ll use for our selected Skew-T.¶
now = datetime.utcnow()
curr_year = now.year
curr_month = now.month
curr_day = now.day
curr_hour = now.hour
print(f"Current time is: {now}")
print(f"Current year, month, date, hour: {curr_year}, {curr_month}, {curr_day}, {curr_hour}")
if (curr_hour > 13) :
raob_hour = 12
hr_delta = curr_hour - 12
elif (curr_hour > 1):
raob_hour = 0
hr_delta = curr_hour - 0
else:
raob_hour = 12
hr_delta = curr_hour + 12
raob_time = now - timedelta(hours=hr_delta)
raob_year = raob_time.year
raob_month = raob_time.month
raob_day = raob_time.day
raob_hour = raob_time.hour
print(f"Time of RAOB is: {raob_year}, {raob_month}, {raob_day}, {raob_hour}")
Current time is: 2024-04-30 00:59:22.308818
Current year, month, date, hour: 2024, 4, 30, 0
Time of RAOB is: 2024, 4, 29, 12
Construct a datetime object that we will use in our query to the data server. Note what it looks like.¶
query_date = datetime(raob_year,raob_month,raob_day,raob_hour)
query_date
datetime.datetime(2024, 4, 29, 12, 0)
If desired, we can choose a past date and time.¶
#current = True
current = False
if (current == False):
query_date = datetime(2012, 10, 29, 12)
raob_timeStr = query_date.strftime("%Y-%m-%d %H UTC")
raob_timeTitle = query_date.strftime("%H00 UTC %-d %b %Y")
print(raob_timeStr)
print(raob_timeTitle)
2012-10-29 12 UTC
1200 UTC 29 Oct 2012
Select our station and query the data server.¶
station = 'ALB'
df = WyomingUpperAir.request_data (query_date,station)
Possible data retrieval errors:
- If you get a
No data available for YYYY-MM-DD HHZ for station XXX, you may have specified an incorrect date or location, or there just might not have been a successful RAOB at that date and time. - If you get a
503 Server Error, try requesting the sounding again. The Wyoming server does not handle simultaneous requests well. Usually the request will eventually work.
What does the returned data file look like? Well, it looks like a Pandas Dataframe!¶
df
| pressure | height | temperature | dewpoint | direction | speed | u_wind | v_wind | station | station_number | time | latitude | longitude | elevation | pw | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 993.0 | 96 | 12.6 | 7.6 | 0.0 | 8.0 | -0.000000 | -8.000000e+00 | ALB | 72518 | 2012-10-29 12:00:00 | 42.7 | -73.83 | 96.0 | 16.8 |
| 1 | 990.0 | 122 | 12.0 | 7.0 | 2.0 | 10.0 | -0.348995 | -9.993908e+00 | ALB | 72518 | 2012-10-29 12:00:00 | 42.7 | -73.83 | 96.0 | 16.8 |
| 2 | 968.6 | 305 | 10.4 | 6.5 | 20.0 | 25.0 | -8.550504 | -2.349232e+01 | ALB | 72518 | 2012-10-29 12:00:00 | 42.7 | -73.83 | 96.0 | 16.8 |
| 3 | 955.0 | 423 | 9.4 | 6.1 | 26.0 | 26.0 | -11.397650 | -2.336865e+01 | ALB | 72518 | 2012-10-29 12:00:00 | 42.7 | -73.83 | 96.0 | 16.8 |
| 4 | 933.8 | 610 | 8.1 | 5.7 | 35.0 | 27.0 | -15.486564 | -2.211711e+01 | ALB | 72518 | 2012-10-29 12:00:00 | 42.7 | -73.83 | 96.0 | 16.8 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 138 | 10.4 | 30480 | -57.3 | NaN | 295.0 | 22.0 | 19.938771 | -9.297602e+00 | ALB | 72518 | 2012-10-29 12:00:00 | 42.7 | -73.83 | 96.0 | 16.8 |
| 139 | 10.0 | 30720 | -58.1 | NaN | 280.0 | 23.0 | 22.650578 | -3.993908e+00 | ALB | 72518 | 2012-10-29 12:00:00 | 42.7 | -73.83 | 96.0 | 16.8 |
| 140 | 9.7 | 30911 | -58.9 | NaN | 277.0 | 26.0 | 25.806200 | -3.168603e+00 | ALB | 72518 | 2012-10-29 12:00:00 | 42.7 | -73.83 | 96.0 | 16.8 |
| 141 | 9.0 | 31394 | -58.2 | NaN | 270.0 | 34.0 | 34.000000 | 6.245699e-15 | ALB | 72518 | 2012-10-29 12:00:00 | 42.7 | -73.83 | 96.0 | 16.8 |
| 142 | 8.7 | 31595 | -57.9 | NaN | NaN | NaN | NaN | NaN | ALB | 72518 | 2012-10-29 12:00:00 | 42.7 | -73.83 | 96.0 | 16.8 |
143 rows × 15 columns
Examine one of the columns
df['pressure']
0 993.0
1 990.0
2 968.6
3 955.0
4 933.8
...
138 10.4
139 10.0
140 9.7
141 9.0
142 8.7
Name: pressure, Length: 143, dtype: float64
The University of Wyoming sounding reader includes a units attribute, which produces a Python dictionary of units specific to Wyoming sounding data. What does it look like?¶
df.units
{'pressure': 'hPa',
'height': 'meter',
'temperature': 'degC',
'dewpoint': 'degC',
'direction': 'degrees',
'speed': 'knot',
'u_wind': 'knot',
'v_wind': 'knot',
'station': None,
'station_number': None,
'time': None,
'latitude': 'degrees',
'longitude': 'degrees',
'elevation': 'meter',
'pw': 'millimeter'}
Attach units to the various columns using MetPy’s pandas_dataframe_to_unit_arrays function, as follows:¶
df = pandas_dataframe_to_unit_arrays(df)
Now that units are attached, how does the pressure column look?
df['pressure']
| Magnitude | [993.0 990.0 968.6 955.0 933.8 925.0 900.0 876.0 867.1 850.0 833.0 818.0 |
|---|---|
| Units | hectopascal |
As with any Pandas Dataframe, we can select columns, so we do so now for all the relevant variables for our Skew-T.¶
P = df['pressure']
Z = df['height']
T = df['temperature']
Td = df['dewpoint']
u = df['u_wind']
v = df['v_wind']
Now that u and v are units aware, we can use one of MetPy’s diagnostic routines to calculate the wind speed from the vector’s components.
wind_speed = mpcalc.wind_speed(u,v)
wind_speed
| Magnitude | [8.0 10.0 25.0 26.0 27.0 27.999999999999996 29.0 33.0 34.0 36.0 39.0 |
|---|---|
| Units | knot |
Thermodynamic Calculations¶
Often times we will want to calculate some thermodynamic parameters of a sounding. The MetPy calc module has many such calculations already implemented!
Lifting Condensation Level (LCL) - The level at which an air parcel’s relative humidity becomes 100% when lifted along a dry adiabatic path.
Parcel Path - Path followed by a hypothetical parcel of air, beginning at the surface temperature/pressure and rising dry adiabatically until reaching the LCL, then rising moist adiabatially.
Calculate the LCL. Note that the units will appear at the end when we print them out!¶
lclp, lclt = mpcalc.lcl(P[0], T[0], Td[0])
print (f'LCL Pressure = ,{lclp}, LCL Temperature = ,{lclt}')
LCL Pressure = ,920.7394732770757 hectopascal, LCL Temperature = ,6.497694364794597 degree_Celsius
Calculate the parcel profile. Here, we pass in the array of all pressure values, and specify the first (index 0, remember Python’s zero-based indexing) index of the temperature and dewpoint arrays as the starting point of the parcel path.¶
parcel_prof = mpcalc.parcel_profile(P, T[0], Td[0])
/tmp/ipykernel_1786263/3676909138.py:1: UserWarning: Duplicate pressure(s) [31.7] hPa provided. Output profile includes duplicate temperatures as a result.
parcel_prof = mpcalc.parcel_profile(P, T[0], Td[0])
What does the parcel profile look like? Well, it’s just an array of all the temperatures of the parcel at each of the pressure levels as it ascends.¶
parcel_prof
| Magnitude | [285.75 285.5030782457374 283.7260147751601 282.5820425453529 |
|---|---|
| Units | kelvin |
Set the limits for the x- and y-axes so they can be changed and easily re-used when desired.¶
pTop = 100
pBot = 1050
tMin = -60
tMax = 40
Since we are plotting a skew-T, define the log-p axis.¶
logTop = np.log10(pTop)
logBot = np.log10(pBot)
interval = np.logspace(logTop,logBot) * units('hPa')
This interval will tend to oversample points near the top of the atmosphere relative to the bottom. One effect of this will be that wind barbs will be hard to read, especially at higher levels. Use a resampling method from NumPy to deal with this.¶
idx = mpcalc.resample_nn_1d(P, interval)
Now we are ready to create our skew-T. Plot it and save it as a PNG file.¶
Note:
Annoyingly, when specifying the color of the dry/moist adiabats, and mixing ratio lines, the parameter must be colors, not color. In other Matplotlib-related code lines in the cell below where we set the color, such as skew.ax.axvline , color must be singular.
# Create a new figure. The dimensions here give a good aspect ratio
fig = plt.figure(figsize=(9, 9))
skew = SkewT(fig, rotation=30)
# Plot the data using normal plotting functions, in this case using
# log scaling in Y, as dictated by the typical meteorological plot
skew.plot(P, T, 'r', linewidth=3)
skew.plot(P, Td, 'g', linewidth=3)
skew.plot_barbs(P[idx], u[idx], v[idx])
skew.ax.set_ylim(pBot, pTop)
skew.ax.set_xlim(tMin, tMax)
# Plot LCL as black dot
skew.plot(lclp, lclt, 'ko', markerfacecolor='black')
# Plot the parcel profile as a black line
skew.plot(P, parcel_prof, 'k', linestyle='--', linewidth=2)
# Shade areas of CAPE and CIN
#skew.shade_cin(p, T, parcel_prof)
#skew.shade_cape(p, T, parcel_prof)
# Shade areas between a certain temperature (in this case, the DGZ)
skew.shade_area(P.m,-18,-12,color=['purple'])
# Plot a zero degree isotherm
skew.ax.axvline(0, color='c', linestyle='-', linewidth=2)
# Add the relevant special lines
skew.plot_dry_adiabats()
skew.plot_moist_adiabats(colors='#58a358', linestyle='-')
skew.plot_mixing_lines()
plt.title(f"RAOB for {station} valid at: {raob_timeTitle}")
# Create a hodograph
# Create an inset axes object that is 30% width and height of the
# figure and put it in the upper right hand corner.
ax_hod = inset_axes(skew.ax, '30%', '30%', loc=1)
h = Hodograph(ax_hod, component_range=100.)
h.add_grid(increment=30)
h.plot_colormapped(u, v, df['height'])
plt.savefig (f"station{raob_timeStr}_skewt.png")
Here’s how the plot would appear if we hadn’t done the resampling. Note that the wind barbs are plotted at every level, and are harder to read:
# Create a new figure. The dimensions here give a good aspect ratio
fig = plt.figure(figsize=(9, 9))
skew = SkewT(fig, rotation=30)
# Plot the data using normal plotting functions, in this case using
# log scaling in Y, as dictated by the typical meteorological plot
skew.plot(P, T, 'r', linewidth=3)
skew.plot(P, Td, 'g', linewidth=3)
skew.plot_barbs(P, u, v)
skew.ax.set_ylim(pBot, pTop)
skew.ax.set_xlim(tMin, tMax)
# Plot LCL as black dot
skew.plot(lclp, lclt, 'ko', markerfacecolor='black')
# Plot the parcel profile as a black line
skew.plot(P, parcel_prof, 'k', linestyle='--', linewidth=2)
# Shade areas of CAPE and CIN
skew.shade_cin(P, T, parcel_prof)
skew.shade_cape(P, T, parcel_prof)
# Plot a zero degree isotherm
skew.ax.axvline(0, color='c', linestyle='-', linewidth=2)
# Add the relevant special lines
skew.plot_dry_adiabats()
skew.plot_moist_adiabats(colors='#58a358', linestyle='-')
skew.plot_mixing_lines()
plt.title(f"RAOB for {station} valid at: {raob_timeTitle}")
# Create a hodograph
# Create an inset axes object that is 30% width and height of the
# figure and put it in the upper right hand corner.
ax_hod = inset_axes(skew.ax, '30%', '30%', loc=1)
h = Hodograph(ax_hod, component_range=100.)
h.add_grid(increment=30)
h.plot_colormapped(u, v, df['height']) # Plot a line colored by a parameter
<matplotlib.collections.LineCollection at 0x146528d62350>
What’s next:
Try different RAOB locations and dates (for example, draw a skew-t relevant to your class project).
Here is a link to a map of NWS RAOB sites.
Depending on the situation, you may want to omit plotting the DGZ, LCL, parcel path, CAPE/CIN, freezing line, and/or move the location of the hodograph.