Advection Correction

[syedhamidali]

Hi,

I'm working on interpolating gridded NEXRAD radar data with advection correction but encountering unexpected wobbling in the interpolated results (not the edge effects), as shown in the attached GIF. I've prepared a notebook to help diagnose the issue, which you can download and run directly from the attached link. Any assistance would be greatly appreciated.

Thank you! Hamid

https://github.com/syedhamidali/dda_workflow/blob/main/Radar_Interpolation_Adv_Corr.ipynb

https://github.com/pySTEPS/pysteps/assets/35923822/c0426845-af9f-48d7-ba44-2cc13adc8598

[dnerini]

hi @syedhamidali, thanks for reporting this interesting case! Can you be a bit more specific in what you're trying to do, what data you're using (in particular their spatial and temporal resolution) and the expected outcome?

In general, the quality of the advection correction is directly linked to the quality of the estimated flow between successive frames. Some of the artifacts might therefore be a consequence of errors in the optical flow estimation. Have you already tried tuning the parameters of the optical flow method?

[syedhamidali]

Hi @dnerini,

I’m working on obtaining regularly spaced (1-min or 2-min intervals) data by interpolating NEXRAD radar data. First, I use Py-ART to grid the data into 3d volumes. Then, I perform advection correction at the 500m altitude of the dataset, as demonstrated below. Also, I have attached a link to the notebook, where I have applied it . I tried it on different resolutions, for example, here, I am doing it on 500 m resolution. I have attached the two animations where the first one is advection corrected, and the other one is simply linearly interpolated to 1min. The second one, you can see the smoothness in the flow, but in the first one, it looks it is wobbling a little. I think it advecting more than it should and that is why it is wobbling, but I am not sure.

Advection Correction

Expected (Simply resampled to 1min)

def advection_correction_ds(radar_ds, tintv_obs, tintv,
                            base_field_name='reflectivity_masked', 
                            use_base_field=False, method="LK"):
    # Evaluate advection
    oflow_method = motion.get_method(method)
    fd_kwargs = {"buffer_mask": 10}  # avoid edge effects

    # Decide whether to use a single optical flow field for all fields in the Dataset
    # (using the "base" field) or to use separate ones for each field
    base_field = radar_ds[base_field_name]
    if use_base_field:
        oflow_field = oflow_method(base_field, fd_kwargs=fd_kwargs)
        oflow_field_list = [oflow_field] * len(radar_ds.items())
    else:
        oflow_field_list = []
        for field_name, field_da in radar_ds.items():
            oflow_field = oflow_method(field_da, fd_kwargs=fd_kwargs)
            oflow_field_list.append(oflow_field)

    # Perform temporal interpolation on all variables in Dataset using the flow field derived from either 
    # the "base" field (by default, reflectivity), or each field individually

    tbgn = base_field[0].coords['time_seconds'].values.item()   # Need to get scalar value, not 0-d
                                                                # numpy array
    print(tbgn)
    print(tintv)
    radar_ds_list = []
    x, y = np.meshgrid(
        np.arange(base_field[0].shape[1], dtype=float), np.arange(base_field[0].shape[0], dtype=float),
    )

    new_time = tbgn
    for i in np.arange(tintv, tintv_obs + tintv, tintv):
        new_time = new_time + tintv
        field_interp_list = []
        for field_da, oflow_field in zip(radar_ds.values(), oflow_field_list):
            pos1 = (y - i / tintv_obs * oflow_field[1], 
                    x - i / tintv_obs * oflow_field[0])
            pos2 = (y + (tintv_obs - i) / tintv_obs * oflow_field[1], 
                    x + (tintv_obs - i) / tintv_obs * oflow_field[0])

            fieldt1 = map_coordinates(field_da[0], pos1, order=1)
            fieldt2 = map_coordinates(field_da[1], pos2, order=1)

            field_interp = field_da.isel(time=[0]).copy()
            field_interp[:] = ((tintv_obs - i) * fieldt1 + i * fieldt2) / tintv_obs
            try:
                field_interp.coords['time'] = field_interp['time'] + np.timedelta64(int(new_time - tbgn), 's')
            except TypeError:
                field_interp.coords['time'] = field_interp['time'] + timedelta(seconds=int(new_time - tbgn))
            field_interp.coords['time_seconds'] = new_time
            field_interp_list.append(field_interp)

        radar_ds_interp = xr.merge(field_interp_list)
        radar_ds_list.append(radar_ds_interp)

    return radar_ds_list

def advection_correction(arr, tintv_obs, tintv):
    """
    R = np.array([qpe_previous, qpe_current])
    T = time between two observations (5 min)
    t = interpolation timestep (1 min)
    """

    # Evaluate advection
    oflow_method = motion.get_method("LK")
    fd_kwargs = {"buffer_mask": 10}  # avoid edge effects
    V = oflow_method(arr, fd_kwargs=fd_kwargs)

    # Perform temporal interpolation
    # arr_d = np.zeros((arr[0].shape))
    arr_list = []
    x, y = np.meshgrid(
        np.arange(arr[0].shape[1], dtype=float), np.arange(arr[0].shape[0], dtype=float),
    )
    for i in np.arange(tintv, tintv_obs + tintv, tintv):

        pos1 = (y - i / tintv_obs * V[1], x - i / tintv_obs * V[0])
        R1 = map_coordinates(arr[0], pos1, order=1)

        pos2 = (y + (tintv_obs - i) / tintv_obs * V[1], x + (tintv_obs - i) / tintv_obs * V[0])
        R2 = map_coordinates(arr[1], pos2, order=1)

        arr_interp = ((tintv_obs - i) * R1 + i * R2) / tintv_obs
        arr_list.append(arr_interp)

    return arr_list

[dnerini]

Hi @syedhamidali thanks for the additional details. I don't see any obvious problem with your code. The output seems reasonably good to me. There are as you say some artifacts in your results that I suspect are due to edge effects during the optical flow estimation. You could try to mitigate them by increasing the buffer_mask argument. You could also try to threshold the reflectivity field with some higher values when computing the optical flow, to reduce the sensitivity to noisy, low intensity echoes.

[syedhamidali]

Hi @dnerini thank you for your comments, will try to implement these.