Revisit MetPy for GOES, NEXRAD stuff

Right now MetPy is at version 0.10.0 and it's a little too rough around the edges for me. They have support for some things that I really want, like the NWSReflectivity colortable and the baked-in US state & county lines that things like Cartopy and py-art don't have.

In the future, when MetPy gets closer to a 1.0-like release, it'd be good to go back and revisit things. For the NEXRAD data in particular right now, I'm not enthused enough to switch away from py-art. Consider:

radar = pyart.io.read_nexrad_archive(filename)

# Pull out the identifiers
site = radar.metadata['instrument_name']
siteLat = radar.latitude['data'][0]
siteLon = radar.longitude['data'][0]
dprod = radar.metadata['original_container']

display = RadarMapDisplay(qcradar)
fig = plt.figure(figsize=figsize, dpi=100)
display.plot_ppi_map('reflectivity_masked',
                                   mask_outside=True,
                                   vmin=-20, vmax=75,
                                   min_lon=lonMin, max_lon=lonMax,
                                   min_lat=latMin, max_lat=latMax,
                                   projection=crs,
                                   fig=fig,
                                   cmap=cmap,
                                   lat_0=siteLat,
                                   lon_0=siteLon,
                                   embelish=False,
                                   colorbar_flag=True,
                                   title_flag=False,
                                   ticklabs=[],
                                   ticks=[],
                                   lat_lines=[],
                                   lon_lines=[],
                                   raster=True)

display.plot_point(siteLon, siteLat, symbol='^', color='orange')
plt.show()
plt.close()

to this mess in MetPy that I'm struggling thru:

dat = metpy.io.Level2File(filename)

# Most of this is from the NEXRAD Level 2 example in the MetPy docs
# https://unidata.github.io/MetPy/latest/examples/formats/NEXRAD_Level_2_File.html

# Pull data out of the file
sweep = 0

# First item in ray is header, which has azimuth angle
az = np.array([ray[0].az_angle for ray in dat.sweeps[sweep]])

# 5th item is a dict mapping a var name (byte string) to a tuple
# of (header, data array). "REF" is the base reflectivity
refHdr = dat.sweeps[sweep][0][4][b'REF'][0]
refRange = np.arange(refHdr.num_gates) * refHdr.gate_width + refHdr.first_gate
ref = np.array([ray[4][b'REF'][1] for ray in dat.sweeps[sweep]])

# RHO is the cross correlation ratio
rhoHdr = dat.sweeps[sweep][0][4][b'RHO'][0]
rhoRange = (np.arange(rhoHdr.num_gates + 1) - 0.5) * rhoHdr.gate_width + rhoHdr.first_gate
rho = np.array([ray[4][b'RHO'][1] for ray in dat.sweeps[sweep]])

# Plot the data (just the reflectivity for now)
fig, ax = plt.subplots(1, 1, figsize=(7, 7))

# Turn into an array, then mask
data = np.ma.array(ref)
data[np.isnan(data)] = np.ma.masked

# Convert az,range to x,y
xlocs = refRange * np.sin(np.deg2rad(az[:, np.newaxis]))
ylocs = refRange * np.cos(np.deg2rad(az[:, np.newaxis]))

cmap = colortables.get_colortable('NWSReflectivity')
ax.pcolormesh(xlocs, ylocs, data, cmap=cmap, vmin=0, vmax=75)
ax.set_aspect('equal', 'datalim')
ax.set_xlim(-40, 20)
ax.set_ylim(-30, 30)

plt.show()
plt.close()

the above one isn't even projected, yet, it's still in raw data units (distance from radar in km). Once MetPy reaches a more stable/streamlined interface it's really worth taking another look.

LowellObservatory / NightShift

Revisit MetPy for GOES, NEXRAD stuff #5