DIFAX Replication
This example replicates the traditional DIFAX images for upper-level observations.
By: Kevin Goebbert
Observation data comes from Iowa State Archive, accessed through the Siphon package. Contour data comes from the GFS 0.5 degree analysis. Classic upper-level data of Geopotential Height and Temperature are plotted.
import urllib.request
from datetime import datetime, timedelta
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
from matplotlib.ticker import FixedLocator
import metpy.calc as mpcalc
import numpy as np
import xarray as xr
from metpy.plots import StationPlot
from metpy.units import units
from siphon.simplewebservice.iastate import IAStateUpperAir
Plotting High/Low Symbols
A helper function to plot a text symbol (e.g., H, L) for relative maximum/minimum for a given field (e.g., geopotential height).
def plot_maxmin_points(lon, lat, data, extrema, nsize, symbol, color='k',
plotValue=True, transform=None):
"""
This function will find and plot relative maximum and minimum for a 2D grid. The function
can be used to plot an H for maximum values (e.g., High pressure) and an L for minimum
values (e.g., low pressue). It is best to used filetered data to obtain a synoptic scale
max/min value. The symbol text can be set to a string value and optionally the color of the
symbol and any plotted value can be set with the parameter color.
Parameters
----------
lon : 2D array
Plotting longitude values
lat : 2D array
Plotting latitude values
data : 2D array
Data that you wish to plot the max/min symbol placement
extrema : str
Either a value of max for Maximum Values or min for Minimum Values
nsize : int
Size of the grid box to filter the max and min values to plot a reasonable number
symbol : str
Text to be placed at location of max/min value
color : str
Name of matplotlib colorname to plot the symbol (and numerical value, if plotted)
plot_value : Boolean (True/False)
Whether to plot the numeric value of max/min point
Return
------
The max/min symbol will be plotted on the current axes within the bounding frame
(e.g., clip_on=True)
"""
from scipy.ndimage import maximum_filter, minimum_filter
if (extrema == 'max'):
data_ext = maximum_filter(data, nsize, mode='nearest')
elif (extrema == 'min'):
data_ext = minimum_filter(data, nsize, mode='nearest')
else:
raise ValueError('Value for hilo must be either max or min')
if lon.ndim == 1:
lon, lat = np.meshgrid(lon, lat)
mxx, mxy = np.where(data_ext == data)
for i in range(len(mxy)):
ax.text(lon[mxx[i], mxy[i]], lat[mxx[i], mxy[i]], symbol, color=color, size=36,
clip_on=True, horizontalalignment='center', verticalalignment='center',
transform=transform)
ax.text(lon[mxx[i], mxy[i]], lat[mxx[i], mxy[i]],
'\n' + str(np.int64(data[mxx[i], mxy[i]])),
color=color, size=12, clip_on=True, fontweight='bold',
horizontalalignment='center', verticalalignment='top', transform=transform)
ax.plot(lon[mxx[i], mxy[i]], lat[mxx[i], mxy[i]], marker='o', markeredgecolor='black',
markerfacecolor='white', transform=transform)
ax.plot(lon[mxx[i], mxy[i]], lat[mxx[i], mxy[i]],
marker='x', color='black', transform=transform)
Station Information
A helper function for obtaining radiosonde station information (e.g., latitude/longitude) requried to plot data obtained from each station. Original code by github user sgdecker.
def station_info(stid):
r"""Provide information about weather stations.
Parameters
----------
stid: str or iterable object containing strs
The ICAO or IATA code(s) for which station information is requested.
with_units: bool
Whether to include units for values that have them. Default True.
Returns
-------
info: dict
Information about the station(s) within a dictionary with these keys:
'state': Two-character ID of the state/province where the station is located,
if applicable
'name': The name of the station
'lat': The latitude of the station [deg]
'lon': The longitude of the station [deg]
'elevation': The elevation of the station [m]
'country': Two-character ID of the country where the station is located
Modified code from Steven Decker, Rutgers University
"""
# Provide a helper function for later usage
def str2latlon(s):
deg = float(s[:3])
mn = float(s[-3:-1])
if s[-1] == 'S' or s[-1] == 'W':
deg = -deg
mn = -mn
return deg + mn / 60.
# Various constants describing the underlying data
url = 'https://www.aviationweather.gov/docs/metar/stations.txt'
# file = 'stations.txt'
state_bnds = slice(0, 2)
name_bnds = slice(3, 19)
icao_bnds = slice(20, 24)
iata_bnds = slice(26, 29)
lat_bnds = slice(39, 45)
lon_bnds = slice(47, 54)
z_bnds = slice(55, 59)
cntry_bnds = slice(81, 83)
# Generalize to any number of IDs
if isinstance(stid, str):
stid = [stid]
# Get the station dataset
infile = urllib.request.urlopen(url)
data = infile.readlines()
state = []
name = []
lat = []
lon = []
z = []
cntry = []
for s in stid:
s = s.upper()
for line_bytes in data:
line = line_bytes.decode('UTF-8')
icao = line[icao_bnds]
iata = line[iata_bnds]
if len(s) == 3 and s in iata or len(s) == 4 and s in icao:
state.append(line[state_bnds].strip())
name.append(line[name_bnds].strip())
lat.append(str2latlon(line[lat_bnds]))
lon.append(str2latlon(line[lon_bnds]))
z.append(float(line[z_bnds]))
cntry.append(line[cntry_bnds])
break
else:
state.append('NA')
name.append('NA')
lat.append(np.nan)
lon.append(np.nan)
z.append(np.nan)
cntry.append('NA')
infile.close()
return {'state': np.array(state), 'name': np.array(name), 'lat': np.array(lat),
'lon': np.array(lon), 'elevation': np.array(z), 'country': np.array(cntry),
'units': {'lat': 'deg', 'lon': 'deg', 'z': 'm'}}
Observation Data
Set a date and time for upper-air observations (should only be 00 or 12 UTC for the hour).
Request all data from Iowa State using the Siphon package. The result is a pandas DataFrame containing all of the sounding data from all available stations.
# Set date for desired UPA data
today = datetime.utcnow()
# Go back one day to ensure data availability
date = datetime(today.year, today.month, today.day, 0) - timedelta(days=1)
# Request data using Siphon request for data from Iowa State Archive
data = IAStateUpperAir.request_all_data(date)
Subset Observational Data
From the request above will give all levels from all radisonde sites available through the service. For plotting a pressure surface map there is only need to have the data from that level. Below the data is subset and a few parameters set based on the level chosen. Additionally, the station information is obtained and latitude and longitude data is added to the DataFrame.
level = 500
if (level == 925) | (level == 850) | (level == 700):
cint = 30
def hght_format(v): return format(v, '.0f')[1:]
elif level == 500:
cint = 60
def hght_format(v): return format(v, '.0f')[:3]
elif level == 300:
cint = 120
def hght_format(v): return format(v, '.0f')[:3]
elif level < 300:
cint = 120
def hght_format(v): return format(v, '.0f')[1:4]
# Create subset of all data for a given level
data_subset = data.pressure == level
df = data[data_subset]
# Get station lat/lon from look-up file; add to Dataframe
stn_info = station_info(list(df.station.values))
df.insert(10, 'latitude', stn_info['lat'])
df.insert(11, 'longitude', stn_info['lon'])
---------------------------------------------------------------------------
HTTPError Traceback (most recent call last)
Cell In[5], line 21
18 df = data[data_subset]
20 # Get station lat/lon from look-up file; add to Dataframe
---> 21 stn_info = station_info(list(df.station.values))
22 df.insert(10, 'latitude', stn_info['lat'])
23 df.insert(11, 'longitude', stn_info['lon'])
Cell In[3], line 52, in station_info(stid)
49 stid = [stid]
51 # Get the station dataset
---> 52 infile = urllib.request.urlopen(url)
53 data = infile.readlines()
55 state = []
File ~/miniconda3/envs/cookbook-dev/lib/python3.10/urllib/request.py:216, in urlopen(url, data, timeout, cafile, capath, cadefault, context)
214 else:
215 opener = _opener
--> 216 return opener.open(url, data, timeout)
File ~/miniconda3/envs/cookbook-dev/lib/python3.10/urllib/request.py:525, in OpenerDirector.open(self, fullurl, data, timeout)
523 for processor in self.process_response.get(protocol, []):
524 meth = getattr(processor, meth_name)
--> 525 response = meth(req, response)
527 return response
File ~/miniconda3/envs/cookbook-dev/lib/python3.10/urllib/request.py:634, in HTTPErrorProcessor.http_response(self, request, response)
631 # According to RFC 2616, "2xx" code indicates that the client's
632 # request was successfully received, understood, and accepted.
633 if not (200 <= code < 300):
--> 634 response = self.parent.error(
635 'http', request, response, code, msg, hdrs)
637 return response
File ~/miniconda3/envs/cookbook-dev/lib/python3.10/urllib/request.py:563, in OpenerDirector.error(self, proto, *args)
561 if http_err:
562 args = (dict, 'default', 'http_error_default') + orig_args
--> 563 return self._call_chain(*args)
File ~/miniconda3/envs/cookbook-dev/lib/python3.10/urllib/request.py:496, in OpenerDirector._call_chain(self, chain, kind, meth_name, *args)
494 for handler in handlers:
495 func = getattr(handler, meth_name)
--> 496 result = func(*args)
497 if result is not None:
498 return result
File ~/miniconda3/envs/cookbook-dev/lib/python3.10/urllib/request.py:643, in HTTPDefaultErrorHandler.http_error_default(self, req, fp, code, msg, hdrs)
642 def http_error_default(self, req, fp, code, msg, hdrs):
--> 643 raise HTTPError(req.full_url, code, msg, hdrs, fp)
HTTPError: HTTP Error 308: Permanent Redirect
Gridded Data
Obtain GFS gridded output for contour plotting. Specifically, geopotential height and temperature data for the given level and subset for over North America. Data are smoothed for aesthetic reasons.
# Get GFS data and subset to North America for Geopotential Height and Temperature
ds = xr.open_dataset('https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/Global_0p5deg_ana/'
'GFS_Global_0p5deg_ana_{0:%Y%m%d}_{0:%H}00.grib2'.format(
date)).metpy.parse_cf()
# Geopotential height and smooth
hght = ds.Geopotential_height_isobaric.metpy.sel(
vertical=level*units.hPa, time=date, lat=slice(70, 15), lon=slice(360-145, 360-50))
smooth_hght = mpcalc.smooth_n_point(hght, 9, 10)
# Temperature, smooth, and convert to Celsius
tmpk = ds.Temperature_isobaric.metpy.sel(
vertical=level*units.hPa, time=date, lat=slice(70, 15), lon=slice(360-145, 360-50))
smooth_tmpc = (mpcalc.smooth_n_point(tmpk, 9, 10)).metpy.convert_units('degC')
Create DIFAX Replication
Plot the observational data and contours on a Lambert Conformal map and add features that resemble the historic DIFAX maps.
# Set up map coordinate reference system
mapcrs = ccrs.LambertConformal(
central_latitude=45, central_longitude=-100, standard_parallels=(30, 60))
# Set up station locations for plotting observations
point_locs = mapcrs.transform_points(
ccrs.PlateCarree(), df['longitude'].values, df['latitude'].values)
# Start figure and set graphics extent
fig = plt.figure(1, figsize=(17, 15))
ax = plt.subplot(111, projection=mapcrs)
ax.set_extent([-125, -70, 20, 55])
# Add map features for geographic reference
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), edgecolor='grey')
ax.add_feature(cfeature.LAND.with_scale('50m'), facecolor='white')
ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='grey')
# Plot plus signs every degree lat/lon
plus_lat = []
plus_lon = []
other_lat = []
other_lon = []
for x in hght.lon.values[::2]:
for y in hght.lat.values[::2]:
if (x % 5 == 0) | (y % 5 == 0):
plus_lon.append(x)
plus_lat.append(y)
else:
other_lon.append(x)
other_lat.append(y)
ax.scatter(other_lon, other_lat, s=2, marker='o',
transform=ccrs.PlateCarree(), color='lightgrey', zorder=-1)
ax.scatter(plus_lon, plus_lat, s=30, marker='+', transform=ccrs.PlateCarree(),
color='lightgrey')
# Add gridlines for every 5 degree lat/lon
ax.gridlines(linestyle='solid', ylocs=range(15, 71, 5), xlocs=range(-150, -49, 5))
# Start the station plot by specifying the axes to draw on, as well as the
# lon/lat of the stations (with transform). We also the fontsize to 10 pt.
stationplot = StationPlot(ax, df['longitude'].values, df['latitude'].values, clip_on=True,
transform=ccrs.PlateCarree(), fontsize=10)
# Plot the temperature and dew point to the upper and lower left, respectively, of
# the center point.
stationplot.plot_parameter('NW', df['temperature'], color='black')
stationplot.plot_parameter('SW', df['dewpoint'], color='black')
# A more complex example uses a custom formatter to control how the geopotential height
# values are plotted. This is set in an earlier if-statement to work appropriate for
# different levels.
stationplot.plot_parameter('NE', df['height'], formatter=hght_format)
# Add wind barbs
stationplot.plot_barb(df['u_wind'], df['v_wind'], length=7, pivot='tip')
# Plot Solid Contours of Geopotential Height
cs = ax.contour(hght.lon, hght.lat, smooth_hght,
range(0, 20000, cint), colors='black', transform=ccrs.PlateCarree())
clabels = plt.clabel(cs, fmt='%d', colors='white', inline_spacing=5, use_clabeltext=True)
# Contour labels with black boxes and white text
for t in cs.labelTexts:
t.set_bbox({'facecolor': 'black', 'pad': 4})
t.set_fontweight('heavy')
# Plot Dashed Contours of Temperature
cs2 = ax.contour(hght.lon, hght.lat, smooth_tmpc, range(-60, 51, 5),
colors='black', transform=ccrs.PlateCarree())
clabels = plt.clabel(cs2, fmt='%d', colors='black', inline_spacing=5, use_clabeltext=True)
# Set longer dashes than default
for c in cs2.collections:
c.set_dashes([(0, (5.0, 3.0))])
# Contour labels with black boxes and white text
for t in cs.labelTexts:
t.set_bbox({'facecolor': 'black', 'pad': 4})
t.set_fontweight('heavy')
# Plot filled circles for Radiosonde Obs
ax.scatter(df['longitude'].values, df['latitude'].values, s=10,
marker='o', color='black', transform=ccrs.PlateCarree())
# Use definition to plot H/L symbols
plot_maxmin_points(hght.lon, hght.lat, smooth_hght.values, 'max', 50,
symbol='H', color='black', transform=ccrs.PlateCarree())
plot_maxmin_points(hght.lon, hght.lat, smooth_hght.values, 'min', 25,
symbol='L', color='black', transform=ccrs.PlateCarree())
# Add titles
plt.title(f'Upper-air Observations at {level}-hPa Analysis Heights/Temperature',
loc='left')
plt.title(f'Valid: {date}', loc='right');
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/pandas/core/indexes/base.py:3805, in Index.get_loc(self, key)
3804 try:
-> 3805 return self._engine.get_loc(casted_key)
3806 except KeyError as err:
File index.pyx:167, in pandas._libs.index.IndexEngine.get_loc()
File index.pyx:196, in pandas._libs.index.IndexEngine.get_loc()
File pandas/_libs/hashtable_class_helper.pxi:7081, in pandas._libs.hashtable.PyObjectHashTable.get_item()
File pandas/_libs/hashtable_class_helper.pxi:7089, in pandas._libs.hashtable.PyObjectHashTable.get_item()
KeyError: 'longitude'
The above exception was the direct cause of the following exception:
KeyError Traceback (most recent call last)
Cell In[7], line 7
2 mapcrs = ccrs.LambertConformal(
3 central_latitude=45, central_longitude=-100, standard_parallels=(30, 60))
5 # Set up station locations for plotting observations
6 point_locs = mapcrs.transform_points(
----> 7 ccrs.PlateCarree(), df['longitude'].values, df['latitude'].values)
9 # Start figure and set graphics extent
10 fig = plt.figure(1, figsize=(17, 15))
File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/pandas/core/frame.py:4102, in DataFrame.__getitem__(self, key)
4100 if self.columns.nlevels > 1:
4101 return self._getitem_multilevel(key)
-> 4102 indexer = self.columns.get_loc(key)
4103 if is_integer(indexer):
4104 indexer = [indexer]
File ~/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/pandas/core/indexes/base.py:3812, in Index.get_loc(self, key)
3807 if isinstance(casted_key, slice) or (
3808 isinstance(casted_key, abc.Iterable)
3809 and any(isinstance(x, slice) for x in casted_key)
3810 ):
3811 raise InvalidIndexError(key)
-> 3812 raise KeyError(key) from err
3813 except TypeError:
3814 # If we have a listlike key, _check_indexing_error will raise
3815 # InvalidIndexError. Otherwise we fall through and re-raise
3816 # the TypeError.
3817 self._check_indexing_error(key)
KeyError: 'longitude'
