filtering around aliased regions

[joshua-wx]

Hi meteocat,

I'm looking at implementing your excellent dual-PRF filter library for processing Australian weather radar data, but I've ran into an issues for aliased regions. From my understanding, cmean_sc and cmean techniques should not be removing the edges of aliased region? Here's an example showing the issue. I'm using the same implementation of cmean and cmean_sc as the notebook example:

vcor.correct_dualprf(radar=radar, two_step=True, radar_ffn=radar_ffn,
                 method_det='cmean_sc', kernel_det=np.ones((7, 7)),
                 method_cor='median', kernel_cor=np.ones((7, 7)),
                 vel_field=vel_field, new_field='vcor_cmean_sc')   
vcor.correct_dualprf(radar=radar, two_step=True, radar_ffn=radar_ffn,
                 method_det='cmean', kernel_det=np.ones((7, 7)),
                 method_cor='cmean', kernel_cor=np.ones((7, 7)),
                 vel_field=vel_field, new_field='vcor_cmean')

I've had to create a custom reader for our Australian metadata in the odimh5 format to create the correct information fields about high/low PRF. See the function _get_prf_pars_odimh5 in my fork

Here's the odim file from the figure above: https://www.dropbox.com/s/e5egjltch8rl7nk/54_20141014_105402.pvol.h5?dl=0

Any help would be greatly appreciated! Thanks! Joshua

[paltube]

Hi Joshua, You're right, it certainly shouldn't be removing those edges. I'll have a look at the data and process asap and come back to you. Thank you for using the code! It hasn't been tested on many different data and this will definitely help improve :)

[joshua-wx]

Great, thank you! The Australian odimh5 data is a little different. pyart can read it fine using pyart.aux_io.read_odim_h5(name, file_field_names=True), but it doesn't read the dualPRF metadata, that's why I had to read it directly using h5py. rapic_UNFOLDING provides the unfolding ratio, rapic_HIPRF indicates whether the high prf rays are even or odd rays, and there's also a highprf field.

[joshua-wx]

Hi paltube,

Just checking if you got a chance to look into this. We are really keen to use your work for our operational weather radars!

Cheers, Joshua

[paltube]

Hi @joshua-wx,

Please, excuse me for the delay. I'm currently reviewing the issue and I hope to be able to find a clue in the next hours/day. Thanks a lot for the patience!

[paltube]

I've done some testing and the problem seems comes from an algorithm-error in the current version of the code that does not happen in the very early versions that we have. This is good news since it means I should be able to trace back the error, but its gonna take me some more time. I'm gonna work mostly on it during the next days, so I'll keep you updated. Cheers!

[joshua-wx]

That's fanatic news! Thanks very much for the update and sorry to bother you about it. Let me know when I can start testing new versions :+1:

[paltube]

Hi Joshua,

I finally found the problem and fixed it (I hope :P). It was basically a problem with the 'fold_circular' function. The new version is currently only in my fork; I've tested it but I think it is safer if you test it too before merging with the origin in meteocat.

I've included your function for retrieving the prf settings of the odim data. I've changed it a bit for its output to match the 'instrument_parameters' dictionary that we get by default when applying the arm-pyart function 'pyart.io.read' in our IRIS data. Now, you can call this function (_instrument_parametersodim5) right after reading the data, and then apply the correction(s):

# Load volume data
radar = pyart.aux_io.read_odim_h5(file_in, file_field_names=True)

# Add instrument parameters
radar = vcor.instrument_parameters_odim5(radar=radar, odim_file=file_in)

# CORRECT CMEAN
vcor.correct_dualprf(radar=radar, two_step=True,
                      method_det='cmean', kernel_det=np.ones((7, 11)),
                      method_cor='median', kernel_cor=np.ones((7, 11)),
                      vel_field='VRADH', new_field='vcor_cmean')

 # CORRECT CMEAN_SC
 vcor.correct_dualprf(radar=radar, two_step=True,
                      method_det='cmean_sc', kernel_det=np.ones((7, 11)),
                      method_cor='median', kernel_cor=np.ones((7, 11)),
                     vel_field='VRADH', new_field='vcor_cmean_sc')

[joshua-wx]

Hi Patricia,

Thanks for working on this fix! I've tested your fork for two other cases and the issues seems to be resolved (while still correcting dual PRF artifacts). The odimh5 reader is very useful too :+1:

Which technique do you recommend that balances processing time with correction skill? We are planning to process several decades of Doppler data using this technique so we need something rather efficient. I found this configuration is reasonably fast:

# CORRECT CMEAN
vcor.correct_dualprf(radar=radar, two_step=True,
                      method_det='cmean', kernel_det=np.ones((7, 11)),
                      method_cor='cmean', kernel_cor=np.ones((7, 11)),
                      vel_field='VRADH', new_field='vcor_cmean')

Cheers, Joshua

[paltube]

Hi Joshua,

It's great news that it works! I'll pull the new version to the original repository.

About your question, I think that what you propose may be the best option. Some considerations:

