"DE-BUGGING" NEXRAD LEVEL 2 PLOTTING¶
This tutorial will showcase some simple ways to display NEXRAD LEVEL 2 data in Python. We will start by requesting data from OSDC Griffin's data ID Service. From there we'll plot the data for a given moment in time and then make a time series animation from multiple plots. We'll finish up by separating out "bioscatter" (bugs and birds) from weather using a simple threshold filter.¶
In [1]:
# this will allow our plots to display in this notebook
%matplotlib notebook
# display animations inline
%pylab inline
# our plotting tool
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.image as mpimg
# will need numpy to work with image arrays
import numpy as np
from functools import partial
# we will be plotting NEXRAD LEVEL 2 data with pyart
import pyart
# need to get file location from ID Service JSON
import requests
import json
import os
# to verify hash of downloaded files
import hashlib
# silence annoying deprecation warnings
import warnings
warnings.filterwarnings('ignore', category=DeprecationWarning, module='.*/IPython/.*')
Reading data from ID Service¶
We ran a query in the ID service that returned a list of ARKs. Those ARKs were then stored in the ‘mayfly_arks.txt’ file reference here.¶
In [2]:
# read in .txt file of ark IDs
with open('mayfly_arks.txt', 'r') as f:
file_lines = f.readlines()
id_service_arks = [line.strip() for line in file_lines]
# show part of ark IDs file
!head mayfly_arks.txt
The ID Service stores the data repositories based on the hash of the specified file. We want to download files from the repositories the ID Service points us to. We also want to confirm that the hash of the downloaded file is the same as what the ID Service has stored.¶
In [3]:
# hash provided in signpost should match locally calculated hash
def confirm_hash(hash_algo, file_,actual_hash):
with open(file_) as f:
computed_hash = hash_algo(f.read()).hexdigest()
if computed_hash == actual_hash:
return True
else:
return False
# download, validate(optional) NEXRAD L2 data
def download_from_arks(id_service_arks, intended_dir, hash_confirmation = True, pref_repo='https://griffin-objstore.opensciencedatacloud.org/'):
hash_algo_dict = {'md5':hashlib.md5, 'sha1':hashlib.sha1, 'sha256':hashlib.sha256}
for ark_id in id_service_arks:
signpost_url = 'https://signpost.opensciencedatacloud.org/alias/' + ark_id
resp = requests.get(signpost_url,
proxies={'http':'http://cloud-proxy:3128','https':'http://cloud-proxy:3128'}
)
# make JSON response into dictionary
signpost_dict = resp.json()
# get repository URLs
repo_urls = data_url = signpost_dict['urls']
for url in repo_urls:
# if preferred repo exists, will opt for that URL
if pref_repo in url:
break
# otherwise, will use last url provided
# wow! we can run this bash command from Jupyter!
!sudo wget -P $intended_dir $url
# need file path for hash validation
file_name = url.split('/')[-1]
file_path = os.path.join(intended_dir, file_name)
if hash_confirmation:
# get dict of hash type: hash
hashes = signpost_dict['hashes']
# iterate though list of (hash type, hash) tuples
for hash_tup in hashes.items():
# get proper hash algorithm function
hash_algo = hash_algo_dict[hash_tup[0]]
# fail if not the downloaded file has diff. hash
assert confirm_hash(hash_algo, file_path, hash_tup[1]), '%s hash calculated does not match hash in metadata' % file_path
In [4]:
download_from_arks(id_service_arks, 'mayfly_data')
Let's plot raw reflectivity data¶
In [5]:
# let's wrap our plotting commands in a function for later use
def plot_unfiltered_nexrad(nexrad_file, radar_type='reflectivity'):
# initialize figure
fig = plt.figure(figsize=(7, 10))
ax = fig.add_subplot(111)
radar = pyart.io.read_nexrad_archive(nexrad_file)
# set up pyart display of data
display = pyart.graph.RadarMapDisplay(radar)
# plot data over map of area
display.plot_ppi_map(radar_type,0,
min_lat=radar.latitude['data'][0]-2, max_lat=radar.latitude['data'][0]+2,
min_lon=radar.longitude['data'][0]-2, max_lon=radar.longitude['data'][0]+2,
ax=ax #attach to fig
)
# let's check it out!
plt.show()
return fig
Great! Now let's make multiple plots and save them locally¶
In [6]:
def save_figures(plot_function, radar_files, plot_file_prefix): #plot_file_prefix should be unique or will write over other plots
plot_files = []
count = 1 #we'll use this to make plot file names
for nexrad_file in radar_files:
plot_fig = plot_function(nexrad_file)
plot_file_name = plot_file_prefix+str(count)+'.png'
plot_fig.savefig(plot_file_name, format='png') #save plot as png
count += 1 #update count
plot_files.append(plot_file_name) #add to list of plots
return plot_files
In [7]:
# prepend directory name to file so our function can find
nexrad_file_list = ['mayfly_data/'+file_name for file_name in os.listdir('mayfly_data') if '.gz' in file_name]
# order chronologically
nexrad_file_list.sort()
In [8]:
# run it! will take a few minutes...
plot_files = save_figures(plot_unfiltered_nexrad, nexrad_file_list, 'plots/unfiltered_plot_')
Awesome! Let's take it one step further and animate these plots chronologically¶
In [15]:
# initialize image for animation
def init():
imobj.set_data(np.zeros((100, 100)))
return imobj,
# function to read single plot, display image
def animate(i, plot_files):
file_name = plot_files[i]
img = mpimg.imread(file_name)[-1::-1]
imobj.set_data(img)
plt.axis('off')
return imobj, #needs to be sequence, not single; hence the comma
HOWEVER, animation won't automatically display in Jupyter Notebooks. We'll have to inject HTML with encoded animation to make displaying a non-static visualization possible.¶
In [10]:
# Don't fret! Let's do some fun hacky stuff to circumvent this issue!
from tempfile import NamedTemporaryFile
from IPython.display import HTML
import matplotlib.animation as animation
In the next cell we will 1.) initiate and produce the animation, and 2.) encode the animation and render it in HTML¶
In [14]:
# initialize
fig = plt.figure(figsize=(9, 12))
ax = plt.gca()
imobj = ax.imshow(np.zeros((100, 100)), origin='lower', alpha=1.0, zorder=1, aspect=1)
# produce animation!
anim = animation.FuncAnimation(fig, partial(animate,plot_files=plot_files), init_func=init, repeat = True,
frames=range(len(plot_files)),
interval=10000, blit=True, repeat_delay=1000)
# Big thanks to Jake VanderPlas (@jakevdp) for
# his blog post on animating in Python notebooks:
# http://jakevdp.github.io/blog/2013/05/12/embedding-matplotlib-animations/
# The next two functions were written by him
# HTML video tag we will use to insert animation
VIDEO_TAG = """<video controls>
<source src="data:video/x-m4v;base64,{0}" type="video/mp4">
Your browser does not support the video tag.
</video>"""
# save animation as temporary file, encode
def anim_to_html(anim):
if not hasattr(anim, '_encoded_video'):
with NamedTemporaryFile(suffix='.mp4') as f:
anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264'])
video = open(f.name, "rb").read()
anim._encoded_video = video.encode("base64")
return VIDEO_TAG.format(anim._encoded_video)
# use above function to embed animation in HTML
def display_animation(anim):
plt.close(anim._fig)
return HTML(anim_to_html(anim))
# lights, camera, action!
display_animation(anim)
Out[14]:
Can we filter the reflectivity data for 'bioscatter'?¶
As you can see above, there is something "growing" out of the Wisconsin border. This is in fact a mayfly hatch from the Mississippi River (reported here: http://krocam.com/2015-mayfly-hatch-recorded-on-radar-again/). Additionally the circular spread of "noise" surrounding the center of the map is caused by this similar "bioscatter" phenomena
Using a simple quality control measure loosely based on the "Simple Classifier" described in http://journals.ametsoc.org/doi/abs/10.1175/JTECH-D-13-00073.1 and recommendations of Zachary Flamig of the University of Oklahoma, we will mask data if any of the following are false:
(1.) reflectivity >= 3
(2.) cross correlation ratio >= .99
(3.) differential reflectivity < 2.3In [16]:
# takes 3 numpy arrays, returns array
# where 2 of 3 arrays must have True for
# returned array to have True
def two_true(a,b,c):
return np.logical_or(
np.logical_and(a, np.logical_or(b,c)),
np.logical_and(c, np.logical_or(a,b)),
np.logical_and(b, np.logical_or(c,a))
)
In [17]:
# new function to make plot with new QC reflectivity
# bioscatter variable allows to filter out weather and keep bioscatter (True)
# or filter out bioscatter and keep weather (False)
def plot_filtered_nexrad(nexrad_file, bioscatter):
fig = plt.figure(figsize=(9, 12))
ax = fig.add_subplot(111)
radar = pyart.io.read_nexrad_archive(nexrad_file)
refl_grid = radar.get_field(0, 'reflectivity')
rhohv_grid = radar.get_field(0, 'cross_correlation_ratio')
zdr_grid = radar.get_field(0, 'differential_reflectivity')
# filter the data by making numpy arrays of boolean variables
refl_high = np.greater_equal(refl_grid, 12)
rhohv_high = np.greater_equal(rhohv_grid, .99)
zdr_low = np.less(np.abs(zdr_grid), 2.3)
this_is_weather = two_true(refl_high, rhohv_high,zdr_low)
if bioscatter:
mask = this_is_weather
else:
mask = np.logical_not(this_is_weather)
# mask the data
QC_refl_grid = np.ma.masked_where(mask,refl_grid)
# extract weather sweeps from original radar data
qc = radar.extract_sweeps([0])
# replace reflectivity in original weather sweep with QC reflectivity
qc.add_field_like('reflectivity', 'reflectivityqc', QC_refl_grid)
# set up pyart display of data
display = pyart.graph.RadarMapDisplay(qc)
# plot data over map of area
display.plot_ppi_map('reflectivityqc',0,
min_lat=radar.latitude['data'][0]-2.5, max_lat=radar.latitude['data'][0]+2.5,
min_lon=radar.longitude['data'][0]-2.5, max_lon=radar.longitude['data'][0]+2.5,
lat_0=radar.latitude['data'][0], lon_0=radar.longitude['data'][0],
ax=ax #attach to fig
)
plt.show()
return fig
In [18]:
# save plots with only weather
only_weather_plots = save_figures(partial(plot_filtered_nexrad, bioscatter=False), nexrad_file_list, 'plots/filtered_weather_plot_')
In [19]:
# saves plots with only bioscatter
bioscatter_plots = save_figures(partial(plot_filtered_nexrad, bioscatter=True), nexrad_file_list, 'plots/filtered_bioscatter_plot_')
Let's animate the plots that filter out the bioscatter using the same functions as before¶
In [20]:
fig = plt.figure(figsize=(9, 12))
# 'get current axes' to make for easy plotting
ax = plt.gca()
imobj = ax.imshow( np.zeros((100, 100)), origin='lower', alpha=1.0, zorder=1, aspect=1 )
anim_2 = animation.FuncAnimation(fig, partial(animate,plot_files=only_weather_plots) , init_func=init, repeat = True,
frames=range(len(only_weather_plots)),
interval=10000, blit=True, repeat_delay=1000)
display_animation(anim_2)
Out[20]:
Given our rather rudimentary filter, we did a decent job displaying only the weather. One noticeable and interesting defect is that we filtered out the outskirts of the thunderstorm cell (will be more obvious in next animation). This is likely hail or very heavy rain that our filter incorrectly classifies as bioscatter.¶
Finally, let's animate the plots that only filter out the weather, retaining the bioscatter¶
In [21]:
fig = plt.figure(figsize=(9, 12))
# 'get current axes' to make for easy plotting
ax = plt.gca()
imobj = ax.imshow( np.zeros((100, 100)), origin='lower', alpha=1.0, zorder=1, aspect=1 )
anim_3 = animation.FuncAnimation(fig, partial(animate,plot_files=bioscatter_plots) , init_func=init, repeat = True,
frames=range(len(bioscatter_plots)),
interval=10000, blit=True, repeat_delay=1000)
display_animation(anim_3)
Out[21]: