Initial commit

plot_nexrad_level3-2panels.py

+#!/usr/bin/python
+import matplotlib as mpl
+mpl.use('Agg')
+import matplotlib.pyplot as plt
+import numpy as np
+import sys
+import time
+import datetime
+
+import cartopy
+import cartopy.crs as ccrs
+import cartopy.feature as cfeature
+import warnings
+
+from metpy.cbook import get_test_data
+from metpy.io import Level3File
+from metpy.plots import add_timestamp, colortables
+from metpy.plots import USCOUNTIES
+
+# SOURCE1: https://unidata.github.io/MetPy/latest/examples/formats/NEXRAD_Level_2_File.html
+# SOURCE2: http://nbviewer.jupyter.org/gist/dopplershift/356f2e14832e9b676207
+# SOURCE3: https://unidata.github.io/MetPy/latest/examples/formats/NEXRAD_Level_3_File.html
+
+#################################################
+main_directory='/store/sfcnet/datasets/nexrad/projects/KSGF/level3'
+in_directory=main_directory+'/source'
+out_directory=main_directory+'/results'
+
+dpi=100
+city_lat=36.6437
+city_lon=-93.2185
+input_box=[-96, -90, 35, 39] # KSGF
+
+#################################################
+# BEGIN PROGRAM
+start=time.time()
+#################################################
+
+# Read in Arguments 
+if len(sys.argv) < 3:
+  sys.exit("USAGE: <INFILE1> <INFILE2>")
+
+infile1=in_directory+'/'+str(sys.argv[1])
+infile2=in_directory+'/'+str(sys.argv[2])
+
+plt.figure(figsize=(15, 8))
+plt.suptitle('Springfield, MO Radar (KSGF) For July 19th, 2018',size=20)
+
+figcounter=1
+for filName, ctable in zip((infile1, infile2), ('NWSReflectivity', 'NWSVelocity')):
+    # Open the file
+    name = filName
+    print(name)
+    f = Level3File(name)
+
+    radar_lat=f.lat
+    radar_lon=f.lon
+
+    #print("CITY LAT/LON",city_lat,city_lon)
+    #print("RADAR LAT/LON",radar_lat,radar_lon)
+
+    # Pull the data out of the file object
+    datadict = f.sym_block[0][0]
+
+    # Turn into an array, then mask
+    data = np.ma.array(datadict['data'])
+    data[data == 0] = np.ma.masked
+
+    # 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)
+    rng=rng*1000
+
+    #####################
+    # PLOT
+    print("PLOT")
+
+    # Plot the data
+    norm, cmap = colortables.get_with_steps(ctable, 16, 16)
+
+    # Open Figure (Centered on Radar Lat/Lon)
+    proj = cartopy.crs.LambertConformal(central_longitude=radar_lon, central_latitude=radar_lat)
+    ax = plt.subplot(1, 2, figcounter, projection=proj)        
+
+    # Convert Radar Lat/Lon Point to LCC projection
+    radar_x, radar_y = proj.transform_points(ccrs.Geodetic(), np.array([radar_lon]), np.array([radar_lat]))[0,:-1] 
+
+    # Convert az,range to x,y....using Radar Lat/Lon Projection as (0,0)
+    xlocs = radar_x + rng * np.sin(np.deg2rad(az))[:,None]
+    ylocs = radar_y + rng * np.cos(np.deg2rad(az))[:,None]
+
+    # Convert LCC projection x/y points to lat/lon points
+    lons, lats = ccrs.Geodetic().transform_points(proj,xlocs,ylocs)[:,:,:-1].T
+    lons=lons.T
+    lats=lats.T
+
+    #print(np.min(lons),np.max(lons))
+    #print(np.min(lats),np.max(lats))
+
+    # Add Boundaries
+    ax.coastlines('50m', 'black', linewidth=2, zorder=2)
+    state_borders = cartopy.feature.NaturalEarthFeature(
+                    category='cultural', name='admin_1_states_provinces_lines',
+                    scale='10m', facecolor='none')
+    ax.add_feature(state_borders, edgecolor='black', linewidth=1, zorder=3)
+    ax.add_feature(USCOUNTIES.with_scale('500k'),linewidth=0.25,zorder=3)
+
+    # Set limits in lat/lon space
+    ax.set_extent(input_box)
+
+    # Add ocean and land background
+    ocean = cartopy.feature.NaturalEarthFeature('physical', 'ocean', scale='50m',edgecolor='face',facecolor=cartopy.feature.COLORS['water'])
+    land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='50m',edgecolor='face',facecolor=cartopy.feature.COLORS['land'])
+    ax.add_feature(ocean, zorder=-1)
+    ax.add_feature(land, zorder=-1)
+    ax.add_feature(cfeature.BORDERS,linewidth=0.5)
+
+    # Plot the Radar Data 
+    cs = ax.pcolormesh(lons, lats, data,transform=ccrs.PlateCarree(),cmap=cmap,norm=norm,zorder=0) #Map, no data
+    ax.set_aspect('equal', 'datalim')
+
+    # Plot Individual Points
+    plt.plot(city_lon,city_lat, transform=ccrs.PlateCarree(), marker='*',markerfacecolor='red',markeredgecolor='black',markersize=20,label="Table Rock Lake")
+    plt.scatter(radar_lon, radar_lat, transform=ccrs.PlateCarree(), marker='*',c='black',s=100,zorder=0,label="Radar Location")
+
+    # Plot Legend (isn't working for some reason)
+    #plt.legend(loc=0)
+    #legend = ax.legend(loc='best', shadow=True, fontsize='x-large')
+    #legend.get_frame().set_facecolor('C0')
+
+    # Plot Time
+    current_time=sys.argv[figcounter]
+    #print(current_time)
+    current_date=datetime.datetime(int(current_time[19:23]),int(current_time[23:25]),int(current_time[25:27]),int(current_time[27:29]),int(current_time[29:31])) - datetime.timedelta(hours=5)  
+    current_day= "%04i%02i%02i" % (current_date.year,current_date.month,current_date.day)                                   
+    current_hour= "%02i:%02ipm" % (current_date.hour-12,current_date.minute) 
+    plt.annotate('Made By Jared Rennie (@jjrennie) | Time='+str(current_day)+' '+str(current_hour)+' CDT',xy=(0.995, 0.01), xycoords='axes fraction', fontsize=10,backgroundcolor='black',color='white',horizontalalignment='right', verticalalignment='bottom',zorder=4)
+
+    # Plot Color Bar
+    cbar = plt.colorbar(cs, orientation='horizontal', ax=ax, pad=0.01)
+    cbar.set_label(ctable)
+
+    figcounter+=1
+
+plt.savefig(out_directory+'/lv3-2panels-'+current_time[19:31]+'.png', format='png', dpi=dpi,bbox_inches='tight')
+plt.clf()
+####################
+# DONE
+####################
+print("DONE!")
+end=time.time()
+print ("Runtime: %8.1f seconds." % (end-start))
+sys.exit()
\ No newline at end of file