Using Python to get the latest NEXRAD Composite

Unidata Python Workshop


Objective: Visualize the latest available reflectivity data composited data

Steps involved:

  • Construct the appropriate URL to get the latest data file
  • Open the URL using netCDF4-python
  • Read the basic metadata
  • Create the appropriate CartoPy projection and plot the Radar Reflectivity

In [1]:
# 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 siphon.catalog import TDSCatalog
from metpy.plots import ctables

Get the latest data URL, grab the metadata, and request the data


In [2]:
# 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
names = list(cat.datasets.keys())
print(names[:15])


['Level3_Composite_n0r_1km_20180302_1730.gini', 'Level3_Composite_n0r_1km_20180302_1725.gini', 'Level3_Composite_n0r_1km_20180302_1720.gini', 'Level3_Composite_n0r_1km_20180302_1715.gini', 'Level3_Composite_n0r_1km_20180302_1710.gini', 'Level3_Composite_n0r_1km_20180302_1705.gini', 'Level3_Composite_n0r_1km_20180302_1700.gini', 'Level3_Composite_n0r_1km_20180302_1655.gini', 'Level3_Composite_n0r_1km_20180302_1650.gini', 'Level3_Composite_n0r_1km_20180302_1645.gini', 'Level3_Composite_n0r_1km_20180302_1640.gini', 'Level3_Composite_n0r_1km_20180302_1635.gini', 'Level3_Composite_n0r_1km_20180302_1630.gini', 'Level3_Composite_n0r_1km_20180302_1625.gini', 'Level3_Composite_n0r_1km_20180302_1620.gini']

In [3]:
# sort, so that the last dataset is the latest
names.sort()
latest = names[-1]
print(latest)


Level3_Composite_n0r_1km_20180302_1730.gini

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

In [5]:
# Construct a NetCDF dataset using the OPeNDAP access URL
dataset = Dataset(latestDs.access_urls['OPENDAP'])
print(list(dataset.variables.keys()))


['time', 'x', 'y', 'LambertConformal', 'Reflectivity']

In [6]:
dataset.variables['LambertConformal'].ncattrs()


Out[6]:
['grid_mapping_name',
 'latitude_of_projection_origin',
 'longitude_of_central_meridian',
 'standard_parallel',
 'earth_radius',
 '_CoordinateTransformType',
 '_CoordinateAxes']

In [7]:
##################
# 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
#   [standard 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)


40.0 -100.0 40.0 6371229.0

Grab the data


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)

Create the Plot


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,
                                   standard_parallels=(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([-129., -63., 22., 49.], cartopy.crs.Geodetic())

# Add some map features
ax.stock_img()
ax.coastlines(resolution='50m')
ax.add_feature(cartopy.feature.STATES)
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 = ctables.registry.get_colortable('NWSReflectivityExpanded')
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);

Exercise

Using what was done above, plot the digital hybrid reflectivity (DHR):


In [ ]: