Csu radartools tutorial
From Lrose Wiki
Contents
Overview
This tutorial will walk through the steps to read and convert a raw radar file, quality control the data (removing unwanted artifacts and echoes, unfold radial velocity and calculate Kdp), then grid the data, and finally process the gridded data to end up with hydrometeor identification, precipitation rates, and DSD information. This is generally focused on surface polarimetric X, C, or S-band radar. We will apply a series of corrections
- Unfold radial velocity
- Calculate Kdp
- Apply thresholds to remove non-meteorological data
- Calculate the attenuation and differential attenuation
- Despeckle and remove 2nd trip
Modules available from GitHub
Scientific Background
- Kdp calculation
Hubbert, J., and V. N. Bringi. "An iterative filtering technique for the analysis of copolar differential phase and dual-frequency radar measurements." Journal of Atmospheric and Oceanic Technology 12.3 (1995): 643-648.
Step-by-step Instructions
The Example file is from the C-band polarimetric SEA-POL radar. The raw data are in sigmet native format.
Python Imports
import numpy as np from copy import deepcopy import os import pyart from CSU_RadarTools.csu_radartools import csu_kdp from CSU_RadarTools.csu_radartools import csu_misc
Convert to CFRADIAL
Convert the data to cfradial data to be more easily transferable to other programs.
We can do this using RadxConvert from a command line, but we can also do it from python with the os command.
mydir = '/Users/bdolan/scratch/LROSE/SEA-POL/ppi/' os.system(f'RadxConvert -f {mydir}SEA20190917_071005 -outdir {mydir}cfrad')
This will put the output in a directory for the date in directory 'cfrad'.
Now we have a file called:
cfrad.20190917_071006.545_to_20190917_071353.149_SEAPOL_SUR.nc
We can read this into radar with pyart:
radar = pyart.io.read(f'{mydir}cfrad/cfrad.20190917_071006.545_to_20190917_071353.149_SEAPOL_SUR.nc')
There are a lot of variables in this file, but the ones of interest are:
radar.fields.keys()
DBZ: Reflectivity
ZDR: Differential Reflectivity
UNKNOWN_ID_82: Correlation Coefficient
VEL: Radial velocity
SQI: Signal Quality Index
PHIDP: Differential phase
SNR: Signal to noise ratio
QC with CSU_Radartools and PyART
Do some quality control using PyART and CSU_Radartools.
Here is the original reflectivity field:
Special handling for SEA-POL
These are PECULIARS to the SEA-POL data. YMMV
- Blanked sector over the wheelhouse
We will look for large blocks of azimuthal jumps to mask them. Get the difference between each azimuth and look for jumps larger than 30º. If so, mask the data so it plots correctly.
az_diff = np.diff(radar.azimuth['data']) jumps = np.where(np.abs(az_diff) >= 30.0)[0] if len(jumps): for f in radar.fields.keys(): for j in jumps: radar.fields[f]['data'][j].mask = True radar.fields[f]['data'][j+1].mask = True
The new reflectivity field is now masked over the wheelhouse:
- The correlation coefficient is incorrect, but we can calculate it from UNKNONW_ID_82 field.
Now correct the UNKNOWN_ID_82 field to transform into the correlation coefficient (Again, this is only needed for this data).
ccorr = deepcopy(radar.fields['UNKNOWN_ID_82']['data']) ccorr[ccorr<0] += 65535 ccorr = (ccorr-1)/65533.0 radar.add_field_like('RHOHV','CC',ccorr)
Calculate Kdp
This method uses the Hubbert and Bringi FIR filter method. We note that the window / gate spacing must be an even number or it will fail. it will also calculate the standard deviation of the phase which can be used to remove non-meteorological data.
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)
dz_um = extract_unmasked_data(radar,'DBZ') dp_um = extract_unmasked_data(radar, 'PHIDP')
Range needs to be supplied as a variable, and it needs to be the same shape as dzN, etc.
rng2d, az2d = np.meshgrid(radar.range['data'], radar.azimuth['data']) kdN, fdN, sdN = csu_kdp.calc_kdp_bringi(dp=dp_um, dz=dz_um, rng=rng2d/1000.0, thsd=12, gs=250.0, window=7, nfilter=2)
kdN is the specific differential phase, fdN is the filtered differential phase, and sdN is the standard deviation of the phase
Add this back into the radar structure.
field_dict = {'data': kdN, 'units': 'deg/km', 'long_name': 'specific_differential_phase', 'standard_name': 'Kdp', '_FillValue': '-32768'} radar.add_field('Kdp', field_dict, replace_existing=True)
Thresholding
Thresholding
Sometimes the sdp goes nuts in high reflectivity regions, and pair this mask with relatively low dBZ regions. Here I am using a threshold of 10, but I usually do a sensitivity study to figure out a good value, typically ranging from 10 to 25. However, remember that things like hail can have high sdN.
sdp_mask = csu_misc.differential_phase_filter(sdN, thresh_sdp=10) sdp_mask = np.logical_and(sdp_mask, radar.fields['DBZ']['data'] <= 10)
old_var_mask = deepcopy(radar.fields['DBZ']['data'].mask) new_var_mask = np.logical_or(old_var_mask, sdp_mask) radar.fields['DBZ']['data'].mask = new_var_mask
Well that didn't get us vary far (this is not unexpected given the RF interference with this particular dataset). But we have other tools at our disposal, and this is usually paired with another filter like rho_hv or SQI.