================================ Forecast Model Vertical Sounding ================================ `Notebook `_ Python-AWIPS Tutorial Notebook -------------- Objectives ========== - Use python-awips to connect to an edex server - Request data using the `ModelSounding class `__ in addition to using the normal `DataAccess class `__ - Create and compare vertical sounding from different AWIPS model data with isobaric levels - Use `Shapely Point geometry `__ to define a point - Convert between units when necessary - Use MetPy to create `SkewT `__ and `Hodograph `__ plots -------------- Table of Contents ----------------- | `1 Imports `__\ | `2 EDEX Connection `__\ | `3 Define Useful Variables `__\ | `4 Function: get_surface_data() `__\ | `5 Function: get_levels_data() `__\ | `6 Function: plot_skewT() `__\ | `7 Retrieve Necessary Plotting Data `__\ | `8 Skew-T/Log-P `__\ | `9 Model Sounding Comparison `__\ | `10 See Also `__\ |     `10.1 Related Notebooks `__\ |     `10.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, ModelSounding import matplotlib.pyplot as plt import numpy as np from metpy.plots import SkewT, Hodograph from metpy.units import units from mpl_toolkits.axes_grid1.inset_locator import inset_axes from math import sqrt from shapely.geometry import Point `Top `__ -------------- 2 EDEX Connection ----------------- First we establish a connection to Unidata’s public EDEX server. This sets the proper server on the **DataAccessLayer**, which we will use numerous times throughout the notebook. .. code:: ipython3 server = 'edex-cloud.unidata.ucar.edu' DataAccessLayer.changeEDEXHost(server) `Top `__ -------------- 3 Define Useful Variables ------------------------- The plotting in this notebook needs a model name, a location point (defined by latitude and longitude), and the most recent time range with the initial forecast run. .. code:: ipython3 # Note the order is Lon,Lat and not Lat,Lon point = Point(-104.67,39.87) model="NAM40" # Get latest forecast cycle run timeReq = DataAccessLayer.newDataRequest("grid") timeReq.setLocationNames(model) cycles = DataAccessLayer.getAvailableTimes(timeReq, True) times = DataAccessLayer.getAvailableTimes(timeReq) fcstRun = DataAccessLayer.getForecastRun(cycles[-2], times) timeRange = [fcstRun[0]] print("Using " + model + " forecast time " + str(timeRange)) .. parsed-literal:: Using NAM40 forecast time [] `Top `__ -------------- 4 Function: get_surface_data() ------------------------------ This function is used to get the initial forecast model data for surface height. This is done separately from the rest of the heights to determine the surface pressure. It uses the `ModelSounding `__ data object from python-awips to retrieve all the relevant information. This function takes the model name, location, and time as attributes, and returns arrays for pressure, temperature, dewpoint, and the u and v wind components. .. code:: ipython3 def get_surface_data(modelName, location, time): """ model name, location, and timeRange desire """ # request data and sort response pressure,temp,dpt,ucomp,vcomp = [],[],[],[],[] use_parms = ['T','DpT','uW','vW','P'] use_level = "0.0FHAG" sndObject = ModelSounding.getSounding(modelName, use_parms, [use_level], location, time) if len(sndObject) > 0: for time in sndObject._dataDict: pressure.append(float(sndObject._dataDict[time][use_level]['P'])) temp.append(float(sndObject._dataDict[time][use_level]['T'])) dpt.append(float(sndObject._dataDict[time][use_level]['DpT'])) ucomp.append(float(sndObject._dataDict[time][use_level]['uW'])) vcomp.append(float(sndObject._dataDict[time][use_level]['vW'])) print("Found surface record at " + "%.1f" % pressure[0] + "MB") else: raise ValueError("sndObject returned empty for query [" + ', '.join(str(x) for x in (modelName, use_parms, point, use_level)) +"]") # return information for plotting return pressure,temp,dpt,ucomp,vcomp `Top `__ -------------- 5 Function: get_levels_data() ----------------------------- This function is similar to *get_surface_data()*, except it gets data values for presure heights above the surface. It uses the `ModelSounding `__ data object from python-awips to retrieve all the relevant information. It takes the model name, location, and time (similar to the other function), and also takes the instantiated pressure, temperature, dewpoint, and wind vector arrays. It returns the fully populated pressure, temperature, dewpoint, u-component, v-component, and computed wind arrays. .. code:: ipython3 def get_levels_data(modelName, location, time, pressure, temp, dpt, ucomp, vcomp): # Get isobaric levels with our requested parameters parms = ['T','DpT','uW','vW'] levelReq = DataAccessLayer.newDataRequest("grid", envelope=point) levelReq.setLocationNames(model) levelReq.setParameters(*(parms)) availableLevels = DataAccessLayer.getAvailableLevels(levelReq) # Clean levels list of unit string (MB, FHAG, etc.) levels = [] for lvl in availableLevels: name=str(lvl) if 'MB' in name and '_' not in name: # If this level is above (less than in mb) our 0.0FHAG record if float(name.replace('MB','')) < pressure[0]: levels.append(lvl) # Get Sounding sndObject = ModelSounding.getSounding(modelName, parms, levels, location, time) if not len(sndObject) > 0: raise ValueError("sndObject returned empty for query [" + ', '.join(str(x) for x in (model, parms, point, levels)) +"]") for time in sndObject._dataDict: for lvl in sndObject._dataDict[time].levels(): for parm in sndObject._dataDict[time][lvl].parameters(): if parm == "T": temp.append(float(sndObject._dataDict[time][lvl][parm])) elif parm == "DpT": dpt.append(float(sndObject._dataDict[time][lvl][parm])) elif parm == 'uW': ucomp.append(float(sndObject._dataDict[time][lvl][parm])) elif parm == 'vW': vcomp.append(float(sndObject._dataDict[time][lvl][parm])) else: print("WHAT IS THIS") print(sndObject._dataDict[time][lvl][parm]) # Pressure is our requested level rather than a returned parameter pressure.append(float(lvl.replace('MB',''))) # convert to numpy.array() pressure = np.array(pressure, dtype=float) temp = (np.array(temp, dtype=float) - 273.15) * units.degC dpt = (np.array(dpt, dtype=float) - 273.15) * units.degC ucomp = (np.array(ucomp, dtype=float) * units('m/s')).to('knots') vcomp = (np.array(vcomp, dtype=float) * units('m/s')).to('knots') wind = np.sqrt(ucomp**2 + vcomp**2) print("Using " + str(len(levels)) + " levels between " + str("%.1f" % max(pressure)) + " and " + str("%.1f" % min(pressure)) + "MB") return pressure,temp,dpt,ucomp,vcomp,wind `Top `__ -------------- 6 Function: plot_skewT() ------------------------ Since we’re plotting many different models for comparison, all that code was used to create this function. The function takes the model name, reference time, and the pressure, temperature, dewpoint, u-component, v-component, and wind arrays. It plots a skewT and hodograph using metpy. .. code:: ipython3 def plot_skewT(modelName, pressure, temp, dpt, ucomp, vcomp, wind, refTime): plt.rcParams['figure.figsize'] = (12, 14) # Skew-T skew = SkewT(rotation=45) skew.plot(pressure, temp, 'r', linewidth=2) skew.plot(pressure, dpt, 'g', linewidth=2) skew.plot_barbs(pressure, ucomp, vcomp) skew.plot_dry_adiabats() skew.plot_moist_adiabats() skew.plot_mixing_lines(linestyle=':') skew.ax.set_ylim(1000, np.min(pressure)) skew.ax.set_xlim(-50, 40) # Title plt.title(modelName + " (" + str(point) + ") " + str(refTime)) # Hodograph ax_hod = inset_axes(skew.ax, '40%', '40%', loc=2) h = Hodograph(ax_hod, component_range=max(wind.magnitude)) h.add_grid(increment=20) h.plot_colormapped(ucomp, vcomp, wind) # Dotted line at 0C isotherm l = skew.ax.axvline(0, color='c', linestyle='-', linewidth=1) plt.show() `Top `__ -------------- 7 Retrieve Necessary Plotting Data ---------------------------------- First we get the initial data at surface level using the get_surface_data function, and then pass those initial data arrays onto the get_levels_data request to finish populating for additional heights needed for Skew-T plots. We want to keep track of the pressure, temeperature, dewpoint, u-component, v-component, and wind arrays so we store them in variables to use later on. .. code:: ipython3 p,t,d,u,v = get_surface_data(model,point,timeRange) p,t,d,u,v,w = get_levels_data(model,point,timeRange,p,t,d,u,v) .. parsed-literal:: Found surface record at 833.2MB Using 32 levels between 833.2 and 50.0MB `Top `__ -------------- 8 Skew-T/Log-P -------------- Here we use our plot_skewT function to generate our skewT & hodograph charts for the data we retreived so far. This is where the pressure, temperature, dewpoint, and wind data is needed. .. code:: ipython3 plot_skewT(model, p, t, d, u, v, w, timeRange[0]) .. image:: Forecast_Model_Vertical_Sounding_files/Forecast_Model_Vertical_Sounding_25_0.png `Top `__ -------------- 9 Model Sounding Comparison --------------------------- Now that we know how to retreive and plot the data for one model, we can run a loop to retreive data for various models and plot them for comparison. In this example we’ll also plot RAP13 and GFS20 data to compare with NAM40. This is also where our functions become so important, because we can easily recall all that logic and keep this for-loop fairly simple. .. code:: ipython3 models = ["RAP13", "GFS20", "NAM40"] for modelName in models: timeReq = DataAccessLayer.newDataRequest("grid") timeReq.setLocationNames(modelName) cycles = DataAccessLayer.getAvailableTimes(timeReq, True) times = DataAccessLayer.getAvailableTimes(timeReq) fr = DataAccessLayer.getForecastRun(cycles[-1], times) print("Using " + modelName + " forecast time " + str(fr[0])) tr = [fr[0]] p,t,d,u,v = get_surface_data(modelName,point,tr) p,t,d,u,v,w = get_levels_data(modelName,point,tr,p,t,d,u,v) # Skew-T plot_skewT(modelName,p,t,d,u,v,w,tr[0]) .. parsed-literal:: Using RAP13 forecast time 2023-07-25 19:00:00 Found surface record at 839.4MB Using 32 levels between 839.4 and 100.0MB .. image:: Forecast_Model_Vertical_Sounding_files/Forecast_Model_Vertical_Sounding_28_1.png .. parsed-literal:: Using GFS20 forecast time 2023-07-25 12:00:00 Found surface record at 842.5MB Using 32 levels between 842.5 and 100.0MB .. image:: Forecast_Model_Vertical_Sounding_files/Forecast_Model_Vertical_Sounding_28_3.png .. parsed-literal:: Using NAM40 forecast time 2023-07-25 18:00:00 Found surface record at 833.8MB Using 32 levels between 833.8 and 50.0MB .. image:: Forecast_Model_Vertical_Sounding_files/Forecast_Model_Vertical_Sounding_28_5.png `Top `__ -------------- 10 See Also ----------- 10.1 Related Notebooks ~~~~~~~~~~~~~~~~~~~~~~ - `Grid Levels and Parameters `__ - `Upper Air BUFR Soundings `__ - `Model Sounding Data `__ 10.2 Additional Documentation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ **python-awips:** - `awips.DataAccessLayer `__ - `awips.ModelSounding `__ **matplotlib:** - `matplotlib.pyplot `__ **MetPy** - `metpy.wind_speed `__ - `metpy.wind_direction `__ - `metpy.vapor_pressure `__ - `metpy.dewpoint `__ - `metpy.skewt `__ - `metpy.hodograph `__ `Top `__ --------------