MetPy Georeferenced Level 3 Radar Sweep

Question

I am trying to get level 3 data into a georeferenced png Mapbox can use. Mapbox requires a format that im not sure how to implement properly. They use top-right, top-left, bottom-left, and bottom-right coordinates to put a raster image on the map. Such as

"coordinates": [ [-80.425, 46.437], [-71.516, 46.437], [-71.516, 37.936], [-80.425, 37.936] ]

I would prefer to write the information to my metadata file as it already needs to be loaded my the application. Is anyone able to point me in the right direction to form this data in order to get my images rendered properly?

Here is the code I have currently written to get the PNG built, now I just need to georeference it.

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import json
import sys
import os

#from metpy.cbook import get_test_data
from metpy.io import Level3File
from metpy.plots import add_metpy_logo, add_timestamp, ctables
from datetime import datetime

###########################################
#fig, axes = plt.subplots(1, 2, figsize=(15, 8))

radar = 'KCLE'

class metaData:
    def __init__(self, lat,lon,updated):
        self.lat = lat
        self.lon = lon
        self.updated = updated
        
    def toJSON(self):
        return json.dumps(self, default=lambda o: o.__dict__, 
            sort_keys=True, indent=4)



#SYSTEM ARG

radar = sys.argv[1]
product = sys.argv[2]

#END OF SYSTEM ARG

try:
    with open('/mnt/nexrad/' + radar + '/' + product + '/metadata.json', 'r') as file:
        jsonFile = file.read().replace('\n', '')
except IOError:
    dataFile = metaData(0, 0, str(datetime.utcnow()) + 'Z')
    os.makedirs(os.path.dirname('/mnt/nexrad/' + radar + '/' + product + '/'), exist_ok=True)
    with open('/mnt/nexrad/' + radar + '/' + product + '/metadata.json', "w") as outfile:
        outfile.write(dataFile.toJSON())

print(jsonFile)

metaDataObject = json.loads(jsonFile)
    


f = Level3File('/mnt/nexrad/'+radar+'/'+product+'/raw')


dataFile = metaData(f.lat, f.lon, str(f.metadata['prod_time'].utcnow()) + 'Z')

print(dataFile.toJSON())
print(datetime.strptime(metaDataObject['updated'],'%Y-%m-%d %H:%M:%S.%fZ'))

latestUpdate = datetime.strptime(metaDataObject['updated'],'%Y-%m-%d %H:%M:%S.%fZ')
rawUpdateTime = f.metadata['prod_time']

print('Latest Update: ',latestUpdate)
print('File Update  : ',rawUpdateTime)

if rawUpdateTime > latestUpdate:
    print('Updating ' + radar + '...')
    
    fig=plt.figure(figsize=(100,100), dpi=100)
    ax=plt.subplot(1,1,1)
    ax.axis('off')


    datadict = f.sym_block[0][0]
#print(datadict)

# Turn into an array using the scale specified by the file
    data = f.map_data(datadict['data'])

#SHOULD BE ADDED
#lon, lat, _ = pyproj.Geod(ellps='WGS84').fwd(ctr_lon, ctr_lat, azimuth, distance)
#x, y = pyproj.Proj(3857)(lon, lat)

# Grab azimuths and calculate a range based on number of gates
    az = np.array(datadict['start_az'] + [datadict['end_az'][-1]])
    rng = np.linspace(0, f.max_range, data.shape[-1] + 1)

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

# Plot the data
#norm, cmap = colortables.get_with_steps(*ctable)
#cmap="BrBG_r"
    cmap = ctables.registry.get_colortable('NWSStormClearReflectivity')
    norm = mpl.colors.Normalize(vmin=-1, vmax=80)
#ax.pcolormesh(xlocs, ylocs, data, norm=norm, cmap=cmap)
#ax.pcolormesh(xlocs, ylocs, data, norm=Normalize(-25, 75), cmap=cmap)
    ax.pcolor(xlocs, ylocs, data, cmap=cmap, norm=norm)


#ax.set_aspect('auto')

#ax.set_xlim(-320, 320)
#ax.set_ylim(-320, 320)

#add_timestamp(ax, f.metadata['prod_time'], y=0.02, high_contrast=True)
    fig.savefig('/mnt/nexrad/'+radar+'/'+product+'/NOQ.png', transparent=True) #,bbox_inches='tight'
    with open('/mnt/nexrad/'+radar+'/'+product + '/metadata.json', "w") as outfile:
        outfile.write(dataFile.toJSON())
    
    plt.show()

Thanks!

Answer 1

I think the solution here is to carefully set your limits with a maximum range (like the 320km you have commented out). From there, you know that the maximum range at the corner (from pythagorean theorem) is max_range * math.sqrt(2). Once you have the range at the corners, you can use metpy.calc.azimuth_range_to_lat_lon to turn that into lat/lons:

import metpy.calc as mpcalc
from metpy.units import units
import numpy as np

corner_az = np.array([45., 135., 225., 315.]) * units.degrees
corner_range = 320. * np.sqrt(2) * units.km
ctr_lat = 35.25
ctr_lon = -97.1
mpcalc.azimuth_range_to_lat_lon(corner_az, corner_range, ctr_lon, ctr_lat)

where ctr_lon and ctr_lat are the radar longitude and latitude, respectively.