Shown here are plots for Base Reflectivity (N0Q, 94) and Base Velocity (N0U, 99) using AWIPS data rendered with Matplotlib, Cartopy, and MetPy. This example improves upon existing Level 3 Python rendering by doing the following:
Display scaled and labeled colorbar below each figure.
Plot radar radial images as coordinate maps in Cartopy and label with lat/lon.
8 bit Z and V colormap and data scaling added to MetPy from operational AWIPS.
Level 3 data are retrieved from the Unidata EDEX Cloud server (
edex-cloud.unidata.ucar.edu)Raw HDF5 byte data are converted to product values and scaled according to (page 3-34 https://www.roc.noaa.gov/wsr88d/PublicDocs/ICDS/2620001U.pdf)
The threshold level fields are used to describe (up to) 256 levels as follows: halfword 31 contains the minimum data value in m/s*10 (or dBZ*10) halfword 32 contains the increment in m/s*10 (or dBZ*10) halfword 33 contains the number of levels (0 - 255)
According to the ICD for the Product Specification, “the 256 data levels of the digital product cover a range of reflectivity between -32.0 to +94.5 dBZ, in increments of 0.5 dBZ. Level codes 0 and 1 correspond to ‘Below Threshold’ and ‘Range Folded’, respectively, while level codes 2 through 255 correspond to the reflectivity data itself”.
So it’s really 254 color values between -32 and +94.5 dBZ.
The ICD lists 16 specific color levels and directs 256-level reflectivity products to use corresponding colors, leaving it the rendering application to scale and blend between the 16 color values, and to make decisions about discrete color changes, apparently. 
For AWIPS, the National Weather Service uses a mostly-blended color scale with a discrete jump to red at reflectivity values of 50 dBZ:
50 dBZ corresponds to the 16-level color light red (FF6060). Note that FF6060 is not used in the NWS AWIPS color scale, instead RGB value is given as 255,0,0 (hex code FF0000). 60 dBZ is not quite exactly where white starts, but it makes sense that it would. Obviously the AWIPS D2D authors took some liberties with their 256-level rendering, not adhering strictly to “dark red” for dBZ values between 60-65 (white was for 70 dBZ and above on the 16-level colormap). For this exercise we will assume 50 dBZ should be red and 60 dBZ white, and 75 dBZ cyan.
Setup
pip install python-awips matplotlib cartopy metpy
Python Script
Download this script as a Jupyter Notebook.
from awips.dataaccess import DataAccessLayer from awips import ThriftClient, RadarCommon from dynamicserialize.dstypes.com.raytheon.uf.common.time import TimeRange from dynamicserialize.dstypes.com.raytheon.uf.common.dataplugin.radar.request import GetRadarDataRecordRequest from datetime import datetime from datetime import timedelta import matplotlib.pyplot as plt import numpy as np from numpy import ma from metpy.plots import ctables import cartopy.crs as ccrs from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER # set EDEX server and radar site definitions site = 'kmux' DataAccessLayer.changeEDEXHost('edex-cloud.unidata.ucar.edu') request = DataAccessLayer.newDataRequest() request.setDatatype('radar') request.setLocationNames(site) # Get latest time for site datatimes = DataAccessLayer.getAvailableTimes(request) dateTimeStr = str(datatimes[-1]) dateTimeStr = "2017-02-02 03:53:03" buffer = 60 # seconds dateTime = datetime.strptime(dateTimeStr, '%Y-%m-%d %H:%M:%S') # Build timerange +/- buffer beginRange = dateTime - timedelta(0, buffer) endRange = dateTime + timedelta(0, buffer) timerange = TimeRange(beginRange, endRange) # GetRadarDataRecordRequest to query site with timerange client = ThriftClient.ThriftClient('edex-cloud.unidata.ucar.edu') request = GetRadarDataRecordRequest() request.setTimeRange(timerange) request.setRadarId(site) # Map config def make_map(bbox, projection=ccrs.PlateCarree()): fig, ax = plt.subplots(figsize=(12, 12), subplot_kw=dict(projection=projection)) ax.set_extent(bbox) ax.coastlines(resolution='50m') gl = ax.gridlines(draw_labels=True) gl.xlabels_top = gl.ylabels_right = False gl.xformatter = LONGITUDE_FORMATTER gl.yformatter = LATITUDE_FORMATTER return fig, ax # ctable defines the colortable, beginning value, data increment # * For N0Q the scale is -20 to +75 dBZ in increments of 0.5 dBZ # * For N0U the scale is -100 to +100 kts in increments of 1 kt nexrad = {} nexrad["N0Q"] = { 'id': 94, 'unit':'dBZ', 'name':'0.5 deg Base Reflectivity', 'ctable': ['NWSStormClearReflectivity',-20., 0.5], 'res': 1000., 'elev': '0.5' } nexrad["N0U"] = { 'id': 99, 'unit':'kts', 'name':'0.5 deg Base Velocity', 'ctable': ['NWS8bitVel',-100.,1.], 'res': 250., 'elev': '0.5' } grids = [] for code in nexrad: request.setProductCode(nexrad[code]['id']) request.setPrimaryElevationAngle(nexrad[code]['elev']) response = client.sendRequest(request) if response.getData(): for record in response.getData(): # Get record hdf5 data idra = record.getHdf5Data() rdat,azdat,depVals,threshVals = RadarCommon.get_hdf5_data(idra) dim = rdat.getDimension() lat,lon = float(record.getLatitude()),float(record.getLongitude()) radials,rangeGates = rdat.getSizes() # Convert raw byte to pixel value rawValue=np.array(rdat.getByteData()) array = [] for rec in rawValue: if rec<0: rec+=256 array.append(rec) if azdat: azVals = azdat.getFloatData() az = np.array(RadarCommon.encode_radial(azVals)) dattyp = RadarCommon.get_data_type(azdat) az = np.append(az,az[-1]) header = RadarCommon.get_header(record, format, rangeGates, radials, azdat, 'description') rng = np.linspace(0, rangeGates, rangeGates + 1) # Convert az/range to a lat/lon from pyproj import Geod g = Geod(ellps='clrk66') center_lat = np.ones([len(az),len(rng)])*lat center_lon = np.ones([len(az),len(rng)])*lon az2D = np.ones_like(center_lat)*az[:,None] rng2D = np.ones_like(center_lat)*np.transpose(rng[:,None])*nexrad[code]['res'] lons,lats,back=g.fwd(center_lon,center_lat,az2D,rng2D) bbox = [lons.min(), lons.max(), lats.min(), lats.max()] # Create 2d array multiArray = np.reshape(array, (-1, rangeGates)) data = ma.array(multiArray) # threshVals[0:2] contains halfwords 31,32,33 (min value, increment, num levels) data = ma.array(threshVals[0]/10. + (multiArray)*threshVals[1]/10.) if nexrad[code]['unit'] == 'kts': data[data<-63] = ma.masked data *= 1.94384 # Convert to knots else: data[data<=((threshVals[0]/10.)+threshVals[1]/10.)] = ma.masked # Save our requested grids so we can render them multiple times product = { "code": code, "bbox": bbox, "lats": lats, "lons": lons, "data": data } grids.append(product) print("Processed "+str(len(grids))+" grids.") Processed 2 grids. Plot N0Q and N0U with Cartopy
for rec in grids: code = rec["code"] bbox = rec["bbox"] lats = rec["lats"] lons = rec["lons"] data = rec["data"] # Create figure %matplotlib inline fig, ax = make_map(bbox=bbox) # Colortable filename, beginning value, increment ctable = nexrad[code]['ctable'][0] beg = nexrad[code]['ctable'][1] inc = nexrad[code]['ctable'][2] norm, cmap = ctables.registry.get_with_steps(ctable, beg, inc) cs = ax.pcolormesh(lons, lats, data, norm=norm, cmap=cmap) ax.set_aspect('equal', 'datalim') cbar = plt.colorbar(cs, extend='both', shrink=0.75, orientation='horizontal') cbar.set_label(site.upper()+" "+ str(nexrad[code]['res']/1000.) +"km " +nexrad[code]['name']+" ("+code+") " +nexrad[code]['unit']+" " +str(record.getDataTime())) # Zoom to within +-2 deg of center ax.set_xlim(lon-2., lon+2.) ax.set_ylim(lat-2., lat+2.) plt.show() 
compare with the same product scan rendered in AWIPS CAVE (slightly different projections and still some color mapping differences, most noticeable in ground clutter).
Two-panel plot, zoomed in
fig, axes = plt.subplots(ncols=2,figsize=(12,9), subplot_kw=dict(projection=ccrs.PlateCarree())) i=0 for rec,ax in zip(grids, axes): code = rec["code"] bbox = rec["bbox"] lats = rec["lats"] lons = rec["lons"] data = rec["data"] # Create figure ax.set_extent(bbox) ax.coastlines(resolution='50m') gl = ax.gridlines(draw_labels=True) gl.xlabels_top = gl.ylabels_right = False if i>0: gl.ylabels_left = False # hide right-pane left axis label gl.xformatter = LONGITUDE_FORMATTER gl.yformatter = LATITUDE_FORMATTER # Colortable filename, beginning value, increment colorvals=nexrad[code]['ctable'] ctable = nexrad[code]['ctable'][0] beg = nexrad[code]['ctable'][1] inc = nexrad[code]['ctable'][2] norm, cmap = ctables.registry.get_with_steps(ctable, beg, inc) cs = ax.pcolormesh(lons, lats, data, norm=norm, cmap=cmap) ax.set_aspect('equal', 'datalim') cbar = fig.colorbar(cs, orientation='horizontal', ax=ax) cbar.set_label(site.upper()+" "+code+" "+nexrad[code]['unit']+" "+str(record.getDataTime())) plt.tight_layout() # Zoom ax.set_xlim(lon-.1, lon+.1) ax.set_ylim(lat-.1, lat+.1) i+=1 
and again compared to CAVE
I'd love to see more great tutorials like this. Perhaps an "Intro to Python-AWIPS and the Jupyter Notebook Examples" tutorial that goes beyond the python-awips installation steps noted in the documentation. I'm particularly interested in the initial setting up/running of the Jupyter Notebook examples. This was quickly covered in the on-site AWIPS training I attended in Oct 2016, but I'm failing to remember the correct steps for setting it up (including which steps are needed, if any, on an EDEX server and/or CAVE client). Thanks!
Posted by Jason Kaiser on July 07, 2017 at 08:13 AM MDT #