Dimensionality of MetPy NEXRAD level II Fields do not Match. Uncertain how to Conduct A Spatial Analysis.

[msw17002]

What can be better?

I'm trying to pair NEXRAD coordinates into a pandas (then geopandas) dataframe.

To do this, I,

1. Read a NEXRAD level II file using MetPy: radar = Level2File(inn)
2. Calculate the lat/lon coordinates using MetPy (at sweep 0): xlocs, ylocs = azimuth_range_to_lat_lon(az, ref_range, cent_lon, cent_lat)
3. Determine base reflectivity (sweep = 0): The calculation is long, but base reflectivity == data
4. Then flatten xlocs, ylocs, and data into a pandas dataframe: pairs = pd.DataFrame({'Lat': ylocs.flatten(), 'Lon': xlocs.flatten(), 'reflectivity': data.flatten()})

The problem is, data.flatten() does not have the same dimensionality as ylocs.flatten()/xlocs.flatten(). For my case, the non-flattened dimensions are; ylocs = (721, 1833), xlocs = (721, 1833), and data = (720, 1832).

Bottom line, is there a way to pair ylocs/xlocs onto data's grid so that I can index what the coordinates are at cell data(x,y)? Since it was able to plot in pcolormesh, I'm sure there's a way. I just need to find it.

Thank you!

[dopplershift]

Can you share the code that gives you az and ref_range?

If this came from our examples, those arrays are intentionally bigger than the data such that you are working with the coordinates for the corners of the quadrilateral that bounds the area that has a value of data[i, j]; this is what's truly needed to feed to pcolormesh.

[msw17002]

I added my sample script to this email... I did copy it from one of the examples (https://unidata.github.io/MetPy/latest/examples/formats/NEXRAD_Level_2_File.html#sphx-glr-examples-formats-nexrad-level-2-file-py). If the radar file doesn't upload to this reply, you can download it via this link: https://noaa-nexrad-level2.s3.amazonaws.com/2021/12/21/KTBW/KTBW20211221_104123_V06.

This is the important part you're asking for,

#===radar file
radar = Level2File(inn)
#---retrieve base fields 
cent_lon = radar.sweeps[0][0][1].lon
cent_lat = radar.sweeps[0][0][1].lat
#-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 radar.sweeps[sweep]])
diff = np.diff(az)
crossed = diff < -180
diff[crossed] += 360.
avg_spacing = diff.mean()
#-Convert mid-point to edge
az = (az[:-1] + az[1:]) / 2
az[crossed] += 180.
#-Concatenate with overall start and end of data we calculate using the average spacing
az = np.concatenate(([az[0] - avg_spacing], az, [az[-1] + avg_spacing]))
az = units.Quantity(az, 'degrees')
#-5th item is a dict mapping a var name (byte string) to a tuple
# of (header, data array)
ref_hdr = radar.sweeps[sweep][0][4][b'REF'][0]
ref_range = (np.arange(ref_hdr.num_gates + 1) - 0.5) * ref_hdr.gate_width + ref_hdr.first_gate
ref_range = units.Quantity(ref_range, 'kilometers')
ref = np.array([ray[4][b'REF'][1] for ray in radar.sweeps[sweep]])
#-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, ylocs = azimuth_range_to_lat_lon(az, ref_range, cent_lon, cent_lat)

I later attempt to determine the base reflectivity at a specific location given;

#-determine smallest distance
row = np.nanargmin(((xlocs.flatten()-qlon)**2+(ylocs.flatten()-qlat)**2)**0.5)
val = data.flatten()[row]
if np.ma.isMaskedArray(val): val = str(np.nan)
else:                        val = str(int(val))

I'd be satisfied with the above if the dimensions of xlocs==ylocs==data.

example.txt

[dopplershift]

To accomplish what you want would require changing the range calculation, which is intentionally adding 1 to the size to:

ref_range = np.arange(ref_hdr.num_gates) * ref_hdr.gate_width + ref_hdr.first_gate
ref_range = units.Quantity(ref_range, 'kilometers')

Then also avoid modifying the collection of azimuths, but only have this line:

az = np.array([ray[0].az_angle for ray in radar.sweeps[sweep]])

eliminate all the other lines that affect az.