python-awips/_sources/examples/generated/Forecast_Model_Vertical_Sounding.rst.txt

================================
Forecast Model Vertical Sounding
================================
`Notebook <http://nbviewer.ipython.org/github/Unidata/python-awips/blob/master/examples/notebooks/Forecast_Model_Vertical_Sounding.ipynb>`_
Python-AWIPS Tutorial Notebook

--------------

Objectives
==========

-  Use python-awips to connect to an edex server
-  Request data using the `ModelSounding
   class <http://unidata.github.io/python-awips/api/ModelSounding.html>`__
   in addition to using the normal `DataAccess
   class <http://unidata.github.io/python-awips/api/DataAccessLayer.html>`__
-  Create and compare vertical sounding from different AWIPS model data
   with isobaric levels
-  Use `Shapely Point
   geometry <https://shapely.readthedocs.io/en/stable/reference/shapely.Point.html>`__
   to define a point
-  Convert between units when necessary
-  Use MetPy to create
   `SkewT <https://unidata.github.io/MetPy/latest/api/generated/metpy.plots.SkewT.html>`__
   and
   `Hodograph <https://unidata.github.io/MetPy/latest/api/generated/metpy.plots.Hodograph.html>`__
   plots

--------------

Table of Contents
-----------------

| `1
  Imports <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html#imports>`__\ 
| `2 EDEX
  Connection <http://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html#edex-connection>`__\ 
| `3 Define Useful
  Variables <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html#define-useful-variables>`__\ 
| `4 Function:
  get_surface_data() <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html#function-get-surface-data>`__\ 
| `5 Function:
  get_levels_data() <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html#function-get-levels-data>`__\ 
| `6 Function:
  plot_skewT() <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html#function-plot-skewt>`__\ 
| `7 Retrieve Necessary Plotting
  Data <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html#retrieve-necessary-plotting-data>`__\ 
| `8
  Skew-T/Log-P <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html#skew-t-log-p>`__\ 
| `9 Model Sounding
  Comparison <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html#model-sounding-comparison>`__\ 
| `10 See
  Also <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html#see-also>`__\ 
|     `10.1 Related
  Notebooks <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html#related-notebooks>`__\ 
|     `10.2 Additional
  Documentation <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html#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 <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html>`__

--------------

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 <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html>`__

--------------

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 [<DataTime instance: 2023-07-25 12:00:00 >]


`Top <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html>`__

--------------

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 <http://unidata.github.io/python-awips/api/ModelSounding.html>`__
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 <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html>`__

--------------

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 <http://unidata.github.io/python-awips/api/ModelSounding.html>`__
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 <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html>`__

--------------

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 <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html>`__

--------------

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 <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html>`__

--------------

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 <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html>`__

--------------

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 <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html>`__

--------------

10 See Also
-----------

10.1 Related Notebooks
~~~~~~~~~~~~~~~~~~~~~~

-  `Grid Levels and
   Parameters <https://unidata.github.io/python-awips/examples/generated/Grid_Levels_and_Parameters.html>`__
-  `Upper Air BUFR
   Soundings <http://unidata.github.io/python-awips/examples/generated/Upper_Air_BUFR_Soundings.html>`__
-  `Model Sounding
   Data <http://unidata.github.io/python-awips/examples/generated/Model_Sounding_Data.html>`__

10.2 Additional Documentation
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

**python-awips:**

-  `awips.DataAccessLayer <http://unidata.github.io/python-awips/api/DataAccessLayer.html>`__
-  `awips.ModelSounding <https://unidata.github.io/python-awips/api/ModelSounding.html>`__

**matplotlib:**

-  `matplotlib.pyplot <https://matplotlib.org/3.3.3/api/_as_gen/matplotlib.pyplot.html>`__

**MetPy**

-  `metpy.wind_speed <https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.wind_speed.html>`__
-  `metpy.wind_direction <https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.wind_direction.html>`__
-  `metpy.vapor_pressure <https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.vapor_pressure.html>`__
-  `metpy.dewpoint <https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.dewpoint.html>`__
-  `metpy.skewt <https://unidata.github.io/MetPy/latest/api/generated/metpy.plots.SkewT.html>`__
-  `metpy.hodograph <https://unidata.github.io/MetPy/latest/api/generated/metpy.plots.Hodograph.html>`__

`Top <https://unidata.github.io/python-awips/examples/generated/Forecast_Model_Vertical_Sounding.html>`__

--------------