The Chicago area saw significant flooding on April 18, 2013. The NWS 1-Day Observed Precipitation map for that date shows significant precipitation (3+ inches) over much of northern Illinois. Recently a collegue asked me for help in creating a high resolution plot of this data for a poster. Below is instructions on how to make a similar map using Python.
The NWS makes its precipitation data available for download in NetCDF format, specifically we wanted to look at the April 18, 2013 data. Two files are found in the tar.gz file, nws_precip_conus_20130418.nc and nws_precip_pr_20130418.nc for the continential US and Puerto Rico, we are interested in the continential US data.
The netCDF file contains a grid of the 24-hour precipitation in hundredths of a millimeters but, as noted on the download page, the latitude and longitude of the grid is not explicitly stated in the file. Luckily the NWS provides a C-program to convert the netCDF file into a ACSII file which contain the latitude and longitude. Running this program on the data a 20130418.txt file is created. We can read in this file using NumPy's genfromtxt function and extract out the precipitation amount, latitudes and longitudes.
%pylab inline
import numpy as np
data = np.genfromtxt('20130418.txt', names=True, delimiter=',')
lats = data['lat'].copy().reshape(813, 1051)
lons = data['long'].copy().reshape(813, 1051)
precip = data['value'].copy().reshape(813, 1051)
precip_in = np.ma.masked_less(precip, 0.01)
Before we return to the data we need to plot a map of the continential USA. The Basemap package makes this easy, the newer Cartopy package offers similar features. We will make the a function since we will be reusing it often.
from mpl_toolkits.basemap import Basemap
def plot_conus():
m = Basemap(
llcrnrlon=-130.0,
llcrnrlat=20.0,
urcrnrlon=-60.0,
urcrnrlat=55.0,
projection='mill',
resolution='c')
m.drawcoastlines()
m.drawcountries()
m.drawstates()
return m
plot_conus()
Now we can plot the precipitation data on this map.
m = plot_conus()
cax = m.pcolormesh(lons, lats, precip_in, latlon=True)
m.colorbar(cax)
We would really like the rain rate colors more in line with those in the original image. Luckily the folks at Unidata have created a NWS colormap, modifying this colormap slightly.
nws_precip_colors = [
"#04e9e7", # 0.01 - 0.10 inches
"#019ff4", # 0.10 - 0.25 inches
"#0300f4", # 0.25 - 0.50 inches
"#02fd02", # 0.50 - 0.75 inches
"#01c501", # 0.75 - 1.00 inches
"#008e00", # 1.00 - 1.50 inches
"#fdf802", # 1.50 - 2.00 inches
"#e5bc00", # 2.00 - 2.50 inches
"#fd9500", # 2.50 - 3.00 inches
"#fd0000", # 3.00 - 4.00 inches
"#d40000", # 4.00 - 5.00 inches
"#bc0000", # 5.00 - 6.00 inches
"#f800fd", # 6.00 - 8.00 inches
"#9854c6", # 8.00 - 10.00 inches
"#fdfdfd" # 10.00+
]
precip_colormap = matplotlib.colors.ListedColormap(nws_precip_colors)
Recreating the plot using this new colormap with the correct precipitation level boundries.
m = plot_conus()
levels = [0.01, 0.1, 0.25, 0.50, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0,
6.0, 8.0, 10., 20.0]
norm = matplotlib.colors.BoundaryNorm(levels, 15)
cax = m.pcolormesh(lons, lats, precip_in, latlon=True, norm=norm,
cmap=precip_colormap)
m.colorbar(cax)
This plot looks very similar to the one provided by the NWS. But we had to compile a C program and make a very large (50+ Meg) ASCII file, can't we determine the latitudes and longitudes in Python and read the data from the netCDF file? Yes we can, translating the relevant the portion of the C file into Python we can create grids of latitude and longitudes. Note that this code could be vectorized to run much faster if better performance is needed.
HRAP_XOR = 14
HRAP_YOR = 10
def lat_lon_from_hrap(hrap_x, hrap_y):
raddeg = 57.29577951
earthrad = 6371.2
stdlon = 105.
mesh_len = 4.7625
tlat = 60. / raddeg
x = hrap_x - 401.
y = hrap_y - 1601.
rr = x * x + y * y
gi = ((earthrad * (1 + np.sin(tlat))) / mesh_len)
gi = gi * gi
ll_y = np.arcsin((gi - rr) / (gi + rr)) * raddeg
ang = np.arctan2(y, x) * raddeg
if (ang < 0):
ang = ang+360.
ll_x = 270 + stdlon - ang
if (ll_x < 0):
ll_x = ll_x + 360
if (ll_x > 360):
ll_x = ll_x - 360
return ll_x, ll_y
lats = np.empty((813, 1051), dtype='float')
lons = np.empty((813, 1051), dtype='float')
for i in xrange(813):
for j in xrange(1051):
hrap_x = j + HRAP_XOR + 0.5
hrap_y = i + HRAP_YOR + 0.5
lon, lat = lat_lon_from_hrap(hrap_x, hrap_y)
lats[i,j] = lat
lons[i,j] = -lon
Now the data precipitation data can be read in from the netCDF file and converted to precipitation in inches.
import netCDF4
dset = netCDF4.Dataset('nws_precip_conus_20130418.nc')
precip = dset.variables['amountofprecip'][:]
precip_in = np.ma.masked_less(precip / 2540., 0.01)
Replotting this data
m = plot_conus()
levels = [0.01, 0.1, 0.25, 0.50, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0,
6.0, 8.0, 10., 20.0]
norm = matplotlib.colors.BoundaryNorm(levels, 15)
cax = m.pcolormesh(lons, lats, precip_in, latlon=True, norm=norm,
cmap=precip_colormap)
m.colorbar(cax)
We can putting this all together into a script that will create a high resolution png file for a poster or presentation.
#! /usr/bin/env python
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import netCDF4
def plot_conus():
""" Create a plot of the Continential US. """
m = Basemap(
llcrnrlon=-130.0,
llcrnrlat=20.0,
urcrnrlon=-60.0,
urcrnrlat=55.0,
projection='mill',
resolution='c')
m.drawcoastlines()
m.drawcountries()
m.drawstates()
return m
# Colorbar with NSW Precip colors
nws_precip_colors = [
"#04e9e7", # 0.01 - 0.10 inches
"#019ff4", # 0.10 - 0.25 inches
"#0300f4", # 0.25 - 0.50 inches
"#02fd02", # 0.50 - 0.75 inches
"#01c501", # 0.75 - 1.00 inches
"#008e00", # 1.00 - 1.50 inches
"#fdf802", # 1.50 - 2.00 inches
"#e5bc00", # 2.00 - 2.50 inches
"#fd9500", # 2.50 - 3.00 inches
"#fd0000", # 3.00 - 4.00 inches
"#d40000", # 4.00 - 5.00 inches
"#bc0000", # 5.00 - 6.00 inches
"#f800fd", # 6.00 - 8.00 inches
"#9854c6", # 8.00 - 10.00 inches
"#fdfdfd" # 10.00+
]
precip_colormap = matplotlib.colors.ListedColormap(nws_precip_colors)
# latitude and longitudes of HRAP grid
HRAP_XOR = 14
HRAP_YOR = 10
def lat_lon_from_hrap(hrap_x, hrap_y):
""" Calculate the latitude and longitude for a HRAP grid. """
raddeg = 57.29577951
earthrad = 6371.2
stdlon = 105.
mesh_len = 4.7625
tlat = 60. / raddeg
x = hrap_x - 401.
y = hrap_y - 1601.
rr = x * x + y * y
gi = ((earthrad * (1 + np.sin(tlat))) / mesh_len)
gi = gi * gi
ll_y = np.arcsin((gi - rr) / (gi + rr)) * raddeg
ang = np.arctan2(y, x) * raddeg
if (ang < 0):
ang = ang + 360.
ll_x = 270 + stdlon - ang
if (ll_x < 0):
ll_x = ll_x + 360
if (ll_x > 360):
ll_x = ll_x - 360
return ll_x, ll_y
lats = np.empty((813, 1051), dtype='float')
lons = np.empty((813, 1051), dtype='float')
for i in xrange(813):
for j in xrange(1051):
hrap_x = j + HRAP_XOR + 0.5
hrap_y = i + HRAP_YOR + 0.5
lon, lat = lat_lon_from_hrap(hrap_x, hrap_y)
lats[i, j] = lat
lons[i, j] = -lon
# read in the data, convert in inches
dset = netCDF4.Dataset('nws_precip_conus_20130418.nc')
precip = dset.variables['amountofprecip'][:]
precip_in = np.ma.masked_less(precip / 2540., 0.01)
# plot the data
m = plot_conus()
levels = [0.01, 0.1, 0.25, 0.50, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0,
6.0, 8.0, 10., 20.0]
norm = matplotlib.colors.BoundaryNorm(levels, 15)
cax = m.pcolormesh(lons, lats, precip_in, latlon=True, norm=norm,
cmap=precip_colormap)
m.colorbar(cax)
plt.savefig('conus_plot.png', dpi=200)