Objective: Visualize the latest available reflectivity data composited data
Steps involved:
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
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])
In [3]:
# sort, so that the last dataset is the latest
names.sort()
latest = names[-1]
print(latest)
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()))
In [6]:
dataset.variables['LambertConformal'].ncattrs()
Out[6]:
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)
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,
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);
Using what was done above, plot the digital hybrid reflectivity (DHR):
In [ ]: