Initial commit

testfile cartopy.py

@@ -0,0 +1,186 @@
+#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 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
+    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())
+        i = i - 1
+    
+    return radarFrames
+
+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
+    """
+        
+    #find the closest available radar to your datetime
+
+    #closest_datetime = _nearestDate(datetimes, datetime_t)
+    #index = datetimes.index(closest_datetime)
+
+    #localfile = tempfile.NamedTemporaryFile()
+    #keys[index].Object().download_file(localfile.name)
+    #aws_object.download_file(localfile.name)
+    #radar = pyart.io.read(localfile.name)
+    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}')
+    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', 16)
+
+    printdebug('Despeckling')
+    despeckle_gatefilter = pyart.correct.despeckle_field(radar_data, 'reflectivity', gatefilter=gatefilter, size=20) 
+
+    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 = 10, vmax = 64,
+                        cmap = "ChaseSpectral",
+                        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=despeckle_gatefilter, add_grid_lines=False, alpha=0.8)
+
+    display.plot_range_ring(radar_data.range["data"][-1] / 1000.0)
+
+    return fig
+
+fmt = '%Y%m%d_%H%M%S' 
+currentDatetime = datetime.now(timezone.utc)
+num_frames = 10
+secondary_frames = [0] * (num_frames - 1)
+
+radar = get_radar_frames('KVNX', currentDatetime, num_frames=num_frames)
+
+tiler = OSM(cache=True)
+mercator = tiler.crs
+
+firstFrameFig = generate_radar_image(radar[0], tiler, mercator)
+firstFrame = convert_fig_to_img(firstFrameFig)
+i = 0
+for i in range(num_frames - 1):
+    fig = generate_radar_image(radar[i + 1], tiler, mercator)
+    secondary_frames[i] = convert_fig_to_img(fig)
+    # filename = aws_object.key.replace('/', '_')
+
+datetimeString = datetime.strftime(currentDatetime, fmt)
+
+firstFrame.save(f'testradar_{datetimeString}.gif', append_images=secondary_frames, duration=200)