Using Python to get the latest Level II RADAR for a station

Objective: Visualize the latest available reflectivity data from the lowest tilt of a user specified radar site

Steps involved:

Define a color map based on the one used by the National Weather Service
Construct the appropriate URL to get the latest data file
Open the URL using PyDAP
Read the basic metadata
Find the lowest tilt and construct a coordinate system
Create the appropriate Basemap projection and plot the Radar Reflectivity



In [1]:

    
import string
import numpy as np
import matplotlib.pyplot as plt
import datetime as dt
import matplotlib as mpl
from pydap.client import open_url
from mpl_toolkits.basemap import Basemap
from pytds.util import get_latest_dods_url

Define the NWS Color Map for reflectivity

While matplotlib has some really nice, built-in colormaps, nothing quite matches the NWS colormap
Python List of HTML colors (used photoshop color picker on a screenshot of the colorbar off of the NWS radar page)
Converted into matplotlib colormap using ListedColormap



In [2]:

    
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
    ]

    cmap = mpl.colors.ListedColormap(nws_reflectivity_colors)

    return cmap

Which radar site do you want to view?



In [3]:

    
site = "KEVX"

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



In [4]:

    
#
# get latest data url
#

date = dt.date.today().strftime("%Y%m%d")
url = "http://thredds.ucar.edu/thredds/catalog/nexrad/level2/" + site  + "/" + date + "/catalog.xml"
latest_data_url = get_latest_dods_url(url)

# open url using pydap
dataset = open_url(latest_data_url)

# get some basic info from the global metadata
global_attrs = dataset.attributes['NC_GLOBAL']
station_lat = global_attrs["StationLatitude"]
station_lon = global_attrs["StationLongitude"]
station_id = global_attrs["StationName"]

# get data array and metadata
variable = "Reflectivity_HI"
tilts = dataset["elevationR_HI"][:,0]
scan = np.where(tilts == tilts.min())[0][0]
tilt = tilts[scan,0]
data = dataset[variable][scan,::].squeeze()
theta = dataset["azimuthR_HI"][scan,::].squeeze()
r = dataset["distanceR_HI"][:]
data_attrs = dataset[variable].attributes

Apply scale, offest, and mask radar data



In [5]:

    
# check for scale and offset
if data_attrs.has_key("scale_factor"):
    data = data * data_attrs["scale_factor"]

if data_attrs.has_key("add_offset"):
    data = data + data_attrs["add_offset"]
    
# mask the data based on NWS threshold
data = np.ma.masked_array(data, data < 5)

Compute x,y coord system from the native RADAR polar coordinates.



In [6]:

    
theta = theta * np.pi / 180.
d = r * np.cos(tilt * np.pi / 180.)
x = np.array((np.matrix(np.cos(theta)).transpose() * np.matrix(d)))
y = np.array((np.matrix(np.sin(theta)).transpose() * np.matrix(d)))
height = np.abs(x.min()) + np.abs(x.max())
width = np.abs(y.min()) + np.abs(y.max())
x = x + height / 2.
y = y + width / 2.

Create the plot



In [7]:

    
fig = plt.figure(figsize=(8,11))
ax = fig.add_subplot(1,1,1)

m=Basemap(lat_0=station_lat,lon_0=station_lon,resolution='i',
          projection='laea',height=height,width=width,ax=ax)

station_x, station_y = m(station_lon, station_lat)

cmap = radar_colormap()
norm = mpl.colors.Normalize(vmin=-35, vmax=80)

ax.text(station_x, station_y, "+{}".format(site))

cax = m.pcolormesh(y, x, data, cmap=cmap, norm=norm)
m.drawcoastlines()
m.drawstates()
m.drawcountries()
cbar = m.colorbar(cax)
cbar.set_label(data_attrs["units"])
ax.set_title(station_id)
ax.set_xlim(station_x-250000,station_x+250000)
ax.set_ylim(station_y-250000,station_y+250000)
plt.show()

Exercise:

Change the Radar site (use the NWS site as a guide to obtain the radar site name)
Visit the THREDDS catalog page, or print out the pydap dataset, and try to visualize a new variable
How do our Radar maps at high latitudes compare to those of the NWS? Why do you think this is the case?