Legacy plus depth class (#9)

* Added in some functionality for legacy data w/o acc. Added in a class for managing depth timeseries.

* Added legacy data file

* Added in more functionality around which depth profile to choose. Working on adding more sophisticated depth class for lyte profile to use

* Added in ZPFO function. Added in error detection in the profile report card. Added some constraints on the depth data when seeking a surface.

* Added mild variation of the fusing of depths

* Transitioning to depth clases

* Added in depth timeseries classes with a simple test

* Working fusion with depth classes

* Found minor bug in blending technique

* updated tests to account for latest changes in depth computation

* Fine tuning somethings. Narrowed the NIR stop detection

* Added pressure as an option

* Added in margin of error related calculations

* Added more distance metrics such as average distance travelled and distance during motion

* Improving on surface detection

* Added in various test data capturing certain scenarions

* Added in some legacy bits, added logging, added some legacy test. Other tests are still failing though

* Finalized test_detect.py

* Loosened requirements for profile tests to be less sensitive to detection modifications.

* Added test datasets

* Updated docs to actually make an api. Cleaned up some old tests that were commented out.

* Updated the pypi manifest to ignore test data

* Added test for pressure profile property

* Removed plotting of the data

.gitignore

@@ -106,3 +106,5 @@ ENV/
 # IDE settings
 .vscode/
 .idea/
+
+docs/api/*

MANIFEST.in

@@ -4,9 +4,9 @@ include HISTORY.rst
 include LICENSE
 include README.rst
 
-recursive-include tests *
+recursive-exclude tests *
 recursive-exclude * __pycache__
 recursive-exclude * *.py[co]
 recursive-exclude * *.csv
 
-recursive-include docs *.rst conf.py Makefile make.bat *.jpg *.png *.gif
+recursive-exclude docs *.rst conf.py Makefile make.bat *.jpg *.png *.gif

docs/conf.py

@@ -33,6 +33,7 @@
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
 extensions = ['sphinx.ext.autodoc','sphinxcontrib.apidoc', 'sphinx.ext.viewcode', 'nbsphinx', 'sphinx_gallery.load_style','sphinx.ext.autosectionlabel']
 
+apidoc_module_dir = '../study_lyte'
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']
 

study_lyte/adjustments.py

@@ -1,6 +1,6 @@
 import numpy as np
 import pandas as pd
-
+from scipy.signal import lfilter
 
 def get_points_from_fraction(n_samples, fraction):
     """
@@ -27,7 +27,7 @@ def get_directional_mean(arr: np.array, fractional_basis: float = 0.01, directio
     return avg
 
 
-def get_neutral_bias_at_border(series: pd.Series, fractional_basis: float = 0.01, direction='forward'):
+def get_neutral_bias_at_border(series: pd.Series, fractional_basis: float = 0.005, direction='forward'):
     """
     Bias adjust the series data by using the XX % of the data either at the front of the data
     or the end of the .
@@ -65,6 +65,31 @@ def get_normalized_at_border(series: pd.Series, fractional_basis: float = 0.01,
         border_norm = series
     return border_norm
 
+def merge_on_to_time(df_list, final_time):
+    """
+    Reindex the df fram the list onto a final time stamp
+    """
+    # Build dummy result in case no data is passed
+    result = pd.DataFrame()
+
+    # Merge everything else to it
+    for i, df in enumerate(df_list):
+        time_df = df.copy()
+        if df.index.name != 'time':
+            time_df = time_df.set_index('time')
+        else:
+            time_df = df.copy()
+
+        data = time_df.reindex(time_df.index.union(final_time)).interpolate(method='cubic').reindex(final_time)
+        if i == 0:
+            result = data
+        else:
+            result = pd.merge_ordered(result, data, on='time', fill_method='cubic')
+
+    # interpolate the nan's
+    result = result.interpolate(method='nearest', limit_direction='both')
+    return result
+
 
 def merge_time_series(df_list):
     """
@@ -87,7 +112,7 @@ def merge_time_series(df_list):
         if i == 0:
             result = df.copy()
         else:
-            result = pd.merge_ordered(result, df, on='time')
+            result = pd.merge_ordered(result, df, on='time', fill_method='cubic')
 
     # interpolate the nan's
     result = result.interpolate(method='index')
@@ -247,4 +272,19 @@ def convert_force_to_pressure(force, tip_diameter_m, geom_adj=1):
     # convert to pressure in Pascals
     pressure = force.div(area)
     # Adjust for shape and convert to kPa
-    return pressure * geom_adj / 1000
+    return pressure * geom_adj / 1000
+
+def zfilter(series, fraction):
+    """
+    Zero phase filter
+    """
+    window = get_points_from_fraction(len(series), fraction)
+    filter_coefficients = np.ones(window) / window
+
+    # Apply the filter forward
+    filtered_signal = lfilter(filter_coefficients, 1, series)
+
+    # Apply the filter backward
+    filtered_signal = lfilter(filter_coefficients, 1, filtered_signal[::-1])[::-1]
+
+    return filtered_signal

study_lyte/cropping.py

@@ -46,7 +46,7 @@ def crop_to_snow(df: pd.DataFrame, active_col='Sensor3', ambient_col='Sensor2',
         cropped: pd.Dataframe cropped to the time period where a surface was detected to the end
     """
     df['nir'] = remove_ambient(df[active_col],  df[ambient_col])
-    surface = get_nir_surface(df['nir'].values, **kwargs)
+    surface = get_nir_surface(df['nir'], **kwargs)
     cropped = df.iloc[surface:]
     return cropped
 

study_lyte/depth.py

@@ -1,11 +1,11 @@
 import pandas as pd
 from scipy.integrate import cumtrapz
 import numpy as np
+from types import SimpleNamespace
 
 from .decorators import time_series
-from .adjustments import get_neutral_bias_at_border, assume_no_upward_motion
-from .detect import get_acceleration_stop, get_acceleration_start, first_peak, nearest_valley, nearest_peak
-
+from .detect import nearest_peak
+from .adjustments import zfilter
 
 @time_series
 def get_depth_from_acceleration(acceleration_df: pd.DataFrame) -> pd.DataFrame:
@@ -20,7 +20,9 @@ def get_depth_from_acceleration(acceleration_df: pd.DataFrame) -> pd.DataFrame:
     Return:
         position_df: pandas Dataframe containing the same input axes plus magnitude of the result position
     """
-    # Auto gather the x,y,z acceleration columns if they're there.
+    if type(acceleration_df) != pd.DataFrame:
+        acceleration_df = pd.DataFrame(acceleration_df)
+        # Auto gather the x,y,z acceleration columns if they're there.
     acceleration_columns = [c for c in acceleration_df.columns if 'Axis' in c or 'acceleration' == c]
 
     # Convert from g's to m/s2
@@ -87,7 +89,6 @@ def get_constrained_baro_depth(baro_depth, start, stop, method='nanmedian'):
     top_search = baro_depth.iloc[:mid]
     default_top = np.where(top_search == top_search.max())[0][0]
     top = nearest_peak(baro_depth.values, start, default_index=default_top, height=-10, distance=100)
-    # top = nearest_peak(df[baro].values, start, default_index=max_out, height=-0.1, distance=100)
 
     # Find valleys after, select closest to midpoint
     soft_stop = mid + int(0.1 * n_points)
@@ -122,10 +123,144 @@ def get_constrained_baro_depth(baro_depth, start, stop, method='nanmedian'):
     constrained = result.set_index('time')
     #assume_no_upward_motion(result[baro])
     constrained = constrained - constrained.iloc[0]
+
     # from .plotting import plot_constrained_baro
     # pos = get_depth_from_acceleration(df).mul(100)
     # pos = pos.reset_index()
     # plot_constrained_baro(df, result, const, pos, top, bottom, start, stop,
     #                       baro=baro, acc_axis=acc_axis)
 
     return constrained
+
+
+class DepthTimeseries:
+    """
+    Class for managing depth time series data
+    """
+    def __init__(self, series, start_idx=None, stop_idx=None, origin=None):
+        # Hang on to the raw data
+        self.raw = series
+
+        # Keep track of the start stop
+        self.start_idx = start_idx
+        self.stop_idx = stop_idx
+
+        # Establish a zero depth index
+        self.origin = origin
+        if self.origin is None:
+            self.origin = start_idx
+
+
+        # Holder for depth data with at least the origin zeroed out
+        self._depth = None
+
+        # Useful attributes
+        self._avg_velocity = None
+        self._distance_traveled = None
+        self._distance_traveled_during_motion = None
+        self._avg_distance_traveled = None
+
+        self._has_upward_motion = None
+        self._velocity = None
+        self._velocity_range = None
+        self._max_velocity = None
+
+    @property
+    def depth(self):
+        if self._depth is None:
+            self._depth = self.raw - self.raw.iloc[self.origin]
+        return self._depth
+
+    @property
+    def velocity(self):
+        if self._velocity is None:
+            dt = self.depth.index[1] - self.depth.index[0]
+            velocity = np.gradient(self.depth, dt)
+            # Due to rounding issues the index is not evenly space, so filter the velocity
+            velocity = zfilter(velocity, 0.01)
+            self._velocity = pd.Series(velocity, self.depth.index, name='velocity')
+        return self._velocity
+
+    @property
+    def velocity_range(self):
+        """min, max of the absulute probe velocity during motion"""
+        if self._velocity_range is None:
+            minimum = np.min(self.velocity.iloc[self.start_idx:self.stop_idx].abs())
+            self._velocity_range = SimpleNamespace(min=minimum, max=self.max_velocity)
+        return self._velocity_range
+
+    @property
+    def avg_velocity(self):
+        if self._avg_velocity is None:
+            self._avg_velocity = abs(self.velocity.iloc[self.start_idx:self.stop_idx].mean())
+        return self._avg_velocity
+
+    @property
+    def max_velocity(self):
+        """Max velocity between start and stop"""
+        if self._max_velocity is None:
+            self._max_velocity = abs(self.velocity.iloc[self.start_idx:self.stop_idx].min())
+        return self._max_velocity
+
+    @property
+    def distance_traveled(self):
+        """Total distance traveled"""
+        if self._distance_traveled is None:
+            self._distance_traveled = abs(self.depth.max() - self.depth.min())
+        return self._distance_traveled
+
+    @property
+    def avg_distance_traveled(self):
+        """Average distance traveled"""
+        if self._avg_distance_traveled is None:
+            self._avg_distance_traveled = abs(self.depth.iloc[0:self.start_idx].mean() -
+                                              self.depth.iloc[self.stop_idx:].mean())
+        return self._avg_distance_traveled
+
+
+    @property
+    def distance_traveled_during_motion(self):
+        """Total distance traveled between start and stop"""
+        if self._distance_traveled_during_motion is None:
+            self._distance_traveled_during_motion = abs(self.depth.iloc[self.start_idx] - self.depth.iloc[self.stop_idx])
+        return self._distance_traveled_during_motion
+
+    @property
+    def has_upward_motion(self):
+        """Contains upward motion in between the start and stop"""
+        if self._has_upward_motion is None:
+            self._has_upward_motion = False
+            # crop the depth data and downsample for speedy check
+            data = self.raw.iloc[self.start_idx:self.stop_idx]
+            if len(data) > 1000:
+                coarse = data.groupby(data.index // 200).first()
+            else:
+                coarse = data
+            # loop and find any values greater than the current value
+            for i, v in enumerate(coarse):
+                upward = np.any(coarse.iloc[i:] > v + 5)
+                if upward:
+                    self._has_upward_motion = True
+                    break
+
+        return self._has_upward_motion
+
+
+class BarometerDepth(DepthTimeseries):
+    @property
+    def depth(self):
+        if self._depth is None:
+            self._depth = get_constrained_baro_depth(self.raw, self.start_idx, self.stop_idx, method='nanmean')['baro']
+            self._depth = self._depth.reindex(self.raw.index, method='nearest')
+        return self._depth
+
+
+class AccelerometerDepth(DepthTimeseries):
+    @property
+    def depth(self):
+        if self._depth is None:
+            valid = ~np.isnan(self.raw)
+            self._depth = get_depth_from_acceleration(self.raw[valid])[self.raw.name]
+            self._depth.iloc[self.stop_idx:] = self._depth.iloc[self.stop_idx]
+
+        return self._depth

study_lyte/detect.py

@@ -129,7 +129,7 @@ def get_acceleration_start(acceleration, threshold=-0.01, max_threshold=0.02):
     return acceleration_start
 
 
-def get_acceleration_stop(acceleration, threshold=0.1, max_threshold=0.3):
+def get_acceleration_stop(acceleration, threshold=-0.05, max_threshold=0.05):
     """
     Returns the index of the last value that has a relative change greater than the
     threshold of absolute normalized signal
@@ -141,20 +141,20 @@ def get_acceleration_stop(acceleration, threshold=0.1, max_threshold=0.3):
         acceleration_start: Integer of index in array of the first value meeting the criteria
     """
     acceleration = acceleration.values
-    accel_neutral = -1 * acceleration[~np.isnan(acceleration)]
-    n_points = get_points_from_fraction(len(acceleration), 0.005)
-    max_ind = first_peak(accel_neutral[::-1], height=0.3, distance=10)
-    max_ind = len(acceleration) - max_ind - 1
+    n = get_points_from_fraction(len(acceleration), 0.005)
+    if n > 40:
+        n = 20
 
-    acceleration_stop = get_signal_event(accel_neutral[max_ind:], threshold=threshold,
+    max_ind = np.argwhere(acceleration == acceleration.min())[0][0]
+    acceleration_stop = get_signal_event(acceleration[max_ind:], threshold=threshold,
                                          max_threshold=max_threshold,
-                                         n_points=None,
+                                         n_points=n,
                                          search_direction='backward')
     acceleration_stop = acceleration_stop + max_ind
     return acceleration_stop
 
 
-def get_nir_surface(clean_active, threshold=0.05, max_threshold=0.1):
+def get_nir_surface(clean_active, threshold=-1, max_threshold=0.25):
     """
     Using the cleaned active, estimate the index at when the probe was in the snow.
 
@@ -166,20 +166,27 @@ def get_nir_surface(clean_active, threshold=0.05, max_threshold=0.1):
     Return:
         surface: Integer index of the estimated snow surface
     """
-    clean_norm = clean_active / clean_active.max()
-    clean_norm = clean_norm - abs(clean_norm).min()
-    surface = get_signal_event(clean_norm, search_direction='forward', threshold=threshold,
-                               max_threshold=max_threshold)
+    n = get_points_from_fraction(len(clean_active), 0.01)
+    # Normalize by data unaffected by ambient
+    clean_norm = clean_active / clean_active[n:].mean()
+    neutral = get_neutral_bias_at_border(clean_norm)
+
+    # Retrieve a likely candidate under challenging ambient conditions
+
+    max_idx = np.argwhere((neutral == neutral.min()).values)[0][0]
+
+    # Detect likely candidate normal ambient conditions
+    surface = get_signal_event(neutral, search_direction='forward', threshold=threshold,
+                               max_threshold=max_threshold, n_points=n)
     return surface
 
 
-def get_nir_stop(active, fractional_basis=0.05, max_threshold=0.01, threshold=-0.01):
+def get_nir_stop(active, fractional_basis=0.05, max_threshold=0.008, threshold=-0.05):
     """
     Often the NIR signal shows the stopping point of the probe by repeated data.
     This looks at the active signal to estimate the stopping point
     """
-    # Perform a removal of ambient but using the end of the data.
-    n = get_points_from_fraction(len(active), 0.05)
+    n = get_points_from_fraction(len(active), 0.1)
     border_fract = 0.3
     norm_active = get_normalized_at_border(active, fractional_basis=border_fract, direction='backward')
     norm_active = norm_active.rolling(window=n, center=True, closed='both', min_periods=1).mean()
@@ -194,12 +201,5 @@ def get_nir_stop(active, fractional_basis=0.05, max_threshold=0.01, threshold=-0
                             max_threshold=max_threshold, n_points=n_points)
     stop += ind
 
-    # from .plotting import plot_ts, plt
-    # ax = plot_ts(norm_active, events=[('stop', stop), ('avg_tail', len(norm_active) - n)],
-    #              color='red', show=False)
-    # ax.axhline(threshold)
-    # ax.axhline(max_threshold)
-    # plt.show()
-
     return stop