Watch Warning and Advisory Plotting

Notebook Python-AWIPS Tutorial Notebook


  • Create a colorized plot with Warnings, Watches, Advisories and Statements (WWAs)

  • Use python-awips to connect to an EDEX server

  • Create and filter the data request specifically for a warning data type

  • Create and use accurate time filter for data requests

  • Define and use functions

  • Define and use dictionaries

  • Colorize shapes based on a dictionary

  • Overlay warnings, watches, and advisories with state and political maps

Table of Contents

1 Imports

The imports below are used throughout the notebook. The python-awips imports allow us to connect to an EDEX server, use the warning lookup dictionary, and define a TimeRange. The additional imports are for data manipulation and visualization.

from datetime import datetime, timedelta

import numpy as np
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
import as ccrs
import cartopy.feature as cfeature
from cartopy.feature import ShapelyFeature, NaturalEarthFeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from shapely.geometry import MultiPolygon, Polygon

from awips.dataaccess import DataAccessLayer
from awips.tables import vtec
from import TimeRange


2 Function: make_map()

In order to plot more than one image, it’s easiest to define common logic in a function. However, for this notebook we only use it in one place. It is a function you will find in most of our example notebooks.

Here, a new function called make_map is defined. This function uses the matplotlib.pyplot package (plt) to create a figure and axis. The lat/lon grids are added.

