Using Python to Access NCEI Archived NEXRAD Level 2 Data

This notebook shows how to access the THREDDS Data Server (TDS) instance that is serving up archived NEXRAD Level 2 data hosted on Amazon S3. The TDS provides a mechanism to query for available data files, as well as provides access to the data as native volume files, through OPeNDAP, and using its own CDMRemote protocol. Since we're using Python, we can take advantage of Unidata's Siphon package, which provides an easy API for talking to THREDDS servers.

NOTE: Due to data charges, the TDS instance in AWS only allows access to .edu domains. For other users interested in using Siphon to access radar data, you can access recent (2 weeks') data by changing the server URL below to: http://thredds.ucar.edu/thredds/radarServer/nexrad/level2/IDD/

But first! Bookmark these resources for when you want to use Siphon later!

Downloading the single latest volume

Just a bit of initial set-up to use inline figures and quiet some warnings.


In [1]:
import matplotlib
import warnings
warnings.filterwarnings("ignore", category=matplotlib.cbook.MatplotlibDeprecationWarning)
%matplotlib inline

First we'll create an instance of RadarServer to point to the appropriate radar server access URL.


In [2]:
# The S3 URL did not work for me, despite .edu domain
#url = 'http://thredds-aws.unidata.ucar.edu/thredds/radarServer/nexrad/level2/S3/'

#Trying motherlode URL
url = 'http://thredds.ucar.edu/thredds/radarServer/nexrad/level2/IDD/'
from siphon.radarserver import RadarServer
rs = RadarServer(url)

Next, we'll create a new query object to help request the data. Using the chaining methods, let's ask for the latest data at the radar KLVX (Louisville, KY). We see that when the query is represented as a string, it shows the encoded URL.


In [3]:
from datetime import datetime, timedelta
query = rs.query()
query.stations('KLVX').time(datetime.utcnow())


Out[3]:
time=2016-06-20T16%3A50%3A48.016854&stn=KLVX

We can use the RadarServer instance to check our query, to make sure we have required parameters and that we have chosen valid station(s) and variable(s)


In [4]:
rs.validate_query(query)


Out[4]:
True

Make the request, which returns an instance of TDSCatalog; this handles parsing the returned XML information.


In [5]:
catalog = rs.get_catalog(query)

We can look at the datasets on the catalog to see what data we found by the query. We find one volume in the return, since we asked for the volume nearest to a single time.


In [6]:
catalog.datasets


Out[6]:
OrderedDict([('Level2_KLVX_20160620_1633.ar2v',
              <siphon.catalog.Dataset at 0x7f8bf07bb590>)])

We can pull that dataset out of the dictionary and look at the available access URLs. We see URLs for OPeNDAP, CDMRemote, and HTTPServer (direct download).


In [7]:
ds = list(catalog.datasets.values())[0]
ds.access_urls


Out[7]:
{'CdmRemote': 'http://thredds.ucar.edu/thredds/cdmremote/nexrad/level2/IDD/KLVX/20160620/Level2_KLVX_20160620_1633.ar2v',
 'HTTPServer': 'http://thredds.ucar.edu/thredds/fileServer/nexrad/level2/IDD/KLVX/20160620/Level2_KLVX_20160620_1633.ar2v',
 'OPENDAP': 'http://thredds.ucar.edu/thredds/dodsC/nexrad/level2/IDD/KLVX/20160620/Level2_KLVX_20160620_1633.ar2v'}

We'll use the CDMRemote reader in Siphon and pass it the appropriate access URL.


In [8]:
from siphon.cdmr import Dataset
data = Dataset(ds.access_urls['CdmRemote'])

We define some helper functions to make working with the data easier. One takes the raw data and converts it to floating point values with the missing data points appropriately marked. The other helps with converting the polar coordinates (azimuth and range) to Cartesian (x and y).


In [9]:
import numpy as np
def raw_to_masked_float(var, data):
    # Values come back signed. If the _Unsigned attribute is set, we need to convert
    # from the range [-127, 128] to [0, 255].
    if var._Unsigned:
        data = data & 255

    # Mask missing points
    data = np.ma.array(data, mask=data==0)

    # Convert to float using the scale and offset
    return data * var.scale_factor + var.add_offset

def polar_to_cartesian(az, rng):
    az_rad = np.deg2rad(az)[:, None]
    x = rng * np.sin(az_rad)
    y = rng * np.cos(az_rad)
    return x, y

The CDMRemote reader provides an interface that is almost identical to the usual python NetCDF interface. We pull out the variables we need for azimuth and range, as well as the data itself.


In [10]:
sweep = 0
ref_var = data.variables['Reflectivity_HI']
ref_data = ref_var[sweep]
rng = data.variables['distanceR_HI'][:]
az = data.variables['azimuthR_HI'][sweep]

Then convert the raw data to floating point values and the polar coordinates to Cartesian.


In [11]:
ref = raw_to_masked_float(ref_var, ref_data)
x, y = polar_to_cartesian(az, rng)

MetPy is a Python package for meteorology (Documentation: http://metpy.readthedocs.org and GitHub: http://github.com/MetPy/MetPy). We import MetPy and use it to get the colortable and value mapping information for the NWS Reflectivity data.


In [12]:
from metpy.plots import ctables  # For NWS colortable
ref_norm, ref_cmap = ctables.registry.get_with_steps('NWSReflectivity', 5, 5)

Finally, we plot them up using matplotlib and cartopy. We create a helper function for making a map to keep things simpler later.


In [13]:
import matplotlib.pyplot as plt
import cartopy

def new_map(fig, lon, lat):
    # Create projection centered on the radar. This allows us to use x
    # and y relative to the radar.
    proj = cartopy.crs.LambertConformal(central_longitude=lon, central_latitude=lat)

    # New axes with the specified projection
    ax = fig.add_subplot(1, 1, 1, projection=proj)

    # Add coastlines
    ax.coastlines('50m', 'black', linewidth=2, zorder=2)

    # Grab state borders
    state_borders = cartopy.feature.NaturalEarthFeature(
        category='cultural', name='admin_1_states_provinces_lines',
        scale='50m', facecolor='none')
    ax.add_feature(state_borders, edgecolor='black', linewidth=1, zorder=3)
    
    return ax

Download a collection of historical data

This time we'll make a query based on a longitude, latitude point and using a time range.


In [14]:
query = rs.query()
#dt = datetime(2012, 10, 29, 15) # Our specified time
dt = datetime(2016, 6, 8, 18) # Our specified time
query.lonlat_point(-73.687, 41.175).time_range(dt, dt + timedelta(hours=1))


Out[14]:
time_start=2016-06-08T18%3A00%3A00&time_end=2016-06-08T19%3A00%3A00&latitude=41.175&longitude=-73.687

The specified longitude, latitude are in NY and the TDS helpfully finds the closest station to that point. We can see that for this time range we obtained multiple datasets.


In [15]:
cat = rs.get_catalog(query)
cat.datasets


Out[15]:
OrderedDict([('Level2_KOKX_20160608_1828.ar2v',
              <siphon.catalog.Dataset at 0x7f8bec11c7d0>),
             ('Level2_KOKX_20160608_1801.ar2v',
              <siphon.catalog.Dataset at 0x7f8bec11c850>),
             ('Level2_KOKX_20160608_1821.ar2v',
              <siphon.catalog.Dataset at 0x7f8bec11c810>),
             ('Level2_KOKX_20160608_1814.ar2v',
              <siphon.catalog.Dataset at 0x7f8bec11c890>),
             ('Level2_KOKX_20160608_1835.ar2v',
              <siphon.catalog.Dataset at 0x7f8bec11c8d0>),
             ('Level2_KOKX_20160608_1841.ar2v',
              <siphon.catalog.Dataset at 0x7f8bec11c910>),
             ('Level2_KOKX_20160608_1848.ar2v',
              <siphon.catalog.Dataset at 0x7f8bec11c950>),
             ('Level2_KOKX_20160608_1854.ar2v',
              <siphon.catalog.Dataset at 0x7f8bec11c990>),
             ('Level2_KOKX_20160608_1807.ar2v',
              <siphon.catalog.Dataset at 0x7f8bec11c9d0>)])

Grab the first dataset so that we can get the longitude and latitude of the station and make a map for plotting. We'll go ahead and specify some longitude and latitude bounds for the map.


In [16]:
ds = list(cat.datasets.values())[0]
data = Dataset(ds.access_urls['CdmRemote'])
# Pull out the data of interest
sweep = 0
rng = data.variables['distanceR_HI'][:]
az = data.variables['azimuthR_HI'][sweep]
ref_var = data.variables['Reflectivity_HI']

# Convert data to float and coordinates to Cartesian
ref = raw_to_masked_float(ref_var, ref_var[sweep])
x, y = polar_to_cartesian(az, rng)

Use the function to make a new map and plot a colormapped view of the data


In [ ]:
fig = plt.figure(figsize=(10, 10))
ax = new_map(fig, data.StationLongitude, data.StationLatitude)

# Set limits in lat/lon space
ax.set_extent([-77, -70, 38, 42])

# Add ocean and land background
ocean = cartopy.feature.NaturalEarthFeature('physical', 'ocean', scale='50m',
                                            edgecolor='face',
                                            facecolor=cartopy.feature.COLORS['water'])
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='50m',
                                           edgecolor='face',
                                           facecolor=cartopy.feature.COLORS['land'])

ax.add_feature(ocean, zorder=-1)
ax.add_feature(land, zorder=-1)
#ax = new_map(fig, data.StationLongitude, data.StationLatitude)
ax.pcolormesh(x, y, ref, cmap=ref_cmap, norm=ref_norm, zorder=0);

Now we can loop over the collection of returned datasets and plot them. As we plot, we collect the returned plot objects so that we can use them to make an animated plot. We also add a timestamp for each plot.


In [ ]:
meshes = []
for item in sorted(cat.datasets.items()):
    # After looping over the list of sorted datasets, pull the actual Dataset object out
    # of our list of items and access over CDMRemote
    ds = item[1]
    data = Dataset(ds.access_urls['CdmRemote'])

    # Pull out the data of interest
    sweep = 0
    rng = data.variables['distanceR_HI'][:]
    az = data.variables['azimuthR_HI'][sweep]
    ref_var = data.variables['Reflectivity_HI']

    # Convert data to float and coordinates to Cartesian
    ref = raw_to_masked_float(ref_var, ref_var[sweep])
    x, y = polar_to_cartesian(az, rng)

    # Plot the data and the timestamp
    mesh = ax.pcolormesh(x, y, ref, cmap=ref_cmap, norm=ref_norm, zorder=0)
    text = ax.text(0.65, 0.03, data.time_coverage_start, transform=ax.transAxes,
                   fontdict={'size':16})
    
    # Collect the things we've plotted so we can animate
    meshes.append((mesh, text))

Using matplotlib, we can take a collection of Artists that have been plotted and turn them into an animation. With matplotlib 1.5 (1.5-rc2 is available now!), this animation can be converted to HTML5 video viewable in the notebook.


In [ ]:
# Set up matplotlib to do the conversion to HTML5 video
import matplotlib
matplotlib.rcParams['animation.html'] = 'html5'

# Create an animation
from matplotlib.animation import ArtistAnimation
ArtistAnimation(fig, meshes)



IOErrorTraceback (most recent call last)
/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/IPython/core/formatters.pyc in __call__(self, obj)
    341             method = _safe_get_formatter_method(obj, self.print_method)
    342             if method is not None:
--> 343                 return method()
    344             return None
    345         else:

/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/matplotlib/animation.pyc in _repr_html_(self)
    979         fmt = rcParams['animation.html']
    980         if fmt == 'html5':
--> 981             return self.to_html5_video()
    982 
    983 

/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/matplotlib/animation.pyc in to_html5_video(self)
    953                                 bitrate=rcParams['animation.bitrate'],
    954                                 fps=1000. / self._interval)
