AWI bug in get_surface

[ajutila]

AWI data in its current state includes All-NaN traces and random traces where the last value has the highest value. That screws up the function get_surface, which in turn messes up get_bounds. Therefore the radar_subset has wrong boundaries, which allows radar data with NaN values to emerge, which in turn ends up ruining the interface picker.

    def get_surface(self, smooth=True):
        '''
        Simple surface tracker based on maximum
        This should be refined and is largely a place holder 
        TODO: surf_time is broken unless the time axis is interpolated
        
        '''
        surf_bin = np.nanargmax(self.data_radar, axis=0)
        surf_time = np.interp(surf_bin,np.arange(0,len(self.time_fast)),self.time_fast)
        return surf_bin , surf_time

So far my hack has been using the elevation compensation, assigning 0 values for NaNs surf_bin = np.nanargmax(np.nan_to_num(self.data_radar), axis=0), and finding the approximate modal value of the surface bin and using a threshold around it for boundaries

from scipy import stats
surf_bin_mode = stats.mode(np.around(radar_dat.surf_bin,-2)).mode[0]
upper = surf_bin_mode - 1000
lower = surf_bin_mode + 500
radar_subset = radar_dat.data_radar[upper:lower, :]

[kingjml]

Patched in the awi_fix_surface branch.

Smoothing is now available to deal with outliers. It uses a median filter with a prescribed window. Seems to handle the outlier end traces and NaN traces (with additional handling)

pySnowRadar/pySnowRadar/snowradar.py

Lines 191 to 192 in 591df7c

	if smooth:
	surf_bin = signal.medfilt(surf_bin, window_len)

For now we need additional handling where the entire trace is nan. The exception triggers infilling of the radar data with 0 values. This does not modify the snow radar object and must be used in conjunction with the smoothing to create a decent surface. It will displace a warning when this happens.

pySnowRadar/pySnowRadar/snowradar.py

Lines 191 to 192 in 591df7c

	if smooth:
	surf_bin = signal.medfilt(surf_bin, window_len)

This is all a bit too temporary and needs to be addressed in the long term as part of #2

I'd rather not roll these into master until the QA is finalized.

[ajutila]

Great, this indeed partly fixes the special surface issues with AWI data. Unfortunately, there's (at least) a third case, where due to descend/ascend of the plane the surf_bin at some part of the frame reaches fast time bin values where some of the traces are NaNs. I dare you @kingjml to try the frame 20170410_01_002/010, if you still have them. I haven't noticed this earlier, as my hack approach of using elevation compensation prior to choosing the subset with a threshold mostly got it decently right.

[ajutila]

Huh, and I guess smoothing doesn't do really well if the weird surf_bin values are too close to the start/end of the frame? Case in point: Data_img_01_20190410_03_009.mat from the TVC flight, although that frame is otherwise really bad example due to really noisy signal (plane turning around).

[kingjml]

Long term solution for this is to change the way processing.py handles individual radar traces. This will require it to drop np.nan values and evaluate only the valid data. The problem here is that we need a method to track how we are manipulating the arrays so that the picker index returns make sense. We'll tackle this before the alpha is done.

Masking might be the key.

[ajutila]

@kingjml & @m9brady, something that I noticed this afternoon. I guess the null space removal in the function get_bounds triggered by the m_above = None option won't work properly.

pySnowRadar/pySnowRadar/snowradar.py

Lines 187 to 208 in 5a2100d

	def get_bounds(self, m_above=None, m_below=5):
	'''
	Get bin numbers where there is valid data (non-nan)
	A threshold can be supplied
	Arguments:
	m_above: bin padding above the signal (?)
	m_below: bin padding below the signal (?)
	Outputs:
	null_lower: the lower bin number bound for use in data-subsetting
	null_upper: the upper bin number bound for use in data-subsetting
	'''
	if m_above:
	null_lower = self.surf_bin.max() + (m_below / self.dfr).astype(int)
	null_upper = self.surf_bin.min() - (m_above / self.dfr).astype(int)
	else:
	null_space = np.argwhere(np.isnan(self.data_radar))[:,0]
	null_upper = null_space[null_space < self.surf_bin.min()].min()
	null_lower = null_space[null_space > self.surf_bin.max()].max()
	return null_lower, null_upper

That's because of zero-patching all-NaN traces in get_surface introduced in awi_fix_surface branch

pySnowRadar/pySnowRadar/snowradar.py

Lines 176 to 196 in 591df7c

	def get_surface(self, smooth=False, window_len = 5):
	'''
	Simple surface tracker based on maximum
	This should be refined and is largely a place holder
	TODO: surf_time is broken unless the time axis is interpolated
	smooth applies a median filter of window_len length
	'''
	try:
	surf_bin = np.nanargmax(self.data_radar, axis=0)
	except ValueError as e:
	surf_bin = np.nanargmax(np.nan_to_num(self.data_radar), axis=0)
	print('Warning: Radar contains all NaN trace(s), patched with 0s to estimate surface')
	if smooth:
	surf_bin = signal.medfilt(surf_bin, window_len)
	surf_time = np.interp(surf_bin,np.arange(0,len(self.time_fast)),self.time_fast)
	return surf_bin , surf_time

That got me thinking, could we then just ignore the zero surf_bin values when assigning the bounds? Just simply

null_lower = self.surf_bin[self.surf_bin > 0].max() + (m_below / self.dfr).astype(int)
null_upper = self.surf_bin[self.surf_bin > 0].min() - (m_above / self.dfr).astype(int)

This should cover many cases, if not all? Or is this too simple to work? As it seems that backflipping (pun intended) nyquist zone crossings in connection to the EM-Bird ascends/descends in AWI data might cause bigger all-NaN patches, smoothing surf_bin in get_surface won't work anymore either. So going from this (elevation compensated, but not nyquist zone backflipped)

via nyquist zone backflipping (not elevation compensated yet)

to this (now it's both elevation compensated and nyquist zone backflipped)

Of course the picker still doesn't like NaN data...

(Oh and for the time being ignore the radar data quality on the nyquist zone backflipped bit. I guess it should be flipped before coherent noise removal and deconvolution, so yay...)

But in essence, would this make sense? If/when switching to dynamic handling of individual radar traces this becomes obsolete, right?

[kingjml]

This emails an AWI issue with calibration flights. Will move to close as a general issue.