How to isolate valid echoes and grid without griding to an artifical boundary or allow hash to raise background
Scott Collis scollis@anl.gov
The Question: How do we identify cloud and precipitation echoes and isolate them on a moving noise background as the noise floor of reflectivity rises as the square root of the distance from the radar (cone) and make it amenable to griding (ie at a valid return edge what do you grid down to?).. This extends apon previous work covered in the class at the ASR PI meeting by looking at what happens if we replace the "non-significant" gates with the minimum detectable signal instead of just masking before we grid...
import pyart
from matplotlib import pyplot as plt
import numpy as np
from matplotlib import rc
from scipy.optimize import curve_fit
from copy import deepcopy
%matplotlib inline
radar = pyart.io.read('/data/radar/nsaxsaprppiC1.a1.20140201.184802.nc')
#Radar gate spacing is wrong in the file is should be 60m
radar.range['data'] = 60.0*radar.range['data']/50.0
print radar.fields.keys()
print radar.fixed_angle['data']
print radar.range['data'].max()
print radar.longitude['data']
print radar.latitude['data']
Ok, first we need to work out what the minimum detectable signal is.. If you look down the page you will see the line of snow is not covering the whole domain.. since the Time Range array of reflectivity covers all sweeps lets just take the min in each ray (axis = 0 means take the min alon axis 0, ie the time axis)
f = plt.figure(figsize = [10,5])
mds = radar.fields['reflectivity_horizontal']['data'].min(axis = 0)
plt.plot(radar.range['data']/1000.0, mds)
plt.xlabel('Range from antenna (km)')
plt.ylabel('Range minimum reflectivity (dBZ)')
Cool bananas! So besides some omnipresent returns close to the radar we have this nice $Z_{min} = a * \sqrt{r} + Z_{min}\bigg|_{r=0}$ form.. Lets fit.. a=fac below and $Z_{min}\bigg|_{r=0}$ = zr
def func(rng, zr, fac):
return fac*np.sqrt(rng) + zr
popt, pcov = curve_fit(func, radar.range['data'][200::], mds[200::])
print popt
print pcov
So the above is the fitting variables and thier covariences.. Life is good! Note we started at the 200th gate in order to avoid the detritus close to the radar. Now lets investigate our fit..
f = plt.figure(figsize = [10,5])
plt.plot(radar.range['data']/1000.0, mds)
plt.plot(radar.range['data']/1000.0, func(radar.range['data'], popt[0], popt[1]))
plt.xlabel('Range from antenna (km)')
plt.ylabel('Range minimum reflectivity (dBZ)')
Not perfect.. Perhpas around 10km we have some signal that corrupts the fit.. I'll fix later! Lets create a moment like 2d field for later...
mds_2d = func(np.meshgrid(radar.range['data'], radar.time['data'])[0], popt[0], popt[1])
Lets plot out raw Z...
font = {'size' : 16}
rc('font', **font)
f = plt.figure(figsize=[15,8])
myd = pyart.graph.RadarMapDisplay(radar)
myd.plot_ppi_map('reflectivity_horizontal', 1, vmin=-20, vmax=20,
min_lon=-157.1, max_lon=-156, min_lat=71.2, max_lat=71.6,
lon_lines = np.arange(-158, -154, .2),
lat_lines = np.arange(69, 72, .1), resolution = 'h',
auto_range=False)
myd.plot_range_ring(10. * 1000.0, line_style = 'k-')
myd.plot_range_ring(20. * 1000.0, line_style = 'k--')
myd.plot_range_ring(30. * 1000.0, line_style = 'k-')
myd.plot_range_ring(40. * 1000.0, line_style = 'k--')
myd.plot_line_xy([-40000.0, 40000.0], [0.0,0.0], line_style = 'k-')
myd.plot_line_xy([0.0,0.0],[-20000.0, 200000.0], line_style = 'k-')
myd.plot_point(radar.longitude['data'], radar.latitude['data'])
ax1 = plt.gca()
x0,x1 = ax1.get_xlim()
y0,y1 = ax1.get_ylim()
ax1.set_aspect((x1-x0)/(y1-y0) - .05)
#plt.text(xp, yp, 'foo')
Holy Ice Clutter Batman!
lets create some classifications based off purely polarimetric data...
# So there are a variety of things we can classify off.. And the nice thing is Numpy's where
# function allows great filtering basded on logic!
# the first two lines create an array of Trues and Falses
is_coherent = radar.fields['normalized_coherent_power']['data'] > 0.5
is_correlated = radar.fields['cross_correlation_ratio']['data'] > 0.8
# Now we are going to do some array based logic math..
is_correlated_and_coherent = np.logical_and(is_coherent, is_correlated)
is_correlated_and_not_coherent = np.logical_and(np.logical_not(is_coherent), is_correlated)
is_not_correlated_and_coherent = np.logical_and(is_coherent, np.logical_not(is_correlated))
is_not_correlated_and_not_coherent = np.logical_and(np.logical_not(is_coherent),
np.logical_not(is_correlated))
# initialize a blank array
logical_classes = np.zeros_like(radar.fields['normalized_coherent_power']['data'])
# Now we are going to fill in regions (allowing overwriting)
# to create a numerical based classificatiob
logical_classes[np.where(is_correlated_and_coherent)] = 1.0
logical_classes[np.where(is_correlated_and_not_coherent)] = 2.0
logical_classes[np.where(is_not_correlated_and_coherent)] = 3.0
logical_classes[np.where(is_not_correlated_and_not_coherent)] = 4.0
# add the field to the radar object
field_dict = {'data' : logical_classes,
'units' : 'Unitless',
'long_name' : 'Logical Classifications based on RhoHV and NCP/SQ',
'standard_name' : "Logical_Class"}
radar.add_field('logical_class', field_dict)
Ok.. now for the texture calculations.. Basically we are going to calculate:
$texture = \frac{<V_r>}{ RMS( V_r)}$
where
$texture = <(V_r - <V_r>)^2>$
where the angluar brackets denote an 11 point running filter. In this case to make it simpler we are going to use the filter as defined in this line of the Py-ART Source code . Essentially it is an 11 point hanning filter convolution.. The tricky bit I am neglecting to include here involves keeping the length of the array constant.. See if you can follow the in-line comments. At the end of the calculation, as before, we simply create a field dictionary and add this back into the radar object.
# time for some texture fields
SNR_VR = np.zeros_like(radar.fields['mean_doppler_velocity']['data'])
for i in range(SNR_VR.shape[0]):
this_ray_of_data = radar.fields['mean_doppler_velocity']['data'][i,:]
signal = pyart.correct.phase_proc.smooth_and_trim(this_ray_of_data)
noise = pyart.correct.phase_proc.smooth_and_trim(np.sqrt((this_ray_of_data - signal) ** 2))
this_ray_SNR = np.abs(signal) / noise
SNR_VR[i,:] = this_ray_SNR
field_dict = {'data' : SNR_VR,
'units' : 'Unitless',
'long_name' : 'Signal to Noise ratio for mean_doppler_velocity',
'standard_name' : "SNR_VR"}
radar.add_field('SNR_VR', field_dict)
Lets use a threshold for SNR of 5 and we can use Numpy's Masked Array methods to create fields with non-meteo returns flagged. Follow the in-line comments below
# The below will create an arry of the same shape of SNR_VR but with Trues where SNR < 10.0
# and Falses where NOT(SNR <10.0)
is_messy = SNR_VR < 5.0
# Now to Numpy's masked array.. the np.ma module has a variety of methods.. here
# we will use the masked_where method. This is invoked:
# masked_where(condition, data)
# and everywhere where condition = True will be marked as masked in the resultant array
vr_masked_Z = np.ma.masked_where(is_messy, radar.fields['reflectivity_horizontal']['data'])
# Ok.. The below is just the usual creating of a field ready to go back into our radar object.
# instead of manual entry we simply copy the metadata from the original Z field
field_dict = {'data' : vr_masked_Z}
for key in ['_FillValue', 'long_name', 'units', 'valid_min', 'valid_max',
'standard_name', 'coordinates']:
field_dict.update({key : radar.fields['reflectivity_horizontal'][key]})
# and BOOM! Just add it back into our radar object!
radar.add_field('vr_filtered_reflectivity_horizontal', field_dict)
Ok! That was cool! Now lets do exactly the same but using Normalized Coherent Power as our condition..
ncp_masked_Z = np.ma.masked_where(np.logical_not(is_coherent),
radar.fields['reflectivity_horizontal']['data'])
field_dict = {'data' : ncp_masked_Z}
for key in ['_FillValue', 'long_name', 'units', 'valid_min', 'valid_max',
'standard_name', 'coordinates']:
field_dict.update({key : radar.fields['reflectivity_horizontal'][key]})
radar.add_field('ncp_filtered_reflectivity_horizontal', field_dict)
OK! Now as many of you reading may know one of the things that keep me up at night is the mapping of radar moments.. See if you hit the edge of a valid return in Azimuth, elevation and range space when using a radius of influence type objective analysis code (ie Trapp and Doswell, eg Barnes, Cressman) then what do you interpolate down to? Furthermore, if you have isolated "missed" screening out at range, where there is increased azimuthal seperation and you want a larger radius of influence (ROI) in order to avoid artifactes (Collis' JTECH article) the ROI will map these gates to a nice neat sphere.. If however these spheres were encompassed by signal levels at the Minimum Descernable Signal (MDS) Level these would integrate to the MDS level by pure fact that the sample density along a ray is very large. So the next few lines, rather than masking gates we replace the values we consider to be "non significant" with the MDS.. this gives something to interpolate down to that is instrumental based. That is this is what a radar would see, given our goal of non (first trip) meteorological echo radar data.. that is a radar that did not suffer from clutter or second trip echoes.. an impossible goal but something to strive for
# The below will create an arry of the same shape of SNR_VR but with Trues where SNR < 10.0
# and Falses where NOT(SNR <10.0)
is_messy = SNR_VR < 5.0
#instead of masking we replace with MDS
vr_repl_Z = deepcopy(radar.fields['reflectivity_horizontal']['data'])
vr_repl_Z[np.where(is_messy)] = mds_2d[np.where(is_messy)]
# Ok.. The below is just the usual creating of a field ready to go back into our radar object.
# instead of manual entry we simply copy the metadata from the original Z field
field_dict = {'data' : vr_repl_Z}
for key in ['_FillValue', 'long_name', 'units', 'valid_min', 'valid_max',
'standard_name', 'coordinates']:
field_dict.update({key : radar.fields['reflectivity_horizontal'][key]})
# and BOOM! Just add it back into our radar object!
radar.add_field('vr_replaced_reflectivity_horizontal', field_dict)
Now we are going to visualize them.. first the SNR (sorry SG) masked Z
font = {'size' : 16}
rc('font', **font)
def pretty_ppi(radar, moment):
f = plt.figure(figsize=[15,8])
myd = pyart.graph.RadarMapDisplay(radar)
myd.plot_ppi_map(moment, 1, vmin=-20, vmax=20,
min_lon=-157.1, max_lon=-156, min_lat=71.2, max_lat=71.6,
lon_lines = np.arange(-158, -154, .2),
lat_lines = np.arange(69, 72, .1), resolution = 'h',
auto_range=False)
ax1 = plt.gca()
x0,x1 = ax1.get_xlim()
y0,y1 = ax1.get_ylim()
ax1.set_aspect((x1-x0)/(y1-y0) - .15)
pretty_ppi(radar, 'vr_filtered_reflectivity_horizontal')
So we an see that the VR texture mask gets rid of not only non-coherent signals but also nukes a heap of the clutter, including the ice clutter! Nice thing is when we grid the ROI will naturally bring $Z_e$ from higher tilts down, to a limit of course.. Lets see how the NCP filitered $Z_e$ looks:
pretty_ppi(radar, 'ncp_filtered_reflectivity_horizontal')
Hmm.. so yes, since NCP is the ratio of the (phase/velocity space) coherent signal to the non-coherent the clutter, which is, at times, very coherent remains.. So for our mapping experiment lets stick with the texture masked reflectivity..
pretty_ppi(radar,'vr_replaced_reflectivity_horizontal')
So above is our "replaced" Z field .. Everywhere deemed to be a non-valid return in backfilled with our estimated MDS.. As you can see it is a bit low as we used the minimum.. In practice we would use an empty scene, non-transmitting or just calculate from the noise floor, radar equation and number of integrations.. Lets GRID!
grid = pyart.map.grid_from_radars(
(radar,),
grid_shape=(31., 151, 151),
grid_limits=((0, 5000.), (-25000., 25000.), (-25000., 25000.)),
fields=radar.fields.keys(), refl_field='reflectivity_horizontal')
So this returns a Python Grid object, lets write a neat wrapper to the plotting with maps module and slice and dice the original raw reflectity some ways from the radar:
def nice_map(grid, moment, lat = 71.4, lon = -156.8, level = 5):
# panel sizes
map_panel_axes = [0.05, 0.05, .4, .80]
x_cut_panel_axes = [0.55, 0.10, .4, .30]
y_cut_panel_axes = [0.55, 0.50, .4, .30]
colorbar_panel_axes = [0.05, 0.90, .4, .03]
#min_lon=-157.1, max_lon=-156, min_lat=71.2, max_lat=71.6
# parameters
vmin = -20
vmax = 20
fig = plt.figure(figsize=[15, 8])
grid_display = pyart.graph.GridMapDisplay(grid)
ax1 = fig.add_axes(map_panel_axes)
grid_display.plot_basemap(resolution = 'h',
lon_lines = np.arange(-158, -154, .2),
lat_lines = np.arange(69, 72, .1))
grid_display.plot_grid(moment,
level=level, vmin=vmin, vmax=vmax)
grid_display.plot_crosshairs(lon=lon, lat=lat)
# colorbar
cbax = fig.add_axes(colorbar_panel_axes)
grid_display.plot_colorbar(cax=cbax)
# panel 2, longitude slice.
ax2 = fig.add_axes(x_cut_panel_axes)
grid_display.plot_longitude_slice(moment,
vmin=vmin, vmax=vmax, lon=lon, lat=lat)
ax2.set_xlabel('Distance from radar (km)')
# panel 3, latitude slice
ax3 = fig.add_axes(y_cut_panel_axes)
grid_display.plot_latitude_slice(moment,
vmin=vmin, vmax=vmax, lon=lon, lat=lat)
x0,x1 = ax1.get_xlim()
y0,y1 = ax1.get_ylim()
ax1.set_aspect((x1-x0)/(y1-y0))
nice_map(grid, 'reflectivity_horizontal')
Nice morphology! as you can see we get a cone of silence in the constant longitude slice, we also have an extra tilt up way high that is so far from the -2th tilt we can not interpolate between (and should not!). Lets see what happens whem we map the NCP filtered $Z_e$:
nice_map(grid,'ncp_filtered_reflectivity_horizontal')
oooh.. that is not nice! Super high $Z_e$ where there are really no returns! What causes this? It is the occasional clutter or hash noise not filtered out which, due to the wide ROI out at range and lack of points to interpolate down to basically creates a ROI sized paint splodge.. This is the issue with not having a nice low background to smooth too.. Lets see what happens when we map using our replaced field..
nice_map(grid,'vr_replaced_reflectivity_horizontal')
Pretty! Now since we have a background of MDS values when the mapping sphere has a missed gate of junk its influence is minimized by a large number of low valued gates.. As an added benefit Z returns naturally map down to the MDS rather than having a convolution between the last gates and the ROI.. But lets really challenge the map.. lets look at the worst possible place, right near the radar at the lowest surface:
nice_map(grid,'vr_replaced_reflectivity_horizontal',
lat = radar.latitude['data'], lon = radar.longitude['data'],
level = 0)
Wonderful! Great isolation, no sea ice clutter.. But yes.. MDS is a bit low.. So what does the difference look like? Lets look at the same slices with the original Z:
nice_map(grid,'reflectivity_horizontal',
lat = radar.latitude['data'], lon = radar.longitude['data'],
level = 0)
yep! There is our clutter.. Lets look at them side by side:
def nice_map_two(grid, moment,moment2, lat = 71.4, lon = -156.8, level = 5):
# panel sizes
map_panel_axes = [0.05, 0.55, .4, .30]
map2_panel_axes = [0.55, 0.55, .4, .30]
x_cut_panel_axes1 = [0.05, 0.05, .4, .20]
y_cut_panel_axes1 = [0.05, 0.30, .4, .20]
x_cut_panel_axes2 = [0.55, 0.05, .4, .20]
y_cut_panel_axes2 = [0.55, 0.30, .4, .20]
colorbar_panel_axes = [0.05, 0.90, .4, .03]
# parameters
vmin = -20
vmax = 20
fig = plt.figure(figsize=[15, 15])
grid_display = pyart.graph.GridMapDisplay(grid)
ax1 = fig.add_axes(map_panel_axes)
grid_display.plot_basemap(resolution = 'h',
lon_lines = np.arange(-158, -154, .4),
lat_lines = np.arange(69, 72, .1))
grid_display.plot_grid(moment,
level=level, vmin=vmin, vmax=vmax)
grid_display.plot_crosshairs(lon=lon, lat=lat)
# colorbar
cbax = fig.add_axes(colorbar_panel_axes)
grid_display.plot_colorbar(cax=cbax)
ax2 = fig.add_axes(map2_panel_axes)
grid_display.plot_basemap(resolution = 'h',
lon_lines = np.arange(-158, -154, .4),
lat_lines = np.arange(69, 72, .1))
grid_display.plot_grid(moment2,
level=level, vmin=vmin, vmax=vmax)
grid_display.plot_crosshairs(lon=lon, lat=lat)
# panel 3, longitude slice.
ax3 = fig.add_axes(x_cut_panel_axes1)
grid_display.plot_longitude_slice(moment,
vmin=vmin, vmax=vmax, lon=lon, lat=lat)
ax3.set_xlabel('Distance from radar (km)')
plt.xlim([-25,25])
# panel 4, latitude slice
ax4 = fig.add_axes(y_cut_panel_axes1)
grid_display.plot_latitude_slice(moment,
vmin=vmin, vmax=vmax, lon=lon, lat=lat)
plt.xlim([-25,25])
# panel 5, longitude slice.
ax5 = fig.add_axes(x_cut_panel_axes2)
grid_display.plot_longitude_slice(moment2,
vmin=vmin, vmax=vmax, lon=lon, lat=lat)
ax5.set_xlabel('Distance from radar (km)')
plt.xlim([-25,25])
# panel 6, latitude slice
ax6 = fig.add_axes(y_cut_panel_axes2)
grid_display.plot_latitude_slice(moment2,
vmin=vmin, vmax=vmax, lon=lon, lat=lat)
plt.xlim([-25,25])
#set aspect ratio
x0,x1 = ax1.get_xlim()
y0,y1 = ax1.get_ylim()
ax1.set_aspect((x1-x0)/(y1-y0))
x0,x1 = ax2.get_xlim()
y0,y1 = ax2.get_ylim()
ax2.set_aspect((x1-x0)/(y1-y0))
nice_map_two(grid, 'reflectivity_horizontal','vr_replaced_reflectivity_horizontal',
lat = radar.latitude['data'], lon = radar.longitude['data'],
level = 0)
plt.savefig('/Users/scollis/JohnHunter.pdf')
So we can clearly see we need a better MDS formulaton but the mapped moments are already of far higher quality than they were.. So progress towards CMAC2 and MMGC2!
Thus concludes the Tutorial! Questions? Comments? Science Lead: Scott Collis Development lead: Jonathan Helmus.