def make_map(bbox, projection=ccrs.PlateCarree()):
    fig, ax = plt.subplots(figsize=(20,12),
    gl = ax.gridlines(draw_labels=True)
    gl.top_labels = gl.right_labels = False
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    return fig, ax


3 Function: get_color()

Since we’ll be needing to access the color using the vtec lookup table in several places, creating an easily recognizable function is useful.

def get_color(phensig):
    return vtec[phensig]['color']


4 Function get_title()

Similar to the color function just defined, accessing the full name for the phensig will also be necessary, so this function will be helpful.

def get_title(phensig):
    return vtec[phensig]['hdln']


5 Initial Setup

5.1 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 warning, and set the Parameters to phensig and sig.

Note: Remember, to see all available parameters use the DataAccess.getAvailableParameters() method as shown in the Grid Levels and Parameters Notebook.

request = DataAccessLayer.newDataRequest()
params = ["phensig", "sig"]

5.2 Significance (Sig) Constants

The two parameters we’re requesting for our warning objects are phensig and sig where phensig is styled “XX.Y” and sig is “Y”. Phen stands for “Phenomena” and sig stands for “Significance”. A more detailed description of phensigs and how they’re used is provided with this NWS pamphlet.

The constants in this section correlate the sig to what type of message it is (what significance it is).



6 Filter by Time

Here we decide how much data we want to pull from EDEX. By default we’ll request 12 hours, but that value can easily be modified by adjusting the ``timedelta(hours = 12)` <>`__ in line 2. The more data we request, the longer the next section will take to run.

# Get records from the last 12 hours
lastHourDateTime = datetime.utcnow() - timedelta(hours = 12)
start = lastHourDateTime.strftime('%Y-%m-%d %H:%M:%S')
end = datetime.utcnow().strftime('%Y-%m-%d %H:%M:%S')

beginRange = datetime.strptime( start , "%Y-%m-%d %H:%M:%S")
endRange = datetime.strptime( end , "%Y-%m-%d %H:%M:%S")
timerange = TimeRange(beginRange, endRange)


7 Use the Data!

7.1 Get the Data

Now that we have our request and TimeRange timerange objects ready, it’s time to request the data array from EDEX.

Note: Above we set timerange to be 12 hours worth of data. This can return on the order of ~2000 records and can take a little while to run.

# Get response
response = DataAccessLayer.getGeometryData(request, timerange)
print("Using " + str(len(response)) + " records")
Using 1502 records

7.2 Extract Phensigs, Geometries, and Times

In this section we start gathering all the information we’ll need to properly display our data. First we create an array to keep track of unique phensigs. This is useful summary information and will be used to help create the legend which we’ll display along with our plot.

Next, we create arrays for each of the 4 types of significance a statement can have. We will group our records this way, so we can easily toggle which records to display or not.

Then, we create two time variables to keep track of the earliest time from our records and the latest time, and will display that information in the title of our plot.

This section has optional print statements at lines 65 and 85. The first prints out the title, phensig, ref time, and shape for each unique phensig, and the second prints out a sum of how many unique phensigs there are.

We cycle through all the data produced from our response object, access its geometries, and create a new ShapelyFeature with the corresponding color. Then we place this new feature in the appropriate shapes array. During this process we also populate the phensigs array with all unique phensig entries.

Finally, after we’re done looping through all the response data, we create a mapping of phensigs to their corresponding titles. This will be used later to sort the legend alphabetically by titles (which differs from simply sorting by phensig). Ex. Blizzard Warning (BZ.W) would come before Areal Flood Advisory (FA.Y) if we simply sorted by phensig.

# Keep track of unique phensigs, to use in legend
phensigs = []

# Sort the geometries based on their sig
watch_shapes = []
warning_shapes = []
advisory_shapes = []
statement_shapes = []

# Keep track of the earliest and latest reftime for title
# start with the first time from the first object in the response
time_str = str(response[0].getDataTime().getRefTime())
# truncate the decimal seconds for datetime parsing
time_str = time_str[:-4]
# convert to datetime object for easy comparison
first_time = datetime.strptime(time_str, '%Y-%m-%d %H:%M:%S')
last_time = datetime.strptime(time_str, '%Y-%m-%d %H:%M:%S')

for ob in response:

    # get the geometry for the object
    poly = ob.getGeometry()
    # get the reftime for the object
    ref  = ob.getDataTime().getRefTime()

    # do not plot if phensig is blank (SPS)
    if ob.getString('phensig'):

        # look at the reftime
        # convert reftime to a string and parse the decimal seconds
        ref_str = str(ref)
        ref_str = ref_str[:-4]
        # convert reftime to a datetime object for comparison
        ref_time = datetime.strptime(ref_str, '%Y-%m-%d %H:%M:%S')
        # compare current time with first and last times and set if appropriate
        if ref_time is not None:
            if ref_time < first_time:
                first_time = ref_time
            elif ref_time > last_time:
                last_time = ref_time

        # get the phensig and sig values from object
        phensigString = ob.getString('phensig')
        sig = ob.getString('sig')

        # set the geometries based on whether it's a MultiPolygon or Polygon
        if poly.geom_type == 'MultiPolygon':
            geometries = np.array([])
            geometries = np.append(geometries,MultiPolygon(poly))
            geometries = np.array([])
            geometries = np.append(geometries,Polygon(poly))

        for geom in geometries:
            bounds = Polygon(geom)
            intersection = bounds.intersection
            geoms = (intersection(geom) for geom in geometries if bounds.intersects(geom))

        # Store the unique phensigs
        if not phensigString in phensigs:
            # Optional printout of unique Phensigs
#             print(get_title(phensigString) + " (" + phensigString + ")

        # get the corresponding color using the dictionary
        color = get_color(phensigString)
        # create a new shape feature for the object
        shape_feature = ShapelyFeature(geoms,ccrs.PlateCarree(),
                        facecolor=color, edgecolor=color)
        # store the shape feature in the correct array
        if sig is WARN_SIG:
        elif sig is WATCH_SIG:
        elif sig is ADVIS_SIG:
        elif sig is STATEM_SIG:

# Optional printout for the number of unique phensigs
print(len(phensigs), " Unique Phensigs")

# Map phensigs to their titles (used for displaying alphabetically by
# title in legend)
phensig_titles = {}
for phensig in phensigs:
    key = get_title(phensig)
    phensig_titles[key] = phensig
14  Unique Phensigs


8 Plot the Data!

8.1 Create State and Political Boundaries

Define the state and political boundaries that we’ll use in our plot to give more of a frame of reference. These objects are standard method calls in the Cartopy Feature package, using the NaturalEarthFeature function.

# Define the state and political boundaries for the plot
states_provinces = cfeature.NaturalEarthFeature(
political_boundaries = cfeature.NaturalEarthFeature(category='cultural',
                               scale='50m', facecolor='none')

8.2 Draw the Plot and Legend for WWAs

Here is where we finally get ot draw something! The very first few lines of this section are constants that we can manually “switch on and off” for what records we want displayed. By default we have all significance types drawn. If we want to “turn off” any of the significance records, simply set it’s corresponding constant to false, and re-run this cell to see how that plot compares.

The next step involves creating the objects that are used to define the legend. We use the phensig_titles dictionary to loop through all the phensigs in alphabetical (by title) order. Then, we compare if the phensig will be displayed or not based on the display constants from the previous lines. If the significance will be drawn then we create a new Patch object of the corresponding color with the corresponding label and add it to our handles array.

After that we define our bounding box and create our new plot with its figure and axes.

Our next step is to create our Title for our plot. We create a title based on the draw variables to accurately describe what is being drawn in our plot. Here is where we use the first and last times defined in a previous cell.

Finally, we create and show our plot. We add the title to the plot, add all the features to the axes, and add the legend as well.

# Set these variables for which records to draw

# Create handles for legend and add items alphabetically by title and
# only display based on the display values above
handles = []
for title in sorted(phensig_titles):
    phensig = phensig_titles[title]
    # check draw booleans
    if ( "."+ADVIS_SIG in phensig and DRAW_ADVISORY or
        "."+WATCH_SIG in phensig and DRAW_WATCH or
        "."+WARN_SIG in phensig and DRAW_WARNING or
        "."+STATEM_SIG in phensig and DRAW_STATEMENT ):
        entry = mpatches.Patch(color=get_color(phensig), label=title)

# Create the plot
fig, ax = make_map(bbox=bbox)

# Add the title
# Construct the title based on which record types are being displayed
title_string = ""
    title_string += "Watches, "
    title_string += "Warnings, "
    title_string += "Advisories, "
    title_string += "Statements, "
# remove the last comma and space
title_string = title_string[:-2]
# add the time range
title_string += " from " + str(first_time)[:-3] + " to " + str(last_time)[:-3] + " UTC"
# set the title on the plot, give it a bigger font size, and increase
# the vertical padding between the title and the figure
plt.title(title_string, fontsize=20, pad=10)

# Draw all features on the plot
ax.add_feature(states_provinces, edgecolor='black')
ax.add_feature(political_boundaries, edgecolor='black')

# Draw WWAs in order: Advisory -> Watch > Warning > Statement
    for shape in advisory_shapes:
    for shape in watch_shapes:
    for shape in warning_shapes:
    for shape in statement_shapes:

# Draw the legend
# use the handles defined earlier for the color associations to
# phensig titles, set the location to the lower center, give it
# 5 columns so it uses all the horizonatal space, place it under
# the actual figure, and give it a larger fontsize
bottom = 0.12 + (len(handles)//5 *.04)
ax.legend(handles=handles, loc='lower center', ncol=5, bbox_to_anchor=(0.5, -bottom), fontsize=16)

# Show the plot


9 See Also

9.2 Additional Documentation