In [143]:
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 siphon.catalog import get_latest_dods_url
from matplotlib.colorbar import ColorbarBase
%matplotlib inline

In [144]:
from IPython.html import widgets
from IPython.display import clear_output, display, HTML

In [145]:
display_dict = {'Reflectivity' : (plt.get_cmap('jet'), mpl.colors.Normalize(vmin=-35, vmax=80)),
                'RadialVelocity' : (plt.get_cmap('jet'), mpl.colors.Normalize(vmin=-30, vmax=30)),
                'SpectrumWidth' : (plt.get_cmap('jet'), mpl.colors.Normalize(vmin=0, vmax=30)),
                'DifferentialReflectivity' : (plt.get_cmap('jet'), mpl.colors.Normalize(vmin=-5, vmax=10)),
                'DifferentialPhase' : (plt.get_cmap('jet'), mpl.colors.Normalize(vmin=0, vmax=100)),
                'CorrelationCoefficient' : (plt.get_cmap('jet'), mpl.colors.Normalize(vmin=0, vmax=1.0))}
def lookup_display(var):
    var = var.split('_')[0]
    return display_dict.get(var, (plt.get_cmap('jet'), mpl.colors.Normalize(vmin=-35, vmax=80)))

In [158]:
class RadarPlotter(object):
    def __init__(self):
        self.data_plot = None
        self.on_vars = []
        self.on_tilts = []
        self.setSite('KFTG')
        self.setVar('Reflectivity')
        self.setData(0)

    def newFigure(self):
        self._fig = plt.figure(figsize=(12, 8))
        self._ax = self._fig.add_axes((0.08, 0.08, 0.8, 0.9), label="1")
        self._cax = self._fig.add_axes((0.92, 0.08, 0.05, 0.9), label="c1")
        self.cbar = ColorbarBase(ax=self._cax, orientation='vertical')
        
    def setSite(self, siteID):
        self.newFigure()
        self.siteID = siteID
        self.updateDataset()

        # Set up the projection for this site
        self.proj = Basemap(lat_0=self.site_lat, lon_0=self.site_lon, resolution='l',
                            projection='laea', height=460000., width=460000., ax=self._ax)
        self.site_x, self.site_y = self.proj(self.site_lon, self.site_lat)
        self.proj.drawstates()
        self.proj.drawcounties()
        
        self._ax.text(self.site_x, self.site_y, "+{}".format(self.siteID))
        
    def updateDataset(self):
        self.fetchData()
        self.parseDataset()
        for r in self.on_vars:
            r(self.variables, self.var)

    def fetchData(self, date=None):
        #
        # get latest data url
        #
        if not date:
            date = dt.date.today().strftime("%Y%m%d")
        url = "http://thredds.ucar.edu/thredds/catalog/nexrad/level2/" + self.siteID  + "/" + date + "/catalog.xml"
        latest_data_url = get_latest_dods_url(url)
        
        # open url using pydap
        self.dataset = open_url(latest_data_url)

    def parseDataset(self):
        # get some basic info from the global metadata
        global_attrs = self.dataset.attributes['NC_GLOBAL']
        self.site_lat = global_attrs["StationLatitude"]
        self.site_lon = global_attrs["StationLongitude"]
        self.site_name = global_attrs["StationName"]
        self.time = global_attrs["time_coverage_end"]
        self.variables = [var for var in self.dataset.keys() if len(self.dataset[var].dimensions) == 3]
        self.variables = list(set(var[:-3] if var.endswith('_HI') else var for var in self.variables))

    def setVar(self, varName):
        self.var = varName
        self.cmap, self.norm = lookup_display(varName)
        self.cbar.set_norm(self.norm)
        self.cbar.set_cmap(self.cmap)
        self.cbar.draw_all()

        suffix = self.getSuffix(varName)
        elVar = 'elevation' + suffix
        
        if varName + '_HI' in self.dataset.keys():
            self.tilts = self.dataset[elVar + '_HI'][:, 0].squeeze().tolist()
        else:
            self.tilts = []
        self.hi_count = len(self.tilts)
        self.tilts.extend(self.dataset[elVar][:, 0].squeeze().tolist())
        for r in self.on_tilts:
            r(self.tilts, self.cur_scan)

    def getSuffix(self, varName):
        # Figure out what variables to use for spatial parameters
        if 'Velocity' in varName or 'Spectrum' in varName:
            suffix = 'V'
        else:
            suffix = varName[0]
        return suffix
        
    def setData(self, tilt):
        varName = self.var
        suffix = self.getSuffix(varName)
        self.cur_el = self.tilts[tilt]
        if tilt < self.hi_count:
            varName += '_HI'
            suffix += '_HI'
        else:
            tilt -= self.hi_count

        self.cur_scan = tilt
        
        azVar = 'azimuth' + suffix
        rngVar = 'distance' + suffix

        self.theta = self.dataset[azVar][tilt, ::].squeeze()
        self.rng = self.dataset[rngVar][:]
        theta = np.deg2rad(self.theta)
        d = self.rng * np.cos(np.deg2rad(self.cur_el))
        self.x = d * np.sin(theta[:, np.newaxis]) + self.site_x
        self.y = d * np.cos(theta[:, np.newaxis]) + self.site_y

        #cbar = self.proj.colorbar(cax)
        #cbar.set_label(data_attrs["units"])

        # Get data and metadata
        data = self.dataset[varName][self.cur_scan, ::].squeeze()
        data_attrs = self.dataset[varName].attributes
        
        # check for missing values
        mask = np.zeros_like(data)
        if 'missing_value' in data_attrs:
            for val in data_attrs['missing_value']:
                mask = mask | (data == val)
        data = np.ma.array(data, mask=mask)

        # 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"]
        
        self.data = data
        
    def drawData(self):
        if self.data_plot and self.data_plot in self._ax.collections:
            self.data_plot.remove()
        self.data_plot = self.proj.pcolormesh(self.x, self.y, self.data, cmap=self.cmap, norm=self.norm)
        self._ax.set_title("{0} at {1}".format(self.site_name, self.time))

    def refresh(self):
        self.drawData()
        clear_output(wait=True)
        display(self._fig)

    def updateSite(self, site):
        self.setSite(site)
        self.setVar(self.var)
        self.setData(self.cur_scan)
        self.refresh()
    
    def updateVar(self, var):
        self.setVar(var)
        self.setData(self.cur_scan)
        self.refresh()
    
    def updateTilt(self, tilt):
        self.setData(tilt)
        self.refresh()

In [159]:
rp = RadarPlotter()

siteSelect = widgets.DropdownWidget(values=['KFTG', 'KTLX'], value='KFTG', description='Site:')
varSelect = widgets.DropdownWidget(values=rp.variables, value=rp.var, description='Field:')
tiltSelect = widgets.DropdownWidget(description='Tilt:')
container = widgets.ContainerWidget()
display(container)
container.children = [siteSelect, varSelect, tiltSelect]
container.add_class('hbox')

siteSelect.on_trait_change(lambda n,v:rp.updateSite(v), 'value')
varSelect.on_trait_change(lambda n,v:rp.updateVar(v), 'value')
tiltSelect.on_trait_change(lambda n,v:rp.updateTilt(v), 'value')

def on_vars(newVars, var):
    varSelect.values = newVars
    varSelect.value = var
rp.on_vars.append(on_vars)

def on_tilts(newTilts, tilt):
    tiltSelect.values = range(len(newTilts))
    tiltSelect.labels = ['{:.2f}'.format(n) for n in newTilts]
    tiltSelect.value = tilt
rp.on_tilts.append(on_tilts)
on_tilts(rp.tilts, rp.cur_scan)



In [ ]: