Spam, Eggs and Hydrometeors

A blog on Python and Atmospheric science

Plotting NWS Precipitation Data

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.

In [1]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib

In [2]:
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.

In [3]:
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()
Out[3]:
<mpl_toolkits.basemap.Basemap at 0x11dfe3110>

Now we can plot the precipitation data on this map.

In [4]:
m = plot_conus()
cax = m.pcolormesh(lons, lats, precip_in, latlon=True)
m.colorbar(cax)
Out[4]:
<matplotlib.colorbar.Colorbar instance at 0x120f0a878>

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.

In [5]:
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.

In [6]:
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)
Out[6]:
<matplotlib.colorbar.Colorbar instance at 0x11dffbea8>

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.

In [7]:
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.

In [8]:
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

In [9]:
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)
Out[9]:
<matplotlib.colorbar.Colorbar instance at 0x120fcd3f8>

We can putting this all together into a script that will create a high resolution png file for a poster or presentation.

In [10]:
#! /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)

This blog post was written entirely in the IPython Notebook. The full notebook can be downloaded here, or viewed statically here.