-
Notifications
You must be signed in to change notification settings - Fork 151
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Advection Correction #360
Comments
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? |
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. 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 |
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 |
Hi @dnerini thank you for your comments, will try to implement these. |
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
demo.mp4
The text was updated successfully, but these errors were encountered: