Skip to content

Revisit MetPy for GOES, NEXRAD stuff #5

Description

@astrobokonon

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.

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

Labels

No labels
No labels

Type

No type

Projects

No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions