notebook.community

Edit and run



In [8]:

    
import datetime as dt
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
import pytz
from matplotlib import animation
from JSAnimation import IPython_display
from mpl_toolkits.basemap import Basemap
from pydap.client import open_url
from siphon.catalog import TDSCatalog
%matplotlib inline
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



In [24]:

    
today = dt.datetime.utcnow() 
#url = "http://atm.ucar.edu/thredds/catalog/nexrad/composite/gini/n0r/1km/{}/catalog.xml".format(today.strftime("%Y%m%d"))
url = "http://atm.ucar.edu/thredds/catalog/grib/nexrad/composite/unidata/NEXRAD_Unidata_Reflectivity-{}/files/catalog.xml".format(today.strftime("%Y%m%d"))
cat = TDSCatalog(url)
names = cat.datasets.keys()
names.sort()
latest = names[-1]
latestDs = cat.datasets[latest]



In [25]:

    
dataset = open_url(latestDs.accessUrls['OPENDAP'])



In [27]:

    
##################
# 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 [starndard parallels in CDM]
#proj_attributes = dataset['LambertConformal'].attributes
proj_attributes = dataset['LambertConformal_Projection'].attributes
rsphere = proj_attributes['earth_radius']

# lat_0	: center of desired map domain (in degrees) [Basemap]
# CDM : 'latitude_of_projection_origin'
# NCO : where Dx and Dy are defined
# :
lat_0 = proj_attributes['latitude_of_projection_origin']

# lon_0	: center of desired map domain (in degrees) [Basemap]
# CDM : 'longitude_of_central_meridian'
# NCO : Lov
#
# Lov - The orientation of the grid; i.e., the east longitude value
# of the meridian which is parallel to the y-axis (or columns of the
#  grid) along which latitude increases as the y-coordinate increases.
# (Note: The orientation longitude may, or may not, appear within a
# particular grid.)
#
lon_0 = proj_attributes['longitude_of_central_meridian'] # Lov

# lat_1, lat_2 : 1st and second parallels [Basemap]
# CDM : 'standard_parallel' - this attr contains both 1st and 2nd
# NCO : ??? Not sure where this shows up on the nco page
lat_1 = proj_attributes['standard_parallel']
print lat_0, lon_0, lat_1, rsphere

# hardcoded from catalog metadata - will add metadata into catalog class
# to grab this programatically
llcrnrlon = -120.02283
llcrnrlat = 23.0132
urcrnrlon = -63.5105
urcrnrlat = 47.30635









    



40.0 260.0 40.0 6371200.0



In [28]:

    
xstride = 10
ystride = 10

# download x and y coords and convert them from
# km to m, offset for use in basemap
x = dataset['x'][:] * 1000.
y = dataset['y'][:] * 1000.

width = x.max() - x.min()
height = y.max() - y.min()

x = x[0::xstride]
y = y[0::ystride]



In [31]:

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

m = Basemap(projection='lcc', lat_0 = lat_0, lon_0 = lon_0, lat_1 = lat_1,
              llcrnrlon = llcrnrlon, llcrnrlat = llcrnrlat,
              urcrnrlat = urcrnrlat, urcrnrlon = urcrnrlon,
              area_thresh = 1000., rsphere = rsphere, resolution='c')


x = (m.llcrnrx - x.min()) + x
y = (m.llcrnry - y.min()) + y

cbar = None
artists = []
for name in names[-10:]:
    latestDs = cat.datasets[name]
    dataset = open_url(latestDs.accessUrls['OPENDAP'])
    time = dataset["time"][:][0] / 1000.
#    data = np.squeeze(dataset["Reflectivity"][0,0::ystride,0::xstride])
    data = np.squeeze(dataset["Base_reflectivity_surface_layer"][0,0::ystride,0::xstride])
    title = dt.datetime.fromtimestamp(time, pytz.utc).isoformat()
    cmap = radar_colormap()
    norm = mpl.colors.Normalize(vmin=-35, vmax=80)
    title_artist = ax.text(0.4, 1.02, title, transform=ax.transAxes)
    #cax = m.pcolormesh(x, y, data, cmap=cmap, norm=norm)
    cax = m.imshow(data, extent = (x.min(), x.max(), y.min(), y.max()), cmap=cmap, norm=norm, origin="lower")
    artists.append((title_artist, cax))
    m.drawcoastlines()
    m.drawstates()
    m.drawcountries()
    if cbar is None:
        cbar = m.colorbar(cax)

animation.ArtistAnimation(fig, artists)



In [ ]: