======================== Upper Air BUFR Soundings ======================== `Notebook `_ 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 ----------------- | `1 Imports `__\ | `2 EDEX Connection `__\ |     `2.1 Initial EDEX Connection `__\ |     `2.2 Setting Additional Request Parameters `__\ |     `2.3 Available Location Names `__\ |     `2.4 Setting the Location Name `__\ | `3 Filtering by Time `__\ | `4 Get the Data! `__\ | `5 Use the Data! `__\ |     `5.1 Prepare Data Objects `__\ |     `5.2 Convert Units `__\ | `6 Plot the Data! `__\ | `7 See Also `__\ |     `7.1 Related Notebooks `__\ |     `7.2 Additional Documentation `__\ 1 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. .. code:: ipython3 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 units `Top `__ -------------- 2 EDEX Connection ----------------- 2.1 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. .. code:: ipython3 # Set the edex server DataAccessLayer.changeEDEXHost("edex-cloud.unidata.ucar.edu") request = DataAccessLayer.newDataRequest() # Set data type request.setDatatype("bufrua") 2.2 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. .. code:: ipython3 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) 2.3 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. **Note:** It is important to note the location names are listed by their WMO Station ID. Their corresponding location and site identifier can be looked up in `this table from UNdata `__. .. code:: ipython3 locations = DataAccessLayer.getAvailableLocationNames(request) locations.sort() print(locations) .. parsed-literal:: ['21824', '21946', '24266', '24343', '24641', '24688', '24959', '25123', '25703', '25913', '31004', '31088', '31300', '31369', '31510', '31538', '31770', '31873', '32061', '32098', '32150', '32389', '32477', '32540', '32618', '47122', '47138', '47158', '47401', '47412', '47582', '47646', '47678', '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', '71823', '71845', '71867', '71906', '71907', '71909', '71913', '71917', '71924', '71925', '71926', '71934', '71945', '71957', '71964', '72201', '72202', '72206', '72208', '72210', '72214', '72215', '72221', '72230', '72233', '72235', '72240', '72248', '72249', '72250', '72251', '72261', '72265', '72274', '72293', '72305', '72317', '72318', '72327', '72340', '72357', '72363', '72364', '72365', '72376', '72381', '72388', '72393', '72402', '72403', '72426', '72440', '72451', '72456', '72469', '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', '74794', '78016', '78384', '78397', '78486', '78526', '78583', '78866', '78954', '78970', '78988', '80001', '91165', '91212', '91285', '91334', '91348', '91366', '91376', '91408', '91413', '91610', '91643', '91680', '91765', '94120', '94203', '94299', '94332', '94461', '94510', '94578', '94637', '94638', '94653', '94659', '94672', '94711', '94776', '94996'] 2.4 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. .. code:: ipython3 # Set station ID (not name) request.setLocationNames("72562") #KLBF `Top `__ -------------- 3 Filtering 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. .. code:: ipython3 # Get all times datatimes = DataAccessLayer.getAvailableTimes(request) `Top `__ -------------- 4 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. .. code:: ipython3 # 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())) .. parsed-literal:: parms = ['tpSigT', 'prSigT', 'tdSigT'] site = 72562 geom = POINT (-100.7005615234375 41.14971923828125) datetime = 2023-05-25 12:00:00 reftime = May 25 23 12:00:00 GMT fcstHour = 0 period = (May 25 23 12:00:00 , May 25 23 12:00:00 ) `Top `__ -------------- 5 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. 5.1 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. .. code:: ipython3 # 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.nan 5.2 Convert 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. .. code:: ipython3 # 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.degC `Top `__ -------------- 6 Plot the Data! ---------------- Create and display SkewT and Hodograph plots using MetPy. .. code:: ipython3 # 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() .. image:: Upper_Air_BUFR_Soundings_files/Upper_Air_BUFR_Soundings_28_0.png `Top `__ -------------- 7 See Also ---------- 7.1 Related Notebooks ~~~~~~~~~~~~~~~~~~~~~ - `Grid Levels and Parameters `__ - `Model Sounding Data `__ - `Forecast Model Vertical Sounding `__ 7.2 Additional Documentation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ **python-awips:** - `awips.DataAccessLayer `__ - `awips.PyGeometryData `__ **matplotlib:** - `matplotlib.pyplot `__ **MetPy** - `metpy.wind_components `__ - `metpy.lcl `__ (Lifted Condensation Level) - `metpy.parcel_profile `__ - `metpy.skewt `__ - `metpy.hodograph `__ `Top `__ --------------