Welcome to the CSU_RadarTools demonstration IPython notebook. CSU_RadarTools is a collection of open-source tools for weather radar data quality control and analysis, written in the Python language. The package collates a disparate number of tools that have been developed over many years at Colorado State University. The purpose of this notebook is to demonstrate how to use the package to accomplish various tasks.
In order to get started, in addition to making sure you have a robust Python installation with most major tools (e.g., numpy, matplotlib, pandas, etc.) please download and install the following modules:
In [1]:
from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import pyart
import glob
from pyart.io.common import radar_coords_to_cart
from skewt import SkewT
from csu_radartools import (csu_fhc, csu_liquid_ice_mass, csu_blended_rain,
csu_dsd, csu_kdp, csu_misc)
%matplotlib inline
First things first, you see that csu_radartools is actually several different sub-modules:
Let's examine each sub-module in turn. First, however, let's define a simple function to do side-by-side PPI radar plots, powered by Py-ART.
In [2]:
def two_panel_plot(radar, sweep=0, var1='reflectivity', vmin1=0, vmax1=65, cmap1='RdYlBu_r',
units1='dBZ', var2='differential_reflectivity', vmin2=-5, vmax2=5,
cmap2='RdYlBu_r', units2='dB', return_flag=False, xlim=[-150,150],
ylim=[-150,150]):
display = pyart.graph.RadarDisplay(radar)
fig = plt.figure(figsize=(13,5))
ax1 = fig.add_subplot(121)
display.plot_ppi(var1, sweep=sweep, vmin=vmin1, vmax=vmax1, cmap=cmap1,
colorbar_label=units1, mask_outside=True)
display.set_limits(xlim=xlim, ylim=ylim)
ax2 = fig.add_subplot(122)
display.plot_ppi(var2, sweep=sweep, vmin=vmin2, vmax=vmax2, cmap=cmap2,
colorbar_label=units2, mask_outside=True)
display.set_limits(xlim=xlim, ylim=ylim)
if return_flag:
return fig, ax1, ax2, display
This package currently works with X-, C-, and S-band polarimetric radar data. Best results are obtained if you use a temperature sounding along with polarimetric radar data.
In [3]:
#Read in the data
sndfile = '/Users/tjlang/Documents/OVWST/CPOL/soundings/snd_Darwin.txt'
radarfile = '/Users/tjlang/Documents/OVWST/CPOL/output/20060119/' + \
'cfrad.20060119_170029.000_to_20060121_020810.000_CPOL_v1_PPI.nc'
radar = pyart.io.read(radarfile)
print(radar.fields.keys())
sounding = SkewT.Sounding(sndfile)
This volume is from the CPOL C-band polarimteric Doppler radar. The fields of note are corrected reflectivity (ZC), differential reflectivity (ZD), specific differential phase (KD), and correlation coefficient (RH). The sounding is a nearby sounding in the UWyo format (i.e., from http://weather.uwyo.edu/upperair/sounding.html).
The way CSU_RadarTools is designed is it works on Python arrays or scalars. This allows you to use it with any kind of radar data, whether it was ingested via Py-ART or something else, or is in polar coordinate or gridded form. Thus, we first need to extract the relevant fields from the Py-ART radar object.
In [4]:
dz = radar.fields['ZC']['data']
dr = radar.fields['ZD']['data']
kd = radar.fields['KD']['data']
rh = radar.fields['RH']['data']
But we also need to get the sounding data interpolated onto the same structure as the radar data. Here is a way to do that via Py-ART and numpy.interp().
In [5]:
def get_z_from_radar(radar):
"""Input radar object, return z from radar (km, 2D)"""
azimuth_1D = radar.azimuth['data']
elevation_1D = radar.elevation['data']
srange_1D = radar.range['data']
sr_2d, az_2d = np.meshgrid(srange_1D, azimuth_1D)
el_2d = np.meshgrid(srange_1D, elevation_1D)[1]
xx, yy, zz = radar_coords_to_cart(sr_2d/1000.0, az_2d, el_2d)
return zz + radar.altitude['data']
def check_sounding_for_montonic(sounding):
"""
So the sounding interpolation doesn't fail, force the sounding to behave
monotonically so that z always increases. This eliminates data from
descending balloons.
"""
snd_T = sounding.soundingdata['temp'] # In old SkewT, was sounding.data
snd_z = sounding.soundingdata['hght'] # In old SkewT, was sounding.data
dummy_z = []
dummy_T = []
if not snd_T.mask[0]: #May cause issue for specific soundings
dummy_z.append(snd_z[0])
dummy_T.append(snd_T[0])
for i, height in enumerate(snd_z):
if i > 0:
if snd_z[i] > snd_z[i-1] and not snd_T.mask[i]:
dummy_z.append(snd_z[i])
dummy_T.append(snd_T[i])
snd_z = np.array(dummy_z)
snd_T = np.array(dummy_T)
return snd_T, snd_z
def interpolate_sounding_to_radar(sounding, radar):
"""Takes sounding data and interpolates it to every radar gate."""
radar_z = get_z_from_radar(radar)
radar_T = None
snd_T, snd_z = check_sounding_for_montonic(sounding)
shape = np.shape(radar_z)
rad_z1d = radar_z.ravel()
rad_T1d = np.interp(rad_z1d, snd_z, snd_T)
return np.reshape(rad_T1d, shape), radar_z
radar_T, radar_z = interpolate_sounding_to_radar(sounding, radar)
And now we are ready to run the HID algorithm. Behold ...
In [6]:
scores = csu_fhc.csu_fhc_summer(dz=dz, zdr=dr, rho=rh, kdp=kd, use_temp=True, band='C',
T=radar_T)
fh = np.argmax(scores, axis=0) + 1
To enable the ability to find out the second-ranked (or third, etc.) species, csu_fhc_summer() returns the scores for all the different categories, not just the max. So to get the traditional HID category number you have to use numpy.argmax() as above. The summer HID from CSU returns 10 possible categories:
And these are represented as integers in the newly created fh array, which as the same structure as dz, dr, etc. We'd like to plot these data using Py-ART, which means we need to turn fh in a radar object field. Let's do that.
In [7]:
def add_field_to_radar_object(field, radar, field_name='FH', units='unitless',
long_name='Hydrometeor ID', standard_name='Hydrometeor ID',
dz_field='ZC'):
"""
Adds a newly created field to the Py-ART radar object. If reflectivity is a masked array,
make the new field masked the same as reflectivity.
"""
masked_field = np.ma.asanyarray(field)
fill_value = -32768
if hasattr(radar.fields[dz_field]['data'], 'mask'):
setattr(masked_field, 'mask', radar.fields[dz_field]['data'].mask)
fill_value = radar.fields[dz_field]['_FillValue']
field_dict = {'data': masked_field,
'units': units,
'long_name': long_name,
'standard_name': standard_name,
'_FillValue': fill_value}
radar.add_field(field_name, field_dict, replace_existing=True)
return radar
In [8]:
radar = add_field_to_radar_object(fh, radar)
And now let's plot it up! First, however, let's fix Py-ART's colorbars to look nice for HID and other category-style fields.
In [9]:
hid_colors = ['White', 'LightBlue', 'MediumBlue', 'DarkOrange', 'LightPink',
'Cyan', 'DarkGray', 'Lime', 'Yellow', 'Red', 'Fuchsia']
cmaphid = colors.ListedColormap(hid_colors)
cmapmeth = colors.ListedColormap(hid_colors[0:6])
def adjust_fhc_colorbar_for_pyart(cb):
cb.set_ticks(np.arange(1.4, 10, 0.9))
cb.ax.set_yticklabels(['Drizzle', 'Rain', 'Ice Crystals', 'Aggregates',
'Wet Snow', 'Vertical Ice', 'LD Graupel',
'HD Graupel', 'Hail', 'Big Drops'])
cb.ax.set_ylabel('')
cb.ax.tick_params(length=0)
return cb
def adjust_meth_colorbar_for_pyart(cb):
cb.set_ticks(np.arange(1.25, 5, 0.833))
cb.ax.set_yticklabels(['R(Kdp, Zdr)', 'R(Kdp)', 'R(Z, Zdr)', 'R(Z)', 'R(Zrain)'])
cb.ax.set_ylabel('')
cb.ax.tick_params(length=0)
return cb
In [10]:
# Actual plotting done here
lim = [-80, 80]
fig, ax1, ax2, display = two_panel_plot(radar, sweep=5, var1='ZC', var2='FH', vmin2=0,
vmax2=10, cmap2=cmaphid, units2='', return_flag=True,
xlim=lim, ylim=lim)
display.cbs[1] = adjust_fhc_colorbar_for_pyart(display.cbs[1])
The HID algorithm, like most algorithms in CSU_RadarTools, works on scalars too. So if you are just curious what category a set of polarimetric values will get, try the following:
In [11]:
scores = csu_fhc.csu_fhc_summer(dz=45.0, zdr=0.0, kdp=-0.2, rho=0.95, T=-1)
print(np.argmax(scores, axis=0) + 1)
... which is high-density graupel.
In [12]:
#Here's help to see what else you can do with csu_fhc.csu_fhc_summer
help(csu_fhc.csu_fhc_summer)
In this module, we use difference reflectivity (ZDP) to derive separate ice and liquid water contributions to reflectivity, then use simple Z-M relationships to estimate mass contents. We can build on the above by using the same radar volume, sounding, etc.
In [13]:
# Function expects, reflectivity, differential reflectivity,
# and altitude (km MSL) at a minimum.
# Temperature is optional.
mw, mi = csu_liquid_ice_mass.calc_liquid_ice_mass(dz, dr, radar_z/1000.0, T=radar_T)
radar = add_field_to_radar_object(mw, radar, field_name='MW', units='g m-3',
long_name='Liquid Water Mass',
standard_name='Liquid Water Mass')
radar = add_field_to_radar_object(mi, radar, field_name='MI', units='g m-3',
long_name='Ice Water Mass',
standard_name='Ice Water Mass')
In [14]:
two_panel_plot(radar, sweep=5, var1='MW', var2='MI', vmin1=0, vmax1=10, vmin2=0, vmax2=10,
cmap1='YlOrRd', cmap2='YlOrRd', units1='g m-3',
units2='g m-3', xlim=lim, ylim=lim)
In [15]:
#Here's help to see what else you can do with csu_liquid_ice_mass.calc_liquid_ice_mass
help(csu_liquid_ice_mass.calc_liquid_ice_mass)
This module provides access to two distinct CSU blended rainfall algorithm methodologies. The first function is calc_blended_rain(), whcih is the classic version that utilizes ZDP to estimate ice fraction, and the second function is calc_hidro_rain(), which is newer and uses HID to help distinguish regions of ice.
There are also a number of individual rain rate calculation functions (e.g., Z-R, R(Kdp), R(Z,Zdr), etc.) provided with this module. The user can specify the coefficients or just use the defaults, which come from standard literature references.
Output from the main blended functions include the method used at each gate, which the user can also plot up. For the example below we show the results from calc_hidro_rain(). Since we already calculated HID, we can use our previous output from csu_fhc to supply that field.
In [16]:
rain, method = csu_blended_rain.csu_hidro_rain(dz=dz, zdr=dr, kdp=kd, fhc=fh)
radar = add_field_to_radar_object(rain, radar, field_name='rain', units='mm h-1',
long_name='HIDRO Rainfall Rate',
standard_name='Rainfall Rate')
radar = add_field_to_radar_object(method, radar, field_name='method', units='',
long_name='HIDRO Rainfall Method',
standard_name='Rainfall Method')
In [17]:
fig, ax1, ax2, display = two_panel_plot(radar, sweep=0, var1='rain', vmin1=0, vmax1=150,
cmap1='GnBu', var2='method', vmin2=0, vmax2=5,
cmap2=cmapmeth, units2='', return_flag=True,
xlim=lim, ylim=lim, units1='mm h-1')
display.cbs[1] = adjust_meth_colorbar_for_pyart(display.cbs[1])
The method variable is an integer output that matches the following legend:
Method Legend
The last method is only used by calc_blended_rain(). Here is how to use that algorithm. It has additional outputs, which include ZDP and ice fraction (unitless), if you set ice_flag = True.
In [18]:
rain, method, zdp, fi = csu_blended_rain.calc_blended_rain(dz=dz, zdr=dr,
kdp=kd, ice_flag=True)
radar = add_field_to_radar_object(rain, radar, field_name='rain_blend', units='mm h-1',
long_name='Blended Rainfall Rate',
standard_name='Rainfall Rate')
radar = add_field_to_radar_object(method, radar, field_name='method_blend', units='',
long_name='Blended Rainfall Method',
standard_name='Rainfall Method')
radar = add_field_to_radar_object(zdp, radar, field_name='ZDP', units='dB',
long_name='Difference Reflectivity',
standard_name='Difference Reflectivity')
radar = add_field_to_radar_object(fi, radar, field_name='FI', units='',
long_name='Ice Fraction',
standard_name='Ice Fraction')
In [19]:
fig, ax1, ax2, display = two_panel_plot(radar, sweep=0, var1='rain_blend', vmin1=0, vmax1=150,
cmap1='GnBu', var2='method_blend', vmin2=0, vmax2=5,
cmap2=cmapmeth, units2='', return_flag=True, xlim=lim,
ylim=lim, units1='mm h-1')
display.cbs[1] = adjust_meth_colorbar_for_pyart(display.cbs[1])
We can plot ZDP and FI too, to gauge how they are influencing the rainfall method choices seen above.
In [20]:
two_panel_plot(radar, sweep=0, var1='ZDP', units1='dB', vmin1=0, vmax1=65,
cmap1='cubehelix', var2='FI', vmin2=0, vmax2=1,
cmap2='YlOrRd_r', units2='', xlim=lim, ylim=lim)
In [21]:
#And of course, we can get help on all the other stuff!
#Original IDL notes are listed toward the top, with the Python-specific notes below
help(csu_blended_rain)
If you dig into the docstrings above, you see individual functions for rainfall, like R(Kdp) and R(Z,Zdr). So you can just access those on your own, if you like.
In [22]:
rkdp = csu_blended_rain.calc_rain_kdp(kd)
radar = add_field_to_radar_object(rkdp, radar, field_name='rkdp', units='mm h-1',
long_name='Rainfall Rate R(Kdp)',
standard_name='Rainfall Rate')
In [23]:
#Compare R(Kdp) with Kdp!
two_panel_plot(radar, sweep=0, var1='rkdp', vmin1=0, vmax1=150,
cmap1='GnBu', units1='mm h-1',
var2='KD', vmin2=-5, vmax2=5, cmap2='RdYlBu',
units2='deg km-1', xlim=lim, ylim=lim)
This module enables retrieval of common drop-size distribution (DSD) parameters. As with most aspects of CSU_RadarTools, the end user is responsible for proper masking of their data, as DSD retrievals are only valid in pure rain. Remember, garbage in = garbage out! The module is designed to work with C-band or S-band, and the retrievals are different depending on radar frequency. The main function is calc_dsd() and it returns median volume diameter (D0), normalized intercept parameter (Nw), and mu assuming a gamma distribution DSD model.
Let's first check things out for C-band.
In [24]:
d0, Nw, mu = csu_dsd.calc_dsd(dz=dz, zdr=dr, kdp=kd, band='C')
radar = add_field_to_radar_object(d0, radar, field_name='D0', units='mm',
long_name='Median Volume Diameter',
standard_name='Median Volume Diameter')
logNw = np.log10(Nw)
radar = add_field_to_radar_object(logNw, radar, field_name='NW', units='',
long_name='Normalized Intercept Parameter',
standard_name='Normalized Intercept Parameter')
radar = add_field_to_radar_object(mu, radar, field_name='MU', units='',
long_name='Mu', standard_name='Mu')
The C-band retrievals, based on Bringi et al. (2009), don't bother retrieving mu. So CSU_RadarTools will fix it at 3 unless the end user changes the csu_dsd.DEFAULT_MU global variable. So let's check out instead what happened with D0 and Nw.
In [25]:
two_panel_plot(radar, sweep=0, var1='D0', vmin1=0, vmax1=2, cmap1='GnBu', units1='mm',
var2='NW', vmin2=0, vmax2=8, cmap2='cubehelix', units2='log10(Nw)', xlim=lim,
ylim=lim)
Note the differences in D0 & Nw pairs inside and outside of convection. Both here and above with ice fraction, which also depends heavily on Zdr, there are some possible differential attenuation impacts on the retrievals downrange of one of the cores just west of the radar.
What about S-band DSD retrievals? For that we need a different radar. Let's try S-PolKa during the NAME project in 2004.
In [26]:
fdir = '/Users/tjlang/Documents/OVWST/NAME/output/20040805/'
files = sorted(glob.glob(fdir+'*nc'))
radarS = pyart.io.read(files[0])
print(radarS.fields.keys())
In [27]:
dzS = radarS.fields['DZ']['data']
drS = radarS.fields['DR']['data']
kdS = radarS.fields['KD']['data']
With S-band, there are two different methodologies available: Bringi et al. (2013) and Bringi et al. (2004). The latter allows mu to vary, so let's try that.
In [28]:
d0, Nw, mu = csu_dsd.calc_dsd(dz=dzS, zdr=drS, kdp=kdS, band='S', method='2004')
radarS = add_field_to_radar_object(d0, radarS, field_name='D0', units='mm',
long_name='Median Volume Diameter',
standard_name='Median Volume Diameter', dz_field='DZ')
logNw = np.log10(Nw)
radarS = add_field_to_radar_object(logNw, radarS, field_name='NW', units='',
long_name='Normalized Intercept Parameter',
standard_name='Normalized Intercept Parameter',
dz_field='DZ')
radarS = add_field_to_radar_object(mu, radarS, field_name='MU', units='', long_name='Mu',
standard_name='Mu', dz_field='DZ')
In [29]:
limS = [-100, 100]
#Changing the color scale a bit to allow for the bigger drops in this case.
two_panel_plot(radarS, sweep=0, var1='D0', vmin1=0, vmax1=3.5, cmap1='GnBu', units1='mm',
var2='NW', vmin2=0, vmax2=8, cmap2='cubehelix', units2='log10(Nw)', xlim=limS,
ylim=limS)
In [30]:
#Now let's see what mu is up to ...
two_panel_plot(radarS, sweep=0, var1='DZ', vmin1=0, vmax1=65, cmap1='RdYlBu_r', units1='dBZ',
var2='MU', vmin2=0, vmax2=5, cmap2='cubehelix', units2='mu', xlim=limS,
ylim=limS)
So for the most part mu gets set to 3, except there is some variability in the cores and a few other places. Normally, mu is very difficult to retrieve accurately, but at least the user has the option of testing different methodologies with csu_dsd.
In [31]:
#And you can always ...
help(csu_dsd)
This module supplies a simple way of estimating specific differential phase (Kdp) via a methodology developed in the CSU Department of Electrical Engineering, and then subsequently adapted and used in the Department of Atmospheric Science. Filtering of the input differential phase field is based on a finite impulse response (FIR) filter, applied in a moving 3-km window. Then, an adaptive linear fit is applied iteratively to the filtered phase field, where half the slope of the linear fit at any specific gate is the Kdp estimate. The length of the line needed (i.e., number of gates considered) depends on the reflectivity at the gate in question.
Standard deviation of differential phase is estimated and used to remove noisy/bad data. Gaps are filled in the filter window if the holes are less than 20% of the window length. To use this module you need a gate spacing that divdes evenly into 3 km (e.g., 50, 100, 150, 250, 300 m, etc.). Let's get started with a NEXRAD example.
In [32]:
fdirN = '/Users/tjlang/Documents/OVWST/NEXRAD/'
filesN = sorted(glob.glob(fdirN+'KAMX*'))
radarN = pyart.io.read(filesN[0])
print(radarN.fields.keys())
In [33]:
def extract_unmasked_data(radar, field, bad=-32768):
"""Simplify getting unmasked radar fields from Py-ART"""
return radar.fields[field]['data'].filled(fill_value=bad)
In [34]:
# The module is a little slow on full volumes,
# so for this tutorial we'll only examine a single sweep
radarN = radarN.extract_sweeps([0])
dzN = extract_unmasked_data(radarN, 'reflectivity')
dpN = extract_unmasked_data(radarN, 'differential_phase')
# Range needs to be supplied as a variable, and it needs to be the same shape as dzN, etc.
rng2d, az2d = np.meshgrid(radarN.range['data'], radarN.azimuth['data'])
OK, now we have all our needed inputs to calculate Kdp. The function we will call is csu_kdp.calc_kdp_bringi(), and it returns Kdp, filtered differential phase, and standard deviation of differential phase, in that order. Input variables can be 1D (rays) or 2D (azimuth/elevation and rays). The fundamental algorithm works on a ray-by-ray basis.
In [35]:
kdN, fdN, sdN = csu_kdp.calc_kdp_bringi(dp=dpN, dz=dzN, rng=rng2d/1000.0, thsd=12, gs=250.0, window=5)
The gs keyword is the gate spacing that the radar data use, in meters. Window is the length (km) over with to apply the phase filtering. Default is 3 km. The thsd keyword is the threshold to apply to the standard deviation of differential phase (deg). OK, let's see what happened.
In [36]:
radarN = add_field_to_radar_object(kdN, radarN, field_name='KDP', units='deg/km',
long_name='Specific Differential Phase',
standard_name='Specific Differential Phase',
dz_field='reflectivity')
radarN = add_field_to_radar_object(fdN, radarN, field_name='FDP', units='deg',
long_name='Filtered Differential Phase',
standard_name='Filtered Differential Phase',
dz_field='reflectivity')
radarN = add_field_to_radar_object(sdN, radarN, field_name='SDP', units='deg',
long_name='Standard Deviation of Differential Phase',
standard_name='Standard Deviation of Differential Phase',
dz_field='reflectivity')
In [37]:
#First let's see what the original data looked like.
limN = [0, 150]
two_panel_plot(radarN, sweep=0, var1='reflectivity', vmin1=0, vmax1=65.0,
cmap1='RdYlBu_r', units1='dBZ',
var2='differential_phase', vmin2=0, vmax2=180,
cmap2='cubehelix', units2='deg',
xlim=limN, ylim=limN)
In [38]:
#Now let's see the filtered and specific differential phase fields
two_panel_plot(radarN, sweep=0, var1='FDP', vmin1=0, vmax1=180,
cmap1='cubehelix', units1='deg',
var2='KDP', vmin2=-5, vmax2=5, cmap2='RdYlBu', units2='deg/km',
xlim=limN, ylim=limN)
If window is too short, the resultant Kdp field is a little checkerboard-ish. That is a known characteristic of this algorithm. Due to the along-ray alignment of the storm, there is a lot of phase shift and possibly even some backscatter effects, and that is also a contributor. The best and most reliable results are obtained in the heavier cores, as expected.
One very nice aspect of this algorithm is it returns the standard deviation of differential phase. You may notice the removal of substantial clear-air echo in the processed phase fields. This is the standard deviation thresholding (thsd keyword) in effect. More on that in our next tutorial section.
In [39]:
#but first ...
help(csu_kdp.calc_kdp_bringi)
This module is a hodgepodge of various algorithms, most of them focused on quality control. They are a series of functions that return masks. Apply the masks to your radar data to clean them up. Let's see how this works in practice, continuing from the previous section.
In [40]:
drN = extract_unmasked_data(radarN, 'differential_reflectivity')
insect_mask = csu_misc.insect_filter(dzN, drN)
sdp_mask = csu_misc.differential_phase_filter(sdN, thresh_sdp=13)
In [41]:
bad = -32768
dz_insect = 1.0 * dzN
dz_insect[insect_mask] = bad
dz_sdp = 1.0 * dzN
dz_sdp[sdp_mask] = bad
In [42]:
radarN = add_field_to_radar_object(dz_insect, radarN, field_name='DZ_insect', units='dBZ',
long_name='Reflectivity (Insect Filtered)',
standard_name='Reflectivity (Insect Filtered)',
dz_field='reflectivity')
radarN = add_field_to_radar_object(dz_sdp, radarN, field_name='DZ_sdp', units='dBZ',
long_name='Reflectivity (Phase Filtered)',
standard_name='Reflectivity (Phase Filtered)',
dz_field='reflectivity')
In [43]:
limN = [-150, 150]
two_panel_plot(radarN, sweep=0, var1='DZ_insect', vmin1=0, vmax1=65,
cmap1='RdYlBu_r', units1='dBZ',
var2='DZ_sdp', vmin2=0, vmax2=65, cmap2='RdYlBu_r', units2='dBZ',
xlim=limN, ylim=limN)
Of course, the real power will come by combining multiple filters together. And then we can despeckle to remove most remaining spurious echo.
In [44]:
new_mask = np.logical_or(insect_mask, sdp_mask)
dz_qc = 1.0 * dzN
dz_qc[new_mask] = bad
dz_pre_despeck = 1.0 * dz_qc
mask_ds = csu_misc.despeckle(dz_qc, ngates=4)
dz_qc[mask_ds] = bad
radarN = add_field_to_radar_object(dz_qc, radarN, field_name='DZ_qc', units='dBZ',
long_name='Reflectivity (Combo Filtered)',
standard_name='Reflectivity (Combo Filtered)',
dz_field='reflectivity')
radarN = add_field_to_radar_object(dz_pre_despeck, radarN,
field_name='DZ_before', units='dBZ',
long_name='Reflectivity (No Despeckling)',
standard_name='Reflectivity (No Despeckling)',
dz_field='reflectivity')
In [45]:
limN = [-150, 150]
two_panel_plot(radarN, sweep=0, var1='DZ_qc', vmin1=0, vmax1=65,
cmap1='RdYlBu_r', units1='dBZ',
var2='SDP', vmin2=0, vmax2=25, cmap2='YlOrRd', units2='deg',
xlim=limN, ylim=limN)
The way csu_misc.insect_filter works is it looks in various reflectivity ranges. If differential reflectivity exceeds a given threshold in that range, the data are flagged as likely insects (or other biological scatterers). The defaults are some values that were tuned for the NAME project.
In [46]:
print(csu_misc.DEFAULT_DZ_RANGE)
print(csu_misc.DEFAULT_DR_THRESH)
So below 10 dBZ, if Zdr exceeds 1 dB the data are flagged, for Zh = 10-15 dBZ Zdr above 1.3 dB is flagged, etc. The user can make their own thresholds by adjusting the dz_range and dr_thresh keywords in csu_misc.insect_filter(). The dz_range keyword expects a list of 2-element tuples (lower and upper bounds for reflectivity) and dr_thresh expects a list of scalars (Zdr threshold value).
In [47]:
#Here is more info
help(csu_misc.insect_filter)
Let's see how the despeckle filter works by examining before and after behavior.
In [48]:
limN = [-150, 150]
two_panel_plot(radarN, sweep=0, var1='DZ_before', vmin1=0, vmax1=65,
cmap1='RdYlBu_r', units1='dBZ',
var2='DZ_qc', vmin2=0, vmax2=65, cmap2='RdYlBu_r', units2='dBZ',
xlim=limN, ylim=limN)
So, it gets rid of short, standalone "specks" of echo. The user can change how many gates they want to count as a speck via the ngates keyword. Feel free to vary the despeckling length and SDP threshold to suit your desire for elimination of non-meteorological echo. We could also play with the Zh/Zdr thresholds in the insect filter. But by and large, for the settings we used, the worst of the non-meteorological echo is removed. That is good enough for a demonstration!
In [49]:
#More info
help(csu_misc.despeckle)
There are a few other basic filters, but you can find all that out by ...
In [50]:
help(csu_misc)
In [ ]: