In [ ]:
#First out imports
#Lets import some stuff!
from datetime import datetime, timedelta
import os
import tempfile
import cartopy.crs as ccrs
from boto.s3.connection import S3Connection
import cartopy
import matplotlib.patheffects as mpatheffects
import matplotlib.pyplot as plt
from netCDF4 import num2date
import numpy as np
import pyart
import pytz
import xarray
import netCDF4
import cartopy.io.img_tiles as cimgt
import pandas
#from botocore.exceptions import BotoServerError
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
%matplotlib inline
In [ ]:
#Now our nifty fetch script
#Helper function for the search
def _nearestDate(dates, pivot):
return min(dates, key=lambda x: abs(x - pivot))
def get_radar_from_aws(site, datetime_t):
"""
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
"""
# First create the query string for the bucket knowing
# how NOAA and AWS store the data
my_pref = datetime_t.strftime('%Y/%m/%d/') + site
# Connect to the bucket
conn = S3Connection(anon = True)
bucket = conn.get_bucket('noaa-nexrad-level2')
# Get a list of files
bucket_list = list(bucket.list(prefix = my_pref))
# we are going to create a list of keys and datetimes to allow easy searching
keys = []
datetimes = []
# populate the list
for i in range(len(bucket_list)):
this_str = str(bucket_list[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)
datetimes.append(dt)
keys.append(bucket_list[i])
if this_str[-3::] == 'V06':
endme = this_str[-19::]
fmt = '%Y%m%d_%H%M%S_V06'
dt = datetime.strptime(endme, fmt)
datetimes.append(dt)
keys.append(bucket_list[i])
# find the closest available radar to your datetime
closest_datetime = _nearestDate(datetimes, datetime_t)
index = datetimes.index(closest_datetime)
localfile = tempfile.NamedTemporaryFile()
keys[index].get_contents_to_filename(localfile.name)
radar = pyart.io.read(localfile.name)
return radar
In [ ]:
In [ ]:
def x_array_from_aws(site, datedesire):
radar = get_radar_from_aws(site, datedesire)
print('Read')
rain_z = radar.fields['reflectivity']['data'].copy()
z_lin = 10.0**(radar.fields['reflectivity']['data']/10.)
rain_z = (z_lin/300.0)**(1./1.4) #Z=300 R1.4
radar.add_field_like('reflectivity', 'rain_z', rain_z, replace_existing = True)
radar.fields['rain_z']['units'] = 'mm/h'
radar.fields['rain_z']['standard_name'] = 'rainfall_rate'
radar.fields['rain_z']['long_name'] = 'rainfall_rate_from_z'
radar.fields['rain_z']['valid_min'] = 0
radar.fields['rain_z']['valid_max'] = 500
print('Gridding')
grids = pyart.map.grid_from_radars(
(radar,), grid_shape=(1, 1001, 1001),
grid_limits=((0, 17000),(-100000, 100000), (-100000, 100000)),
fields=radar.fields.keys(), gridding_algo="map_gates_to_grid",
weighting_function='BARNES')
print('gridded')
long, lat = grids.get_point_longitude_latitude()
height = grids.point_z['data'][:,0,0]
time = np.array([ netCDF4.num2date(grids.time['data'][0], grids.time['units'])])
ds = xarray.Dataset()
for this_field in list(grids.fields.keys()):
this_data = grids.fields[this_field]['data']
my_data = xarray.DataArray(np.expand_dims(this_data,0),
dims = ('time', 'z', 'y', 'x'),
coords = {'time' : (['time'], time),
'z' : (['z'], height),
'lat' :(['y','x'], lat),
'lon' : (['y','x'],long),
'y' : (['y'],lat[:,0]),
'x' : (['x'],long[0,:])})
for this_meta in list(grids.fields[this_field].keys()):
if this_meta is not 'data':
my_data.attrs.update({this_meta: grids.fields[this_field][this_meta]})
ds[this_field] = my_data
ds.lon.attrs = [('long_name', 'longitude of grid cell center'),
('units', 'degrees_east')]
ds.lat.attrs = [('long_name', 'latitude of grid cell center'),
('units', 'degrees_north')]
ds.z.attrs['long_name'] = "height above sea sea level"
ds.z.attrs['units'] = "m"
ds.z.encoding['_FillValue'] = None
ds.lat.encoding['_FillValue'] = None
ds.lon.encoding['_FillValue'] = None
return ds
In [ ]:
def plot_xradar(ds, filename):
lat_lines = np.arange(np.around(ds.lat.min(), decimals=1),
ds.lat.max(), .2)
lon_lines = np.arange(np.around(ds.lon.min(),decimals=1),
ds.lon.max(), .5)
fig = plt.figure(figsize=[10,7])
my_ax = plt.subplot(projection = ccrs.PlateCarree())
rr = ds.rain_z
pc = ds.rain_z.where(rr > 1)[0].sel(z=500, method='nearest').plot.pcolormesh(transform=ccrs.PlateCarree(), ax=my_ax,
x='lon', y='lat', vmin=1,
vmax=100, cmap=pyart.graph.cm.LangRainbow12,
cbar_kwargs = {'label':'Rain Rate (mm/h)'})
my_ax.set_xticks(lon_lines, crs=ccrs.PlateCarree())
my_ax.set_yticks(lat_lines, crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
my_ax.xaxis.set_major_formatter(lon_formatter)
my_ax.yaxis.set_major_formatter(lat_formatter)
gl = my_ax.gridlines(draw_labels=False,
linewidth=2, color='gray', alpha=0.5, linestyle='--')
political_boundaries = cartopy.feature.NaturalEarthFeature(category='cultural',
name='admin_0_boundary_lines_land',
scale='50m', facecolor='none')
states = cartopy.feature.NaturalEarthFeature(category='cultural',
name='admin_1_states_provinces_lines',
scale='50m', facecolor='none')
coast = cartopy.feature.NaturalEarthFeature(category='physical', scale='10m',
facecolor='none', name='coastline')
my_ax.add_feature(political_boundaries, linestyle='-', edgecolor='black')
my_ax.add_feature(states, linestyle='-', edgecolor='black',linewidth=2)
my_ax.add_feature(coast, linestyle='-', edgecolor='black',linewidth=2)
extent = [ds.lon.min(), ds.lon.max(), ds.lat.min(), ds.lat.max()]
my_ax.set_extent(extent)
request = cimgt.GoogleTiles(style='satellite')
my_ax.add_image(request, 10, zorder=0)
plt.savefig(filename)
In [ ]:
def plot_xradar_zoom(ds, filename):
lat_lines = np.arange(np.around(ds.lat.min(), decimals=1),
ds.lat.max(), .05)
lon_lines = np.arange(np.around(ds.lon.min(),decimals=1),
ds.lon.max(), .05)
fig = plt.figure(figsize=[15,10])
my_ax = plt.subplot(projection = ccrs.PlateCarree())
rr = ds.rain_z
pc = ds.rain_z.where(rr > 1)[0].sel(z=500, method='nearest').plot.pcolormesh(transform=ccrs.PlateCarree(), ax=my_ax,
x='lon', y='lat', vmin=1,
vmax=100, cmap=pyart.graph.cm.LangRainbow12,
cbar_kwargs = {'label':'Rain Rate (mm/h)'}, alpha=0.8)
my_ax.set_xticks(lon_lines, crs=ccrs.PlateCarree())
my_ax.set_yticks(lat_lines, crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
my_ax.xaxis.set_major_formatter(lon_formatter)
my_ax.yaxis.set_major_formatter(lat_formatter)
gl = my_ax.gridlines(draw_labels=False,
linewidth=2, color='gray', alpha=0.5, linestyle='--')
political_boundaries = cartopy.feature.NaturalEarthFeature(category='cultural',
name='admin_0_boundary_lines_land',
scale='50m', facecolor='none')
states = cartopy.feature.NaturalEarthFeature(category='cultural',
name='admin_1_states_provinces_lines',
scale='50m', facecolor='none')
coast = cartopy.feature.NaturalEarthFeature(category='physical', scale='10m',
facecolor='none', name='coastline')
my_ax.add_feature(political_boundaries, linestyle='-', edgecolor='black')
my_ax.add_feature(states, linestyle='-', edgecolor='black',linewidth=2)
#my_ax.add_feature(coast, linestyle='-', edgecolor='black',linewidth=2)
extent = [-70.3, -70.2, 43.6, 43.75]
my_ax.set_extent(extent)
request = cimgt.GoogleTiles(style='satellite')
my_ax.add_image(request, 15, zorder=0)
plt.savefig(filename)
In [ ]:
datetime_start = datetime(2015,9,30,13,0)
for ts in range(24):
time_step_needs_doing = True
my_dt = datetime_start + timedelta(minutes=ts*5)
while time_step_needs_doing:
try:
ds = x_array_from_aws('KGYX', my_dt)
time_step_needs_doing = False
except:
print('Slow Down')
time_step_needs_doing = True
fname = datetime.strftime(pandas.to_datetime(ds.time[0].values), '/data/maine/zoom_grid_%Y%m%d_%H%M%S.jpg')
plot_xradar_zoom(ds, fname)
In [ ]:
In [ ]:
datetime_start = datetime(2014,8,13,0,0)
for ts in range(12):
time_step_needs_doing = True
my_dt = datetime_start + timedelta(minutes=ts*10)
while time_step_needs_doing:
try:
ds = x_array_from_aws('KGYX', my_dt)
time_step_needs_doing = False
except:
print('Slow Down')
time_step_needs_doing = True
fname = datetime.strftime(pandas.to_datetime(ds.time[0].values), '/data/maine/grid_%Y%m%d_%H%M%S.jpg')
plot_xradar(ds, fname)
In [ ]: