Over aggressive filtration

Corrections/Data Quality
1

Hi All,

I want your help with noise removal and dealiasing. I am conducting quality checks and filtration on my data. However, the process also removes good data. I attempted to adjust the values based on the characteristics of the field, but it has yet to be effective. The current method of keeping the VT above a certain threshold still leaves some noise present, while if the threshold is set too low, the filtration becomes too aggressive. I would greatly appreciate any advice or guidance you can give me. Please look at the plot; if you think “it is a tiny patch”, it is possibly a tornado, and I want to perform some dual doppler analysis, and neither can’t afford noise nor lose even that much of the data.
Thank you!

download
download873Ă—717 246 KB 

import pyart
import numpy as np
from pyart.retrieve import get_freq_band
from matplotlib import pyplot as plt
import glob, os

def align_field(field):
    values, counts = np.unique(field['data'], return_counts=True)
    c_value = values[np.argmax(counts)]
    field['data'] = np.array([c_value])
    return field

def align_radar_coords(radar):
    radar.longitude = align_field(radar.longitude)
    radar.latitude = align_field(radar.latitude)
    radar.altitude = align_field(radar.altitude)
    radar.altitude_agl = align_field(radar.altitude_agl)
    return radar

def filter_radar(radar, vel_field="VEL_F", refl_field="DBZHCC_F", 
                 ncp_field="NCP", rhv_field="RHOHV", phi_field="PHIDP"):
    '''Remove noise based on velocity texture and mask all the fields'''
    ##Drop some fields
    fields_to_drop = ["DBMHC", "DBMVC", "VS", "VS_F", "VL", "VL_F", "DBZHCC"]
    for field in fields_to_drop:
        if field in radar.fields:
            del radar.fields[field]

    texture = pyart.retrieve.calculate_velocity_texture(radar, 
                                                        vel_field=vel_field, 
                                                        wind_size=5, 
                                                        check_nyq_uniform=False)
    radar.add_field('VT',texture,replace_existing=True)

    ##create gatefilter
    gf = pyart.filters.GateFilter(radar)
    BAND = get_freq_band(radar.instrument_parameters['frequency']['data'])
    if BAND == "X":
        vt = 5
    if BAND == "C":
        vt = 25
    gf.exclude_above('VT', vt)
    gf.exclude_outside(refl_field, -30, 100)
    gf.exclude_below("SNRHC", -3)
    gf_despeckeld = pyart.correct.despeckle_field(radar, refl_field, gatefilter=gf)
    corr_ZH = radar.fields[refl_field].copy()
    ZH_array = np.ma.masked_where(gf_despeckeld.gate_included == False, radar.fields[refl_field]['data'])
    corr_ZH['data'] = ZH_array
    radar.add_field('DBZH',corr_ZH,replace_existing=True)
  
    ##get the mask
    mask = np.ma.getmask(radar.fields['DBZH']['data'])
    radar.scan_type = b'ppi'
    ##iterate through remaining fields
    skip_fields = ["DBZHC", "VEL"]
    for field in radar.fields.keys():
        if any(field == skip_field for skip_field in skip_fields):
            continue
        ##mask the field
        radar.fields[field]['data'] = np.ma.masked_where(mask, radar.fields[field]['data'])
    return radar

def dealiase(radar, vel_name = "VEL_F"):
    try: 
        gatefilter = pyart.correct.GateFilter(radar)
        corr_vel   = pyart.correct.dealias_region_based(
            radar, vel_field=vel_name, keep_original=False, gatefilter = gatefilter)
        radar.add_field(vel_name, corr_vel, replace_existing=True)
    except:
        None

def qc(file=None, savedir = "COW1/QC/"):
    radar = pyart.io.read(file)
    radar = align_radar_coords(radar)
    print("")
    print(radar.fields.keys())
    print("Coords: ", radar.longitude['data'],
          radar.latitude['data'],
          radar.altitude['data'],)
    radar = filter_radar(radar)
    dealiase(radar)
     ##pyart.io.write(os.path.join(basedir, savedir, file))
    return radar

def plot_radar(radar):
    # Let us view the vleocity with the filter applied.
    fig = plt.figure(figsize=[10, 8])
    display = pyart.graph.RadarMapDisplay(radar)
    ax = plt.subplot(221)
    display.plot_ppi("DBZHC", vmin=-10., vmax=70., cmap='pyart_HomeyerRainbow')
    ax = plt.subplot(222)
    display.plot_ppi("DBZH", vmin=-10., vmax=70., cmap='pyart_HomeyerRainbow', title='filtered')
    ax = plt.subplot(223)
    display.plot_ppi("VEL", sweep=0, cmap="pyart_NWSVel", vmin = -30, vmax = 30,)
    ax = plt.subplot(224)
    display.plot_ppi("VEL_F", sweep=0, cmap="pyart_NWSVel", vmin = -30, vmax = 30, title = "filtered")
    plt.show()

basedir = "/depot/dawson29/data/Projects/PERiLS/obsdata/2022/Illinois_Mobile_Radar/IOP1/"
cow = os.path.join(basedir, "COW1/cmerged/20220322/")

files = glob.glob(os.path.join(cow, "*"))
files.sort()
len(files)

radar = qc(files[112])
plot_radar(radar)

In case you need a data file, here is a link.

2

Yes we have had challenges with this as well. The TVS is creating texture… Have you tried adding an include_inside to force Z > some threshold back in
unless the TVS is right on the edge of a BWER or suchlike it should have a decent Z signature as well…

Or you can use a fuzzy logic approach… This code really needs some refresh but this is how we do it on ARM data cmac2.0/cmac_processing.py at main · EVS-ATMOS/cmac2.0 · GitHub

2 Likes
3

In addition to what Scott suggests after putting the reflectivity greater than ~10 dBZ back in the field I would apply despeckling to get rid of the isolated speckles that can create errors in dealiasing. This should not affect your tornado since it’s embedded in an MCS.

1 Like
4

Thank you, @scollis and @rcjackson. Your recommendations worked; However, I’d like to mention that I had to fine-tune the wind size.

5

Great to know! Would be great if you (eventually) share some code… What would be amazing is a cookbook on this cc @mgrover1

1 Like