• Due to the aliasing, you should use either 'cmean' or 'cmean_sc' in order to keep the edges of the aliased regions. These two methods perform similarly but perhaps the first is slightly more robust (I haven't quantified this, it's just my feeling after testing it with our data for 2 years).

• Currently, all the methods that involve 'median' calculation are a bit slow, due to the implemented median filter routine. It used to be the only one available that accepted Nan values; I will review the libraries to see if currently there's any other that can be applied, but it's unlikely.

• Looking at your data while testing the fix I saw that there were some large clusters of non-matching data (image below). I think these are not aliased regions but seem to be regions of strong azimuthal propagation of the dual-PRF outliers, corresponding to regions of high W and low velocity (probably low SQI). We don't get those artifacts in our data, which has much coarser resolution, so I'm not sure about this but, if I'm right, these regions are difficult to correct. So, the proposal is to use an asymmetrical filtering window, wider in range than azimuth, like the one in the example ((7, 11)). This does not correct the regions completely but improves the results. Another possibility that I can think of is to apply iteratively the 'cmean' method to the data, but I'd say the improvement won't be worth the extra computing time (I've tried it and it requires 4-5 iterations to correct them, and still not completely).

[joshua-wx]

Thanks for the suggestions for the optimal technique and filter window. I've now got those implemented and it seems to be running fine. I suspect those regions to the SW of the Sydney radar might have low SQI due to interference from clutter - this is a common problem inland of this radar.

Thanks again :smile:

[joshua-wx]

One more issue relating to this has appeared today! I tested the your script for an odimh5 volume that had fixed PRF and it got stuck in a couple places when building the instruments_parameter field (which I then use to test whether to skip running vcor). The code block below is for the _get_prf_pars_odimh5 function, I've changed "if prf_ratio != 'None':" to be "if prf_ratio != b'None':" (binary string) and moved "prf_type = d_how['rapic_HIPRF']" after the test for dual PRF metadata, as this field only exists when the radar operates in dual PRF mode.

I should also point out the odimh5 function is specific to the Australian flavour of odimh5. I don't think it will work with other countries odimh5 as the rapic_##### fields are specific to us.

Thanks! Joshua

def _get_prf_pars_odimh5(odim_file, nrays, nsweeps, sw_start_end):
    """
    Retrieves PRF scanning parameters from odim5 file, shaped for 
    building the 'instrument_parameters' dictionary: 
    nyquist velocity, PRF, dual PRF factor and PRF flags for each ray
    (if batch mode dual-PRF).

    Parameters
    ----------
    radar : Radar
        Py-ART radar structure

    Returns
    -------
    ny_array : numpy array (float)
        Nyquist velocity for each ray.
    prt_array : numpy array (float)
        PRT for each ray.
    prt_mode_array: numpy array (str)
        PRT mode for each sweep, 'dual' or 'fixed'.
    prt_ratio_array: numpy array (float)
        PRT ratio for each ray.
    prf_flag_array: numpy array (bool int)
        Ray PRF flag, high (0) or low (1) PRF.    
    """

    ny_array = np.zeros(nrays)
    prt_array = np.zeros(nrays)
    prt_mode_array = np.repeat(b'single', nsweeps)
    prt_ratio_array = np.ones(nrays)
    prf_flag_array = np.zeros(nrays)

    with h5py.File(odim_file, 'r') as hfile:

        for sw in range(0, nsweeps):

            # extract PRF/NI data from odimh5 file
            d_name = 'dataset' + str(sw+1)
            d_how = hfile[d_name]['how'].attrs
            ny = d_how['NI']              # Nyquist
            prf_h = d_how['highprf']
            prf_ratio = d_how['rapic_UNFOLDING'] # the prf ratio (e.g. 2:3 etc) or None

            # extract rays for current sweep
            ray_s, ray_e = sw_start_end(sw) # start and end rays of sweep
            ray_e += 1

            # Assign values
            prt_array[ray_s:ray_e] = 1/prf_h
            ny_array[ray_s:ray_e] = ny

            if prf_ratio != b'None':

                prf_type = d_how['rapic_HIPRF']

                prt_mode_array[sw] = b'dual'

                fact_h = float(prf_ratio.decode('ascii')[0])
                fact_l = float(prf_ratio.decode('ascii')[2])

                prt_ratio_array[ray_s:ray_e] = fact_l/fact_h
                flag_sw = prf_flag_array[ray_s:ray_e]

                if prf_type==b'EVENS':
                    flag_sw[::2] = 1 # with 1 as the first index, 1=1, 2=0, 3=1
                elif prf_type==b'ODDS': #odds=0, evens=1
                    flag_sw[1::2] = 1 #with 1 as the first index, 1=0, 2=1, 3=0
                else:
                    print('error, unknown flag type', prf_type)
            else:
                prt_mode_array[sw] = b'fixed'

    return ny_array, prt_array, prt_mode_array, prt_ratio_array, prf_flag_array.astype(int)