Meteorology in Python
North American Regional Reanalysis 850-mb Temperature Advection Map
# Name: narr_850mb_temp_adv_map.py
#
# Purpose: Extracts numerical weather
# prediction model output from
# North American Regional Reanalysis
# (NARR). Calculates several key
# atmospheric properties from
# output. Plots contours,
# filled-contours, and barbs from output.
# Here I focus on creating a synoptic scale, 850-mb map.
#
# Caveats: Run code with Python 3.6
# or higher. Code requires making
# modification to one line of MetPy's
# cartopy_utils.py code to function
# with later releases of Cartopy
# (see footnote [1] below).
#
# Authors: Source code originally
# developed by University
# Corporation for Atmospheric Research.
# Code modified by Lawrence Burkett
# (2020) for academic use.
# Import the following modules:
import matplotlib
matplotlib.use('Agg')
from matplotlib import pyplot as plt
from datetime import datetime
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
import numpy as np
from metpy.plots import StationPlot
from metpy.units import units
from netCDF4 import Dataset, num2date
from scipy.ndimage import gaussian_filter
from siphon.catalog import TDSCatalog
# Year, month, day, and hour of interest.
year = 2004
month = 5
day = 29
hour = 12
dt = datetime(year, month, day, hour)
# Pressure level of interest in mb.
p0 = 850
# Read NARR data from THREDDS
# server.
base_url = 'https://www.ncei.noaa.gov/thredds/catalog/model-narr-a-files/'
# Generate URL to year, month,
# day, and hour of interest.
cat = TDSCatalog(f'{base_url}{dt:%Y%m}/{dt:%Y%m%d}/catalog.xml')
# Have Siphon find appropriate dataset.
ds = cat.datasets.filter_time_nearest(dt)
# Interface with data through NetCDF
# Subset Service (NCSS).
ncss = ds.subset()
# Create NCSS query with desired
# specifications.
query = ncss.query()
query.lonlat_box(north = 60, south = 18, east = 300, west = 225)
query.all_times()
query.add_lonlat()
query.accept('netcdf')
query.variables('Geopotential_height_isobaric', 'Temperature_isobaric', 'u-component_of_wind_isobaric', 'v-component_of_wind_isobaric')
# Use query to obtain NetCDF data.
data = ncss.get_data(query)
# Extract data and assign units.
# Temperature in K.
h = gaussian_filter(data.variables['Geopotential_height_isobaric'][0], sigma = 1.0) * units.meter
# West-to-east (zonal) wind in m/s.
T = gaussian_filter(data.variables['Temperature_isobaric'][0], sigma = 1.0) * units.K
# Geopotential height in m AMSL.
u = gaussian_filter(data.variables['u-component_of_wind_isobaric'][0], sigma = 1.0) * units('m/s')
# South-to-north (meridional) wind in m/s.
v = gaussian_filter(data.variables['v-component_of_wind_isobaric'][0], sigma = 1.0) * units('m/s')
# Extract coordinate data for plotting (i.e., latitude, longitude, and height).
lat = data.variables['lat'][:]
lon = data.variables['lon'][:]
p = data.variables['isobaric1'][:]
# Extract time for plotting. Check whether
# time2, time1, or time by printing data:
# print(data)
time = data.variables['time2']
vtime = num2date(time[0], units = time.units)
# Calculate dx and dy for calculations.
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
# Model pressure level k in (i, j, k)
# grid index coordinate system.
k = np.where(p == p0)[0][0]
# Geopotential height at k-th
# index in m.
h = h[k]
# Temperature at k-th index
# in K.
T = T[k]
# West-to-east (zonal) wind
# at k-th index in m/s.
u = u[k]
# South-to-north (meridional)
# wind at k-th index in m/s.
v = v[k]
# Temperature advection in C/s.
T_adv = mpcalc.advection(T, [u, v], (dx, dy), dim_order = 'yx').to('degC/s')
# Data projection.
dataproj = ccrs.PlateCarree()
# Plot projection.
plotproj = ccrs.LambertConformal(central_longitude = -100., central_latitude = 40., standard_parallels = [30, 60])
# Define background map.
def create_map_background():
fig = plt.figure(figsize = (14, 12))
ax = plt.subplot(111, projection = plotproj)
ax.set_extent([-125, -73, 25, 50],ccrs.PlateCarree())
ax.coastlines("50m", linewidth = 0.75)
ax.add_feature(cfeature.STATES, linewidth = 0.5)
return fig, ax
# Create figure.
fig, ax = create_map_background()
# First layer are temperatures
# in deg. C contours.
ly1 = ax.contour(lon, lat, T.to('degC'), range(-50, 50, 2), colors = 'grey', linestyles = 'dotted', transform = dataproj)
plt.clabel(ly1, fontsize = 10, inline = 1, inline_spacing = 10, fmt = '%i', rightside_up = True, use_clabeltext = True)
# Second layer are geopotential
# heights in m AMSL contours.
geo_range = np.arange(0, 4000, 30)
ly2 = ax.contour(lon, lat, h, geo_range, colors = 'k', linewidths = 1.0, linestyles = 'solid', transform = dataproj)
plt.clabel(ly2, fontsize = 10, inline = 1, inline_spacing=10, fmt = '%i', rightside_up = True, use_clabeltext = True)
# Third layer are colorfill temperature
# advection in deg. C/h.
adv_range = [-3, -2.2, -2, -1.5, -1, -0.5, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0]
cf = ax.contourf(lon, lat, T_adv*3600, adv_range, cmap = 'bwr', extend = 'both', transform = dataproj)
plt.colorbar(cf, orientation = 'horizontal', pad = 0, aspect = 50, extendrect = True, ticks = adv_range)
# Fourth layer are wind barbs in kt.
ax.barbs(lon, lat, u.to('kts').m, v.to('kts').m, regrid_shape = 15, transform = dataproj)
# Titles
plt.title('NARR 850-mb geo heights (m), temp (\xb0C), temp adv (\xb0C hr$^{-1}$), and wind barbs (kt)', loc = 'left')
plt.title(f'valid: {vtime} utc', loc = 'right')
# Save figure.
plt.savefig('narr-5-29-2004-12utc-850mb-temp-adv-map.png')
Figure 1. NARR 850-mb output map that indicate geopotential height (m), temperature (℃), temperature advection (℃ hr-1 ), and wind barbs (kt) for 12 UTC on May 29, 2004.
- The "ValueError: 20m is not a valid Natural Earth scale" problem with cartopy_utils.py: here.