MetPy地理参考3级雷达扫描

Question

我正在尝试将3级数据获取到地理参考的png Mapbox中。Mapbox需要一种我不确定如何正确实现的格式。它们使用右上、左上、左下和右下坐标在地图上放置光栅图像。如

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

我更愿意将信息写入我的元数据文件，因为它已经需要加载到我的应用程序中。有人能为我指明正确的方向来形成这些数据，以便正确渲染我的图像吗？

这是我目前为构建PNG而编写的代码，现在我只需要对其进行地理参考

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()

谢谢！

Answer 1

我认为这里的解决方案是小心地设置最大续航里程的限制(比如你评论的320公里(。从那里，你知道拐角处的最大范围(根据勾股定理(是max_range * math.sqrt(2)。一旦你在拐角处有了射程，你就可以使用metpy.calc.azimuth_range_to_lat_lon将其转换为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)

其中ctr_lon和ctr_lat分别是雷达经度和纬度。