Regional Surface Obs Plot
Notebook Python-AWIPS Tutorial Notebook
Objectives
Use python-awips to connect to an edex server
Create a plot for a regional area of the United States (Florida)
Define and filter data request for METAR and Synoptic surface obs
Use the maps database to request and draw state boundaries (no use of Cartopy.Feature in this example)
Stylize and plot surface data using Metpy
Table of Contents
Imports
The imports below are used throughout the notebook. Note the first two imports are coming directly from python-awips and allow us to connect to an EDEX server, and define a timrange used for filtering the data. The subsequent imports are for data manipulation and visualization.
from awips.dataaccess import DataAccessLayer
from dynamicserialize.dstypes.com.raytheon.uf.common.time import TimeRange
from datetime import datetime, timedelta
import numpy as np
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.feature import ShapelyFeature
from shapely.geometry import Polygon
import matplotlib.pyplot as plt
from metpy.units import units
from metpy.calc import wind_components
from metpy.plots import simple_layout, StationPlot, StationPlotLayout, sky_cover
import warnings
Function: get_cloud_cover()
Returns the cloud coverage values as integer codes (0 through 8).
def get_cloud_cover(code):
if 'OVC' in code:
return 8
elif 'BKN' in code:
return 6
elif 'SCT' in code:
return 4
elif 'FEW' in code:
return 2
else:
return 0
Function: make_map()
In order to plot more than one image, it’s easiest to define common logic in a function. Here, a new function called make_map is defined. This function uses the matplotlib.pyplot package (plt) to create a figure and axis. The geographic extent is set and lat/lon gridlines are added for context.
def make_map(bbox, proj=ccrs.PlateCarree()):
fig, ax = plt.subplots(figsize=(16,12),subplot_kw=dict(projection=proj))
ax.set_extent(bbox)
gl = ax.gridlines(draw_labels=True, color='#e7e7e7')
gl.top_labels = gl.right_labels = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
return fig, ax
Function: extract_plotting_data()
Grab the simple variables out of the response data we have (attaching correct units), and put them into a dictionary that we will hand the plotting function later:
Get wind components from speed and direction
Convert cloud coverage values to integer codes [0 - 8]
Assign temperature, dewpoint, and sea level pressure the the correct units
Account for missing values (by using
nan)
def extract_plotting_data(arr, datatype):
"""
Extract all necessary data for plotting for either
datatype: 'obs' or 'sfcobs'
"""
data = dict()
data['latitude'] = np.array(arr['latitude'])
data['longitude'] = np.array(arr['longitude'])
tmp = np.array(arr['temperature'], dtype=float)
dpt = np.array(arr['dewpoint'], dtype=float)
direction = np.array(arr['windDir'])
# Suppress nan masking warnings
warnings.filterwarnings("ignore",category =RuntimeWarning)
# Account for missing values
tmp[tmp == -9999.0] = 'nan'
dpt[dpt == -9999.] = 'nan'
direction[direction == -9999.0] = 'nan'
data['air_pressure_at_sea_level'] = np.array(arr['seaLevelPress'])* units('mbar')
u, v = wind_components(np.array(arr['windSpeed']) * units('knots'),
direction * units.degree)
data['eastward_wind'], data['northward_wind'] = u, v
data['present_weather'] = arr['presWeather']
# metars uses 'stationName' for its identifier and temps are in deg C
# metars also has sky coverage
if datatype == "obs":
data['stid'] = np.array(arr['stationName'])
data['air_temperature'] = tmp * units.degC
data['dew_point_temperature'] = dpt * units.degC
data['cloud_coverage'] = [int(get_cloud_cover(x)) for x in arr['skyCover']]
# synoptic obs uses 'stationId', and temps are in Kelvin
elif datatype == "sfcobs":
data['stid'] = np.array(arr['stationId'])
data['air_temperature'] = tmp * units.kelvin
data['dew_point_temperature'] = dpt * units.kelvin
return data
Function: plot_data()
This function makes use of Metpy.StationPlotLayout and Metpy.StationPlot
to add all surface observation data to our plot. The logic is very
similar for both METAR and Synoptic data, so a datatype argument is
used to distinguish between which data is being drawn, and then draws
the appropriate features.
This function plots: - Wind barbs - Air temperature - Dew point temperature - Precipitation - Cloud coverage (for METARS)
def plot_data(data, title, axes, datatype):
custom_layout = StationPlotLayout()
custom_layout.add_barb('eastward_wind', 'northward_wind', units='knots')
custom_layout.add_value('NW', 'air_temperature', fmt='.0f', units='degF', color='darkred')
custom_layout.add_value('SW', 'dew_point_temperature', fmt='.0f', units='degF', color='darkgreen')
custom_layout.add_value('E', 'precipitation', fmt='0.1f', units='inch', color='blue')
# metars has sky coverage
if datatype == 'obs':
custom_layout.add_symbol('C', 'cloud_coverage', sky_cover)
axes.set_title(title)
stationplot = StationPlot(axes, data['longitude'], data['latitude'], clip_on=True,
transform=ccrs.PlateCarree(), fontsize=10)
custom_layout.plot(stationplot, data)
Initial Setup
Connect to an EDEX server and define several new data request objects.
In this example we’re using multiple different datatypes from EDEX, so we’ll create a request object for each of the following: - The states outlines (datatype maps) - The METAR data (datatype obs) - The Synoptic data (datatype sfc)
Some of the request use filters, so we’ll also create several filters than can be used for the various data requests as well.
Initial EDEX Connection
First we establish a connection to Unidata’s public EDEX server.
# EDEX connection
edexServer = "edex-cloud.unidata.ucar.edu"
DataAccessLayer.changeEDEXHost(edexServer)
Maps Request and Response
The maps data request will give us data to draw our state outlines of interest (Florida and its neighboring states). We will retrieve the data response object here so we can create a geographic filter for the METAR and Synoptic data requests.
# Define the maps request
maps_request = DataAccessLayer.newDataRequest('maps')
# filter for multiple states
maps_request.addIdentifier('table', 'mapdata.states')
maps_request.addIdentifier('geomField', 'the_geom')
maps_request.addIdentifier('inLocation', 'true')
maps_request.addIdentifier('locationField', 'state')
maps_request.setParameters('state','name','lat','lon')
maps_request.setLocationNames('FL','GA','MS','AL','SC','LA')
maps_response = DataAccessLayer.getGeometryData(maps_request)
print("Found " + str(len(maps_response)) + " MultiPolygons")
Found 6 MultiPolygons
Define Geographic Filter
The previous EDEX request limited the data by using a parameter for the maps database called state. We can take the results from that filter and get a geographic envelope based on the Florida polygon that was returned from the previous cell.
Warning: Without such a filter you may be requesting many tens of thousands of records.
# Append each geometry to a numpy array
states = np.array([])
for ob in maps_response:
print(ob.getString('name'), ob.getString('state'), ob.getNumber('lat'), ob.getNumber('lon'))
states = np.append(states,ob.getGeometry())
# if this is Florida grab geographic info
if ob.getString('name') == "Florida":
bounds = ob.getGeometry().bounds
fl_lat = ob.getNumber('lat')
fl_lon = ob.getNumber('lon')
if bounds is None:
print("Error, no record found for Florida!")
else:
# buffer our bounds by +/i degrees lat/lon
bbox=[bounds[0]-3,bounds[2]+3,bounds[1]-1.5,bounds[3]+1.5]
# Create envelope geometry
envelope = Polygon([(bbox[0],bbox[2]),(bbox[0],bbox[3]),
(bbox[1], bbox[3]),(bbox[1],bbox[2]),
(bbox[0],bbox[2])])
print(envelope)
Florida FL 28.67402 -82.50934
Georgia GA 32.65155 -83.44848
Louisiana LA 31.0891 -92.02905
Alabama AL 32.79354 -86.82676
Mississippi MS 32.75201 -89.66553
South Carolina SC 33.93574 -80.89899
POLYGON ((-90.63429260299995 23.02105161600002, -90.63429260299995 32.50101280200016, -77.03199876199994 32.50101280200016, -77.03199876199994 23.02105161600002, -90.63429260299995 23.02105161600002))
Define Time Filter
Both the METAR and Synoptic datasets should be filtered by time to avoid requesting an unreasonable amount of data. By defining one filter now, we can use it in both of their data requests to EDEX.
Note: Here we will use the most recent hour as our default filter. Try adjusting the timerange and see the difference in the final plots.
# Filter for the last hour
lastHourDateTime = datetime.utcnow() - timedelta(minutes = 60)
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)
print(timerange)
(Nov 11 22 19:00:54 , Nov 11 22 20:00:54 )
Define Common Parameters for Data Requests
METAR obs and Synoptic obs share several of the same parameters. By defining them here, they can be reused for both of the requests and this makes our code more efficient.
shared_params = ["timeObs", "longitude", "latitude", "temperature",
"dewpoint", "windDir", "windSpeed", "seaLevelPress",
"presWeather", "skyLayerBase"]
print(shared_params)
['timeObs', 'longitude', 'latitude', 'temperature', 'dewpoint', 'windDir', 'windSpeed', 'seaLevelPress', 'presWeather', 'skyLayerBase']
Define METAR Request
To get METAR data we must use the obs datatype. To help limit the amount of data returned, we will narrow the request by using a geographic envelope, setting the request parameters, and using timerange as a time filter.
# New metar request
metar_request = DataAccessLayer.newDataRequest("obs", envelope=envelope)
# metar specifc parameters
metar_params = ["stationName", "skyCover"]
# combine all parameters
all_metar_params = shared_params + metar_params
# set the parameters on the metar request
metar_request.setParameters(*(all_metar_params))
print(metar_request)
DefaultDataRequest(datatype=obs, identifiers={}, parameters=['timeObs', 'longitude', 'latitude', 'temperature', 'dewpoint', 'windDir', 'windSpeed', 'seaLevelPress', 'presWeather', 'skyLayerBase', 'stationName', 'skyCover'], levels=[], locationNames=[], envelope=<dynamicserialize.dstypes.com.vividsolutions.jts.geom.Envelope.Envelope object at 0x13abe40a0>)
Define Synoptic Request
Similar to the request above, we will limit the amount of data returned by using a geographic envelope, setting the request parameters, and using timerange as a time filter.
However, in order to access synoptic observations we will use the sfcobs datatype.
# New sfcobs/SYNOP request
syn_request = DataAccessLayer.newDataRequest("sfcobs", envelope=envelope)
# (sfcobs) uses stationId, while (obs) uses stationName
syn_params = ["stationId"]
# combine all parameters
all_syn_params = shared_params + syn_params
# set the parameters on the synoptic request
syn_request.setParameters(*(all_syn_params))
print(syn_request)
DefaultDataRequest(datatype=sfcobs, identifiers={}, parameters=['timeObs', 'longitude', 'latitude', 'temperature', 'dewpoint', 'windDir', 'windSpeed', 'seaLevelPress', 'presWeather', 'skyLayerBase', 'stationId'], levels=[], locationNames=[], envelope=<dynamicserialize.dstypes.com.vividsolutions.jts.geom.Envelope.Envelope object at 0x105048bb0>)
Get the Data!
We have already obtained our maps data, but we still have to collect our observation data.
Get the EDEX Responses
# METARs data
metar_response = DataAccessLayer.getGeometryData(metar_request,timerange)
# function getMetarObs was added in python-awips 18.1.4
metars = DataAccessLayer.getMetarObs(metar_response)
print("Found " + str(len(metar_response)) + " METAR records")
print("\tUsing " + str(len(metars['temperature'])) + " temperature records")
# Synoptic data
syn_response = DataAccessLayer.getGeometryData(syn_request,timerange)
# function getSynopticObs was added in python-awips 18.1.4
synoptic = DataAccessLayer.getSynopticObs(syn_response)
print("Found " + str(len(syn_response)) + " Synoptic records")
print("\tUsing " + str(len(synoptic['temperature'])) + " temperature records")
Found 4116 METAR records
Using 179 temperature records
Found 259 Synoptic records
Using 63 temperature records
Extract Plotting Data
# Pull all necessary plotting information for the metar data
metars_data = extract_plotting_data(metars, 'obs')
print(str(len(metars_data['stid'])) + " METARs stations")
# Pull all necessary plotting information for the synoptic data
synoptic_data = extract_plotting_data(synoptic, 'sfcobs')
print(str(len(synoptic_data['stid'])) + " Synoptic stations")
179 METARs stations
63 Synoptic stations
Plot the Data
Draw the Region
Here we will draw our region by using the states polygons we retreived from EDEX earlier in this example. To create this plot we use the make_map() function which also adds lines of latitude and longitude for additional context.
# Create the figure and axes used for the plot
fig, ax = make_map(bbox=bbox)
# Create a feature based off our states polygons
shape_feature = ShapelyFeature(states,ccrs.PlateCarree(),
facecolor='none', linestyle="-",edgecolor='#000000',linewidth=2)
ax.add_feature(shape_feature)
<cartopy.mpl.feature_artist.FeatureArtist at 0x13b2ae5e0>
Plot METAR Data
On the same axes (ax) and figure (fig) plot the METAR data.
# Create a title for the plot
title = str(metar_response[-1].getDataTime()) + " | METAR Surface Obs | " + edexServer
# Plot the station information for METARs data
plot_data(metars_data, title, ax, 'obs')
# Display the figure
fig
Plot Synoptic Data
On a new axes and figure (ax_syn, fig_syn) plot the map and synoptic data.
# Create a new figure and axes for the synoptic data
fig_syn, ax_syn = make_map(bbox=bbox)
# Create the states feature from the polygons
shape_feature = ShapelyFeature(states,ccrs.PlateCarree(),
facecolor='none', linestyle="-",edgecolor='#000000',linewidth=2)
ax_syn.add_feature(shape_feature)
# Create a title for the figure
title = str(syn_response[-1].getDataTime()) + " | SYNOP Surface Obs | " + edexServer
# Draw the synoptic data
plot_data(synoptic_data, title, ax_syn, 'sfcobs')
Plot both METAR and Synoptic Data
Add the synoptic data to our first axes and figure (ax, fig) that already contains our map and METARs data.
# Create a title for both the METAR and Synopotic data
title = str(syn_response[-1].getDataTime()) + " | METAR & Synoptic Surface Obs | " + edexServer
# Draw the synoptic on the first axes that already has the metar data
plot_data(synoptic_data, title, ax, 'sfcobs')
# Display the figure
fig
See Also
Additional Documentation
python-awips:
datetime:
numpy:
cartopy:
matplotlib:
metpy: