over aggressive filtration

[syedhamidali]

Dear PyArt Team,

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 basically a tornado, and I want to perform some dual doppler analysis, and neither can't afford noise nor lose even that much amount of the data.

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)

[scollis]

Hey! Can you post this at https://openradar.discourse.group/ ?
This is a great radar discussion.. BTW the answer is the NCP is low because the spectral width is huge in those gates.. Suggest using texture of velocity. Once you move to the Discorse @rcjackson and @zssherman can share some CMAC code.

[syedhamidali]

Thank you, Dr. @scollis! I am using the velocity texture only. Well, let me post it to openradar discourse group.