--> 955                 self.save(f.name, writer=writer)
    956 
    957             # Now open and base64 encode

/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/matplotlib/animation.pyc in save(self, filename, writer, fps, dpi, codec, bitrate, extra_args, metadata, extra_anim, savefig_kwargs)
    808                     # TODO: Need to see if turning off blit is really necessary
    809                     anim._draw_next_frame(d, blit=False)
--> 810                 writer.grab_frame(**savefig_kwargs)
    811 
    812         # Reconnect signal for first draw if necessary

/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/matplotlib/animation.pyc in grab_frame(self, **savefig_kwargs)
    228             # frame format and dpi.
    229             self.fig.savefig(self._frame_sink(), format=self.frame_format,
--> 230                              dpi=self.dpi, **savefig_kwargs)
    231         except RuntimeError:
    232             out, err = self._proc.communicate()

/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/matplotlib/figure.pyc in savefig(self, *args, **kwargs)
   1563             self.set_frameon(frameon)
   1564 
-> 1565         self.canvas.print_figure(*args, **kwargs)
   1566 
   1567         if frameon:

/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/matplotlib/backend_bases.pyc in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, **kwargs)
   2230                 orientation=orientation,
   2231                 bbox_inches_restore=_bbox_inches_restore,
-> 2232                 **kwargs)
   2233         finally:
   2234             if bbox_inches and restore_bbox:

/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in print_raw(self, filename_or_obj, *args, **kwargs)
    517             close = False
    518         try:
--> 519             fileobj.write(renderer._renderer.buffer_rgba())
    520         finally:
    521             if close:

IOError: [Errno 32] Broken pipe
Out[ ]:
<matplotlib.animation.ArtistAnimation at 0x7f8bec11cc90>

In [ ]: