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 cartopy.crs 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 dynamicserialize.dstypes.com.raytheon.uf.common.time import TimeRange
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),
subplot_kw=dict(projection=projection))
ax.set_extent(bbox)
gl = ax.gridlines(draw_labels=True)
gl.top_labels = gl.right_labels = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
return fig, ax
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']
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']
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*.
DataAccessLayer.changeEDEXHost("edex-cloud.unidata.ucar.edu")
request = DataAccessLayer.newDataRequest()
request.setDatatype("warning")
params = ["phensig", "sig"]
request.setParameters(*(params))
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).
WATCH_SIG = 'A'
WARN_SIG = 'W'
ADVIS_SIG = 'Y'
STATEM_SIG = 'S'
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)
Now that we have our request
and TimeRange timerange
objects ready, we're ready to get the data array from EDEX.
# Get response
response = DataAccessLayer.getGeometryData(request, timerange)
print("Using " + str(len(response)) + " records")
Using 1514 records
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))
else:
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:
phensigs.append(phensigString)
# 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:
warning_shapes.append(shape_feature)
elif sig is WATCH_SIG:
watch_shapes.append(shape_feature)
elif sig is ADVIS_SIG:
advisory_shapes.append(shape_feature)
elif sig is STATEM_SIG:
statement_shapes.append(shape_feature)
# 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
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(
category='cultural',
name='admin_1_states_provinces_lines',
scale='50m',
facecolor='none')
political_boundaries = cfeature.NaturalEarthFeature(category='cultural',
name='admin_0_boundary_lines_land',
scale='50m', facecolor='none')
Here is where we finally get to 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 its 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
DRAW_ADVISORY = True
DRAW_WATCH = True
DRAW_WARNING = True
DRAW_STATEMENT = True
# 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)
handles.append(entry)
# Create the plot
bbox=[-127,-64,24,49]
fig, ax = make_map(bbox=bbox)
# Add the title
# Construct the title based on which record types are being displayed
title_string = ""
if DRAW_WATCH:
title_string += "Watches, "
if DRAW_WARNING:
title_string += "Warnings, "
if DRAW_ADVISORY:
title_string += "Advisories, "
if DRAW_STATEMENT:
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(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(states_provinces, edgecolor='black')
ax.add_feature(political_boundaries, edgecolor='black')
# Draw WWAs in order: Advisory -> Watch > Warning > Statement
if DRAW_ADVISORY:
for shape in advisory_shapes:
ax.add_feature(shape)
if DRAW_WATCH:
for shape in watch_shapes:
ax.add_feature(shape)
if DRAW_WARNING:
for shape in warning_shapes:
ax.add_feature(shape)
if DRAW_STATEMENT:
for shape in statement_shapes:
ax.add_feature(shape)
# 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
plt.show()
python-awips
datetime
cartopy
matplotlib