225 lines
8.0 KiB
Python
225 lines
8.0 KiB
Python
#Lets import some stuff!
|
|
import boto3
|
|
from botocore import UNSIGNED
|
|
from botocore.client import Config
|
|
from datetime import timedelta, datetime, timezone
|
|
import os
|
|
import pyart
|
|
from matplotlib import pyplot as plt
|
|
import tempfile
|
|
import numpy as np
|
|
import io
|
|
from PIL import Image
|
|
import threading
|
|
import concurrent.futures
|
|
|
|
import cartopy.crs as ccrs
|
|
from cartopy.io.img_tiles import OSM
|
|
|
|
def printdebug(msg):
|
|
print(f'DEBUG: {msg}')
|
|
|
|
#Helper function for the search
|
|
def _nearestDate(dates, pivot):
|
|
return min(dates, key=lambda x: abs(x - pivot))
|
|
|
|
def _getRadarFileNames(site, datetime_t):
|
|
return datetimes
|
|
|
|
def get_radar_frames(site, datetime_t, num_frames):
|
|
#First create the query string for the bucket knowing
|
|
#how NOAA and AWS store the data
|
|
|
|
currentDateString = datetime_t.strftime('%Y/%m/%d/') + site
|
|
previousDateString = (datetime_t - timedelta(days=1)).strftime('%Y/%m/%d/') + site
|
|
|
|
conn = boto3.resource('s3', config=Config(signature_version=UNSIGNED))
|
|
bucket = conn.Bucket('noaa-nexrad-level2')
|
|
|
|
bucketListCurrentDate = list(bucket.objects.filter(Prefix = currentDateString))
|
|
bucketListPreviousDate = list(bucket.objects.filter(Prefix = previousDateString))
|
|
|
|
# we are going to create a list of keys and datetimes to allow easy searching
|
|
|
|
keys = []
|
|
datetimes = []
|
|
|
|
#populate the list
|
|
printdebug('Getting nearby dates')
|
|
for i in range(len(bucketListCurrentDate)):
|
|
this_str = str(bucketListCurrentDate[i].key)
|
|
if 'gz' in this_str:
|
|
endme = this_str[-22:-4]
|
|
fmt = '%Y%m%d_%H%M%S_V0'
|
|
dt = datetime.strptime(endme, fmt).replace(tzinfo=timezone.utc)
|
|
datetimes.append(dt)
|
|
keys.append(bucketListCurrentDate[i])
|
|
|
|
if this_str[-3::] == 'V06':
|
|
endme = this_str[-19::]
|
|
fmt = '%Y%m%d_%H%M%S_V06'
|
|
dt = datetime.strptime(endme, fmt).replace(tzinfo=timezone.utc)
|
|
datetimes.append(dt)
|
|
keys.append(bucketListCurrentDate[i])
|
|
|
|
#datetimes = _getRadarFileNames(site, datetime_t)
|
|
printdebug('Getting nearest datetime')
|
|
closest_datetime = _nearestDate(datetimes, datetime_t)
|
|
index = datetimes.index(closest_datetime)
|
|
|
|
#radarFrames = [0] * num_frames
|
|
radarFrames = {}
|
|
i = index
|
|
printdebug('Pulling all the radar images requested')
|
|
while i > (index - num_frames) and i >= 0:
|
|
#radarFrames[num_frames - ((index - i) + 1)] = get_radar_from_aws(keys[i].Object())
|
|
filename = keys[i].key.replace('/', '_')
|
|
radarFrames[filename] = get_radar_from_aws(keys[i].Object())
|
|
i = i - 1
|
|
|
|
sortedIndices = sorted(radarFrames)
|
|
sortedRadarFrames = {}
|
|
for key in sortedIndices:
|
|
sortedRadarFrames[key] = radarFrames[key]
|
|
|
|
return sortedRadarFrames
|
|
|
|
def get_radar_from_aws(aws_object):
|
|
"""
|
|
Get the closest volume of NEXRAD data to a particular datetime.
|
|
Parameters
|
|
----------
|
|
site : string
|
|
four letter radar designation
|
|
datetime_t : datetime
|
|
desired date time
|
|
Returns
|
|
-------
|
|
radar : Py-ART Radar Object
|
|
Radar closest to the queried datetime
|
|
"""
|
|
filename = aws_object.key.replace('/', '_')
|
|
if not os.path.isfile(f'cache/rawradar/{filename}'):
|
|
printdebug('Cache Miss')
|
|
aws_object.download_file(f'cache/rawradar/{filename}')
|
|
#radar = pyart.io.read(f'cache/rawradar/{filename}')
|
|
radar = f'cache/rawradar/{filename}'
|
|
return radar
|
|
|
|
def convert_fig_to_img(fig):
|
|
fig.canvas.draw()
|
|
(w,h) = fig.canvas.get_width_height()
|
|
img = Image.frombytes('RGBA', fig.canvas.get_width_height(), np.frombuffer(fig.canvas.buffer_rgba(), dtype=np.uint8).reshape((h,w,4)))
|
|
return img
|
|
|
|
def generate_radar_image(radar_data, tiler, mercator):
|
|
#determine max_lat / max_lon etc from position of radar
|
|
max_lat = radar_data.latitude["data"][0] + 4.5
|
|
min_lat = radar_data.latitude["data"][0] - 4.5
|
|
max_lon = radar_data.longitude["data"][0] + 5.5
|
|
min_lon = radar_data.longitude["data"][0] - 5.5
|
|
|
|
# Used to be .5
|
|
lal = np.arange(min_lat, max_lat, .2)
|
|
lol = np.arange(min_lon, max_lon, .1)
|
|
|
|
printdebug('Creating RadarMapDisplay')
|
|
display = pyart.graph.RadarMapDisplay(radar_data)
|
|
|
|
projection = ccrs.LambertConformal(
|
|
central_latitude=radar_data.latitude["data"][0],
|
|
central_longitude=radar_data.longitude["data"][0],
|
|
)
|
|
|
|
printdebug('Processing data to clean it up a bit')
|
|
|
|
radar_data.fields['reflectivity']['data'][:, -10:] = np.ma.masked
|
|
|
|
gatefilter = pyart.filters.GateFilter(radar_data)
|
|
gatefilter.exclude_transition()
|
|
gatefilter.exclude_masked("reflectivity")
|
|
gatefilter.exclude_below('reflectivity', 5)
|
|
#gatefilter = pyart.filters.moment_based_gate_filter(radar_data)
|
|
|
|
printdebug('Despeckling')
|
|
despeckle_gatefilter = pyart.correct.despeckle_field(radar_data, 'reflectivity', gatefilter=gatefilter)
|
|
|
|
fig = plt.figure(figsize = [10,8])
|
|
fig.patch.set_facecolor('none')
|
|
ax = fig.add_subplot(1, 1, 1, projection=mercator)
|
|
#ax = fig.add_subplot(1, 1, 1)
|
|
ax.set_extent([min_lon, max_lon, min_lat, max_lat], crs=ccrs.PlateCarree())
|
|
ax.add_image(tiler, 6)
|
|
printdebug('Plotting!')
|
|
display.plot_ppi_map('reflectivity', sweep = 0, resolution='10m',
|
|
vmin = -8, vmax = 64,
|
|
#cmap = "ChaseSpectral",
|
|
cmap = "NWSRef",
|
|
min_lat = min_lat, min_lon = min_lon,
|
|
max_lat = max_lat, max_lon = max_lon,
|
|
lat_lines = lal, lon_lines = lol,
|
|
projection=projection, fig=fig, ax=ax, title_flag=False,
|
|
lat_0=radar_data.latitude["data"][0], lon_0=radar_data.longitude["data"][0],
|
|
gatefilter=gatefilter, add_grid_lines=False, alpha=0.8)
|
|
|
|
display.plot_range_ring(radar_data.range["data"][-1] / 1000.0)
|
|
|
|
return fig
|
|
|
|
def _thread_generate_radar_image(radar_data, dataName, tiler, mercator, index, lock):
|
|
if not os.path.isfile(f'cache/radarimages/{dataName}.png'):
|
|
radarData = pyart.io.read(radar_data)
|
|
fig = generate_radar_image(radarData, tiler, mercator)
|
|
frame = convert_fig_to_img(fig)
|
|
frame.save(f'cache/radarimages/{dataName}.png')
|
|
del radarData
|
|
lock.acquire()
|
|
secondary_frames[index] = frame
|
|
#plt.clear()
|
|
lock.release()
|
|
else:
|
|
frame = Image.open(f'cache/radarimages/{dataName}.png')
|
|
lock.acquire()
|
|
secondary_frames[index] = frame
|
|
lock.release()
|
|
|
|
fmt = '%Y%m%d_%H%M%S'
|
|
currentDatetime = datetime.now(timezone.utc)
|
|
num_frames = 20
|
|
secondary_frames = [0] * (num_frames - 1)
|
|
|
|
radar = get_radar_frames('KVNX', currentDatetime, num_frames=num_frames)
|
|
print(list(radar.keys()))
|
|
|
|
tiler = OSM(cache=True)
|
|
mercator = tiler.crs
|
|
|
|
radarKeys = list(radar.keys())
|
|
|
|
if not os.path.isfile(f'cache/radarimages/{radarKeys[0]}.png'):
|
|
radarData = pyart.io.read(radar[radarKeys[0]])
|
|
firstFrameFig = generate_radar_image(radarData, tiler, mercator)
|
|
firstFrame = convert_fig_to_img(firstFrameFig)
|
|
firstFrame.save(f'cache/radarimages/{radarKeys[0]}.png')
|
|
del radarData
|
|
else:
|
|
firstFrame = Image.open(f'cache/radarimages/{radarKeys[0]}.png')
|
|
|
|
i = 0
|
|
lock = threading.Lock()
|
|
threadPool = concurrent.futures.ThreadPoolExecutor(max_workers=4)
|
|
for i in range(num_frames - 1):
|
|
# if not os.path.isfile(f'cache/radarimages/{radarKeys[i + 1]}.png'):
|
|
# fig = generate_radar_image(radar[radarKeys[i + 1]], tiler, mercator)
|
|
# secondary_frames[i] = convert_fig_to_img(fig)
|
|
# secondary_frames[i].save(f'cache/radarimages/{radarKeys[i + 1]}.png')
|
|
# else:
|
|
# secondary_frames[i] = Image.open(f'cache/radarimages/{radarKeys[i + 1]}.png')
|
|
threadPool.submit(_thread_generate_radar_image, radar[radarKeys[i + 1]], radarKeys[i + 1], tiler, mercator, i, lock)
|
|
|
|
threadPool.shutdown(wait=True)
|
|
|
|
datetimeString = datetime.strftime(currentDatetime, fmt)
|
|
|
|
firstFrame.save(f'testradar_{datetimeString}.webp', append_images=secondary_frames, duration=200)
|