Upper Air BUFR Soundings¶
Python-AWIPS Tutorial Notebook
Objectives¶
Retrieve an Upper Air vertical profile from EDEX
Plot a Skew-T/Log-P chart with Matplotlib and MetPy
Understand the bufrua plugin returns separate objects for parameters at mandatory levels and at significant temperature levels
Significant temperature levels are used to plot the pressure, temperature and dewpoint lines
Mandatory levels are used to plot the wind profile
Table of Contents¶
Imports¶
The imports below are used throughout the notebook. Note the first import is coming directly from python-awips and allows us to connect to an EDEX server. The subsequent imports are for data manipulation and visualization.
from awips.dataaccess import DataAccessLayer
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import numpy as np
from metpy.calc import wind_components, lcl, parcel_profile
from metpy.plots import SkewT, Hodograph
from metpy.units import unitsEDEX Connection¶
Initial EDEX Connection¶
First we establish a connection to Unidata’s public EDEX server. With that connection made, we can create a new data request object and set the data type to bufrua, and define additional parameters and an identifier on the request.
# Set the edex server
DataAccessLayer.changeEDEXHost("edex-cloud.unidata.ucar.edu")
request = DataAccessLayer.newDataRequest()
# Set data type
request.setDatatype("bufrua")Setting Additional Request Parameters¶
Here we populate arrays of all the parameters that will be necessary for plotting the Skew-T. The MAN_PARAMS are the mandatory levels and the SIGT_PARAMS are the significant temperature parameters that were both mentioned in the objectives section above.
Also request the station name and elevation to use in the figure title later on.
MAN_PARAMS = set(['prMan', 'wdMan', 'wsMan'])
SIGT_PARAMS = set(['prSigT', 'tpSigT', 'tdSigT'])
request.setParameters("staElev", "staName")
request.getParameters().extend(MAN_PARAMS)
request.getParameters().extend(SIGT_PARAMS)Available Location Names¶
When working with a new data type, it is often useful to investigate all available options for a particular setting. Shown below is how to see all available location names for a data request with type bufrua. This step is not necessary if you already know exactly what the location ID you’re interested in is.
locations = DataAccessLayer.getAvailableLocationNames(request)
locations.sort()
print(locations)['21824', '21946', '24266', '24343', '24641', '24688', '24959', '25123', '25703', '25913', '31004', '31088', '31300', '31369', '31510', '31770', '31873', '32061', '32098', '32389', '32477', '32540', '32618', '47138', '47401', '47412', '47582', '47600', '47646', '47678', '47778', '47807', '47827', '47909', '47918', '47945', '47971', '47991', '70026', '70133', '70200', '70219', '70231', '70261', '70273', '70308', '70316', '70326', '70350', '70361', '70398', '70414', '71043', '71081', '71082', '71109', '71119', '71603', '71722', '71802', '71811', '71815', '71816', '71845', '71867', '71906', '71907', '71909', '71913', '71917', '71924', '71925', '71926', '71934', '71945', '71957', '71964', '72201', '72202', '72206', '72208', '72210', '72215', '72221', '72230', '72233', '72235', '72240', '72248', '72249', '72250', '72251', '72261', '72265', '72274', '72293', '72305', '72317', '72318', '72327', '72340', '72357', '72363', '72364', '72365', '72376', '72388', '72393', '72403', '72426', '72440', '72451', '72456', '72476', '72489', '72493', '72501', '72518', '72520', '72528', '72558', '72562', '72572', '72582', '72597', '72632', '72634', '72645', '72649', '72659', '72662', '72672', '72681', '72694', '72712', '72747', '72764', '72768', '72776', '72786', '72797', '74004', '74005', '74389', '74455', '74560', '74626', '74794', '76225', '76256', '76394', '76405', '76458', '76526', '76612', '76644', '76654', '76679', '76692', '76805', '78016', '78384', '78397', '78526', '78866', '78954', '78970', '78988', '80001', '91165', '91212', '91285', '91334', '91348', '91376', '91408', '91413', '91680', '91765', '94203', '94299', '94326', '94332', '94461', '94510', '94578', '94637', '94653', '94659', '94672', '94711', '94995', '94996']
Setting the Location Name¶
In this case we’re setting the location name to the ID for KLBF which is the North Platte Regional Airport/Lee Bird, Field in Nebraska.
# Set station ID (not name)
request.setLocationNames("72562") #KLBFFiltering by Time¶
Models produce many different time variants during their runs, so let’s limit the data to the most recent time and forecast run.
# Get all times
datatimes = DataAccessLayer.getAvailableTimes(request)Get the Data!¶
Here we can now request our data response from the EDEX server with our defined time filter.
Printing out some data from the first object in the response array can help verify we received the data we were interested in.
# Get most recent record
response = DataAccessLayer.getGeometryData(request,times=datatimes[-1].validPeriod)
obj = response[0]
print("parms = " + str(obj.getParameters()))
print("site = " + str(obj.getLocationName()))
print("geom = " + str(obj.getGeometry()))
print("datetime = " + str(obj.getDataTime()))
print("reftime = " + str(obj.getDataTime().getRefTime()))
print("fcstHour = " + str(obj.getDataTime().getFcstTime()))
print("period = " + str(obj.getDataTime().getValidPeriod()))parms = ['wsMan', 'wdMan', 'prMan']
site = 72562
geom = POINT (-100.7005615234375 41.14971923828125)
datetime = 2026-01-23 18:00:00 (0)
reftime = Jan 23 26 18:00:00 GMT
fcstHour = 0
period = (Jan 23 26 18:00:00 , Jan 23 26 18:00:00 )
Use the Data!¶
Since we filtered on time, and requested the data in the previous cell, we now have a response object we can work with.
Prepare Data Objects¶
Here we construct arrays for each parameter to plot (temperature, dewpoint, pressure, and wind components). After populating each of the arrays, we sort and mask missing data.
# Initialize data arrays
prMan,wdMan,wsMan = np.array([]),np.array([]),np.array([])
prSig,tpSig,tdSig = np.array([]),np.array([]),np.array([])
manGeos = []
sigtGeos = []
# Build arrays
for ob in response:
parm_array = ob.getParameters()
if set(parm_array) & MAN_PARAMS:
manGeos.append(ob)
prMan = np.append(prMan,ob.getNumber("prMan"))
wdMan = np.append(wdMan,ob.getNumber("wdMan"))
wsMan, wsUnit = np.append(wsMan,ob.getNumber("wsMan")), ob.getUnit("wsMan")
continue
if set(parm_array) & SIGT_PARAMS:
sigtGeos.append(ob)
prSig = np.append(prSig,ob.getNumber("prSigT"))
tpSig = np.append(tpSig,ob.getNumber("tpSigT"))
tpUnit = ob.getUnit("tpSigT")
tdSig = np.append(tdSig,ob.getNumber("tdSigT"))
continue
# Sort mandatory levels (but not sigT levels) because of the 1000.MB interpolation inclusion
ps = prMan.argsort()[::-1]
wpres = prMan[ps]
direc = wdMan[ps]
spd = wsMan[ps]
# Flag missing data
prSig[prSig <= -9999] = np.nan
tpSig[tpSig <= -9999] = np.nan
tdSig[tdSig <= -9999] = np.nan
wpres[wpres <= -9999] = np.nan
direc[direc <= -9999] = np.nan
spd[spd <= -9999] = np.nanConvert Units¶
We need to modify the units several of the data parameters are returned in. Here we convert the units for Temperature and Dewpoint from Kelvin to Celsius, convert pressure to milibars, and extract wind for both the u and v directional components in Knots and Radians.
# assign units
p = (prSig/100) * units.mbar
wpres = (wpres/100) * units.mbar
u,v = wind_components(spd * units.knots, np.deg2rad(direc))
if tpUnit == 'K':
T = (tpSig-273.15) * units.degC
Td = (tdSig-273.15) * units.degCPlot the Data!¶
Create and display SkewT and Hodograph plots using MetPy.
# Create SkewT/LogP
plt.rcParams['figure.figsize'] = (10, 12)
skew = SkewT()
skew.plot(p, T, 'r', linewidth=2)
skew.plot(p, Td, 'g', linewidth=2)
skew.plot_barbs(wpres, u, v)
skew.ax.set_ylim(1000, 100)
skew.ax.set_xlim(-60, 30)
title_string = " T(F) Td "
title_string += " " + str(ob.getString("staName"))
title_string += " " + str(ob.getDataTime().getRefTime())
title_string += " (" + str(ob.getNumber("staElev")) + "m elev)"
title_string += "\n" + str(round(T[0].to('degF').item(),1))
title_string += " " + str(round(Td[0].to('degF').item(),1))
plt.title(title_string, loc='left')
# Calculate LCL height and plot as black dot
lcl_pressure, lcl_temperature = lcl(p[0], T[0], Td[0])
skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')
# Calculate full parcel profile and add to plot as black line
prof = parcel_profile(p, T[0], Td[0]).to('degC')
skew.plot(p, prof, 'k', linewidth=2)
# An example of a slanted line at constant T -- in this case the 0 isotherm
l = skew.ax.axvline(0, color='c', linestyle='--', linewidth=2)
# Draw hodograph
ax_hod = inset_axes(skew.ax, '30%', '30%', loc=3)
h = Hodograph(ax_hod, component_range=max(wsMan))
h.add_grid(increment=20)
h.plot_colormapped(u, v, spd)
# Show the plot
plt.show()