Using Python to get the latest NEXRAD Composite

In [ ]:

    

# Set-up for notebook
%matplotlib inline

# Some needed imports
import datetime as dt
import matplotlib.pyplot as plt
import matplotlib as mpl
import cartopy
import numpy as np
from netCDF4 import Dataset
from pyudl.tds import TDSCatalog


    

In [ ]:

    

# Function to create a colortable that matches the NWS colortable
def radar_colormap():
    nws_reflectivity_colors = [
    "#646464", # ND
    "#ccffff", # -30
    "#cc99cc", # -25
    "#996699", # -20
    "#663366", # -15
    "#cccc99", # -10
    "#999966", # -5
    "#646464", # 0
    "#04e9e7", # 5
    "#019ff4", # 10
    "#0300f4", # 15
    "#02fd02", # 20
    "#01c501", # 25
    "#008e00", # 30
    "#fdf802", # 35
    "#e5bc00", # 40
    "#fd9500", # 45
    "#fd0000", # 50
    "#d40000", # 55
    "#bc0000", # 60
    "#f800fd", # 65
    "#9854c6", # 70
    "#fdfdfd" # 75
    ]

    return mpl.colors.ListedColormap(nws_reflectivity_colors)


    

In [ ]:

    

# Get today's date...
today = dt.datetime.utcnow()

# ...and use that to assemble the URL and grab the catalog
url = "http://thredds.ucar.edu/thredds/catalog/nexrad/composite/gini/n0r/1km/{}/catalog.xml".format(today.strftime("%Y%m%d"))
cat = TDSCatalog(url)

# Get the list of names of datasets and sort, so that the last dataset is the latest
names = cat.datasets.keys()
names.sort()
latest = names[-1]

# Grab the dataset for the latest
latestDs = cat.datasets[latest]


    

In [ ]:

    

# Construct a NetCDF dataset using the OPeNDAP access URL
dataset = Dataset(latestDs.accessUrls['OPENDAP'])
print dataset.variables.keys()


    

In [ ]:

    

print dir(dataset.variables['LambertConformal'])


    

In [ ]:

    

##################
# Projection fun #
##################

# get basic info from the file regarding projection attributes
# see the following for more info on projections as described here:
#   http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html#GRID218
#   http://www.wmo.int/pages/prog/www/WDM/Guides/Guide-binary-2.html [see LAMBERT CONFORMAL SECANT OR TANGENT CONE GRIDS]
#   http://www.unidata.ucar.edu/mailing_lists/archives/netcdf-java/2006/msg00006.html [starndard parallels in CDM]
proj = dataset.variables['LambertConformal']
rsphere = proj.earth_radius

# lat_0	: center of desired map domain (in degrees) [Basemap]
# CDM : 'latitude_of_projection_origin'
lat_0 = proj.latitude_of_projection_origin

# lon_0	: center of desired map domain (in degrees) [Basemap]
# CDM : 'longitude_of_central_meridian'
lon_0 = proj.longitude_of_central_meridian

# lat_1, lat_2 : 1st and second parallels [Basemap]
# CDM : 'standard_parallel' - this attr contains both 1st and 2nd
lat_1 = proj.standard_parallel
print lat_0, lon_0, lat_1, rsphere

# hardcoded from WMS/WCS metadata - need to add metadata into catalog class
# to grab this programatically
llcrnrlon = -129.2192094372674
llcrnrlat = 21.742620467908093
urcrnrlon = -63.559241752667674
urcrnrlat = 49.126151231865485


    

In [ ]:

    

# Used to subset the data
xstride = 10
ystride = 10

# download x and y coords and convert them from km to m
x = dataset.variables['x'][::xstride] * 1000.
y = dataset.variables['y'][::ystride] * 1000.

# Grab the reflectivity data. Mask values less than -30 dBz
data = dataset.variables["Reflectivity"][0, 0::ystride, 0::xstride]
data = np.ma.array(data, mask=data<=-30)


    

In [ ]:

    

# Set up the projection for the LambertConformal projection we know we have
lcc = cartopy.crs.LambertConformal(central_longitude=lon_0, central_latitude=lat_0, secant_latitudes=(lat_0, lat_1))

# Create a large figure and axes with this projection
fig = plt.figure(figsize=(24, 12))
ax = fig.add_subplot(1, 1, 1, projection=lcc)

# Limit to the bounds of the data we have
ax.set_extent([llcrnrlon, urcrnrlon, llcrnrlat, urcrnrlat], cartopy.crs.Geodetic())

# Add some map features
ax.stock_img()
ax.coastlines(resolution='50m')
ax.add_feature(cartopy.feature.NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lines', scale='50m',
                                                   facecolor='none'))
ax.add_feature(cartopy.feature.BORDERS, linewidth='2', edgecolor='black')
ax.gridlines()

# Convert the time to text and add as the title
time = dataset.variables["time"][:][0] / 1000.
title = dt.datetime.fromtimestamp(time).isoformat()
ax.set_title(title)

# Plot the data as an image, using the x and y values we have as the extents of the image
# NOTE: This assumes equal-spaced points
cmap = radar_colormap()
norm = mpl.colors.Normalize(vmin=-35, vmax=80)
cax = ax.imshow(data, extent=(x.min(), x.max(), y.min(), y.max()), cmap=cmap, norm=norm, origin="upper",
                transform=lcc)
plt.colorbar(cax)