Elevation compensation

[kingjml]

Changes in aircraft elevation change the location of radar surface peak in the SnowRadar data. This is not a huge issue for deriving snow depth but it becomes an issue for visualizations where the aircraft is moving around. CRESIS has established code to allign the bins written in MATLAB (See here). This needs to ported and generlized in the snowradar class so it can be applied against incoming data.

[ajutila]

Hi there! Are there any updates on the elevation compensation? I'm slowly starting to try out the code with our data, and there seems to be a difference if I pick the interfaces using the whole data frame or only a subset of it (see attached figures). Or it might just be me being inexperienced in python.

I guess using the elevation compensation first might be the trick here, because I don't feel like determining a subset for each individual frame in our dataset...

[kingjml]

It should not be that there are differences depending on the length of processing so this is concerning. Each radar trace is processed on its own. Could you post an example notebook to reproduce the error?

[ajutila]

The only thing that I changed between those figures is the variable radar_sub, from a subset _radar_sub = np.transpose(radar_dat.data_radar)[3000:7000,0:radar_dat.dataradar.shape[0]] to the whole frame _radar_sub = np.transpose(radar_dat.dataradar)

I think I could make an example notebook, but currently I'm not using notebooks. I adapted your AWI demo notebook into a python script, which I run with a bash script on our supercomputer to process several data frames simultaneously.

[kingjml]

I'll try this when I'm back in the office next week. The interface bin picks should be completely independent of the reference frame (ie number of traces). Perhaps there is a smoothing function that is sensitive to the total number of bins evaluated but that would be a bug.

I don't believe that elevation compensation will fix this.

[m9brady]

I got the same result as Arttu after removing the y-axis subsetting... May be a bug with the picklayers function.

A temporary band-aid fix before any changes to picklayers are developed could be to automatically subset in the y-axis (Fast Time Bin?) using the following:

def get_bounds(row, thresh=0.3):
    '''retrieve the "valid data" index ranges of a given 1D array'''
    idxes = np.argwhere(np.abs(row) > thresh)
    return (np.min(idxes), np.max(idxes))

# transpose data array
arr = np.transpose(radar_dat.data_radar)
# for each column of arr, determine the beginning and end of "valid" data
result = np.apply_along_axis(get_bounds, 0, arr)
# from the result of above, retrieve the largest "min" index and the smallest "max" index
idxmin, idxmax = np.max(result[0,:]), np.min(result[1,:])
# subset to this "valid" data range
radar_sub = arr[idxmin:idxmax,:]

Feel free to experiment with values for thresh in get_bounds(). Here are the results when using thresh=0.3:

fig, axes = plt.subplots(ncols=2, figsize=(14,7))
ax_orig, ax_sub = axes
ax_orig.imshow(10*np.log10(arr), cmap='gray', vmin=0, vmax=80)
ax_orig.set_title('original data array')
ax_orig.hlines([idxmin, idxmax], 0, arr.shape[1], color='red', label='"valid" range')
ax_orig.legend(framealpha=1)
ax_sub.imshow(10*np.log10(radar_sub), cmap='gray', vmin=0, vmax=80)
ax_sub.set_title('subset array')
fig.tight_layout();

[ajutila]

Thanks for looking into this and confirming the result! I can process most of our data by limiting the fast time bin axis to ~3000:5000. However, our data has semi-regular ascents and descents between 200 ft and 500 ft to calibrate (/check instrument drift of) the EM-Bird, which is flown together with the snow radar, and that makes choosing a subset slightly funky business on those particular frames, e.g.

Our snow radar data is not yet fully pre-processed but waiting for CReSIS's help to go through coherent noise removal and deconvolution, so this is just me learning to work with the code and with python in general. Anyway, I hope the reason for the interface picking bug is found quickly and easily fixed!

[kingjml]

The first version of the elevation compensation had been added to the repo. The approach only works where the CRESIS data has already been supplied with a 'Surface' field. The radar traces are padded with 0s and adjusted to remove variations in flight altitude. Now the bins are referenced to the WGS-84 surface rather than to the flight altitude.

Some other minor changes: The matrix that holds the radar data is now transposed when its read from file so that the data products are equivalent between OIB and AWI. The snowradar sub-classes also now look for elevation corrections that are available in file and apply them (OIB ONLY). The AWI compensation method needs to be moved into the snowradar child class but for testing in the Dev branch this is enough to apply against your data and report.

However, this will not 'solve' the problem until we apply a standard window to limit the data above and below the surface. After compensation is applied for elevation we could do something simple like limit data to 5 m above and below the surface and then run the picker. If we don't do this the no-data to data transitions get picked as @m9brady has pointed out.

Summary: Elevation is the issue because no-data to data transitions are 0 to X causing the picker to think its a surface. Compensating for elevation will fix this so long as we apply a strategy to remove all the non-radar data regions.

[ajutila]

Hey, I tried now the newest version of the code where the elevation compensation is part of the SnowRadar class. I don't know if it's the best option to window the data using min/max surface_elev, like in the example notebook, because the AWI data has those ascents/descents to calibrate the EM-Bird every ~10 mins. The min/max values are so large that the non-radar data regions end up included anyway. (It might also be that the snow radar data during those ascents/descents is crap and have to be discarded anyway.) I tried a workaround using the mode of surface_elev and sort of a buffer instead of min/max values:

from scipy import stats depth_mode = stats.mode(np.around(radar_dat.surface_elev[:,0],2)).mode[0] depth_thresh = 10 valid_idx = np.argwhere((elev_axis<depth_mode+depth_thresh) & (elev_axis>depth_mode-depth_thresh)).flatten()

[kingjml]

Maybe we should add some windowing functions to the snowradar class.

I think its a good idea to look https://github.com/kingjml/pyWavelet/issues/2 at this point because filtering is needed to remove all data where the aircraft is not in level flight. This is generally not an issue for OIB but is a strong consideration for AWI. Pitch, Yaw, and Roll and be used to evaluate and flag the data. Then flagging can be used to ignore the non-level traces.

[kingjml]

Elevation compensation is fully functional in master with the latest commit. There are some larger issues with the AWI snow radar data itself that need to be resolved to answer some of the questions here. As for now elevation compensation has been added to pySnowRadar.