Merge pull request #13 from dionhaefner/dynamic-windows

Add dynamic windows and directionality index

fowd/constants.py

@@ -23,7 +23,19 @@
 )
 
 # sea state aggregation periods (in minutes)
-SEA_STATE_INTERVALS = (10, 30)
+SEA_STATE_INTERVALS = (10, 30, 'dynamic')
+
+# settings for dynamic sea state aggregation period
+# the shortest and longest permitted dynamic window sizes (in minutes)
+DYNAMIC_WINDOW_LENGTH_BOUNDS = (10, 60)
+# time in minutes after which the dynamic window size is re-computed
+DYNAMIC_WINDOW_UPDATE_FREQUENCY = 30
+# dynamic window calculation is based on this many minutes before and after current wave
+DYNAMIC_WINDOW_REFERENCE_PERIOD = (12 * 60, 0)
+# total number of different window lengths to try within length bounds
+NUM_DYNAMIC_WINDOWS = 11
+# number of times each window is re-evaluated for better statistics
+NUM_DYNAMIC_WINDOW_SAMPLES = 10
 
 # gravitational acceleration in m/s^2
 GRAVITY = 9.81
@@ -32,13 +44,20 @@
 DENSITY = 1024
 
 # QC thresholds
-QC_FLAG_A_THRESHOLD = 25  # maximum allowed zero-crossing period (in seconds)
-QC_FLAG_B_THRESHOLD = 2  # maximum allowed multiple of limit rate of change
-QC_FLAG_C_THRESHOLD = 10  # maximum allowed number of consecutive identical values
-QC_FLAG_D_THRESHOLD = 8  # maximum allowed exceedance of surface elevation MADN
-QC_FLAG_E_THRESHOLD = 2  # maximum float precision to which sampling rate has to be uniform
-QC_FLAG_F_THRESHOLD = 0.05  # maximum allowed ratio of invalid data
-QC_FLAG_G_THRESHOLD = 100  # minimum number of zero-crossing periods in wave history
+# maximum allowed zero-crossing period (in seconds)
+QC_FLAG_A_THRESHOLD = 25
+# maximum allowed multiple of limit rate of change
+QC_FLAG_B_THRESHOLD = 2
+# maximum allowed number of consecutive identical values
+QC_FLAG_C_THRESHOLD = 10
+# maximum allowed exceedance of surface elevation MADN
+QC_FLAG_D_THRESHOLD = 8
+# maximum float precision to which sampling rate has to be uniform
+QC_FLAG_E_THRESHOLD = 2
+# maximum allowed ratio of invalid data
+QC_FLAG_F_THRESHOLD = 0.05
+# minimum number of zero-crossing periods in wave history
+QC_FLAG_G_THRESHOLD = 100
 
 # QC logging options
 # wave height (in units of sig. wave height) for which failed QC is logged

fowd/operators.py

@@ -56,6 +56,61 @@ def find_wave_indices(z, start_idx=0):
             active = True
 
 
+def compute_dynamic_window_size(t, elevation, min_length, max_length, num_windows, num_samples):
+    """Compute best window size from data by minimizing fluctuations around mean energy drift.
+
+    Reference:
+
+        Fedele, Francesco, et al. “Large Nearshore Storm Waves off the Irish Coast.”
+        Scientific Reports, vol. 9, no. 1, 1, Nature Publishing Group, Oct. 2019, p. 15406.
+        www.nature.com, doi:10.1038/s41598-019-51706-8.
+
+    """
+    best_window_length = None
+    best_window_stationarity = float('inf')
+
+    if isinstance(min_length, np.timedelta64):
+        min_length = int(min_length / np.timedelta64(1, 's'))
+
+    if isinstance(max_length, np.timedelta64):
+        max_length = int(max_length / np.timedelta64(1, 's'))
+
+    window_lengths = np.linspace(min_length, max_length, num_windows, dtype='int')
+
+    for window_length in window_lengths:
+        window_stationarity = []
+        window_offsets = np.linspace(0, window_length, num_samples, endpoint=False, dtype='int')
+
+        for window_offset in window_offsets:
+            if window_offset >= len(t):
+                continue
+
+            current_time_idx = window_offset
+            window_stds = []
+            while True:
+                next_time = t[current_time_idx] + np.timedelta64(window_length, 's')
+                if next_time > t[-1]:
+                    break
+
+                next_time_idx = get_time_index(next_time, t)
+                window_stds.append(np.nanstd(elevation[current_time_idx:next_time_idx], ddof=1.5))
+                current_time_idx = next_time_idx
+
+            window_stds = np.asarray(window_stds)
+            window_stationarity.append(np.nanstd(window_stds[1:] / window_stds[:-1] - 1, ddof=1))
+
+        mean_window_stationarity = np.nanmean(window_stationarity)
+
+        if mean_window_stationarity < best_window_stationarity:
+            best_window_length = window_length
+            best_window_stationarity = mean_window_stationarity
+
+    if best_window_length is None:
+        return min_length
+
+    return best_window_length
+
+
 def add_prefix(dic, prefix):
     """Adds a prefix to every key in given dictionary."""
     return {f'{prefix}_{key}': value for key, value in dic.items()}
@@ -145,7 +200,8 @@ def compute_ursell_number(wave_height, wavelength, water_depth):
 
 def integrate(y, x):
     """Computes numerical integral ∫ y(x) dx."""
-    return np.trapz(y, x)
+    mask = np.isfinite(y) & np.isfinite(x)
+    return np.trapz(y[mask], x[mask])
 
 
 def compute_elevation(displacement):
@@ -209,9 +265,8 @@ def compute_spectral_density(elevation, sample_dt):
     nperseg = round(SPECTRUM_WINDOW_SIZE / sample_dt)
     nfft = 2 ** (math.ceil(math.log(nperseg, 2)))  # round to next higher power of 2
     return scipy.signal.welch(
-        elevation, 1 / sample_dt,
+        elevation, 1 / sample_dt, window='hann',
         nperseg=nperseg, nfft=nfft, noverlap=nperseg // 2,
-        window='hann', detrend=False  # data is already detrended
     )
 
 
@@ -297,7 +352,7 @@ def compute_nth_moment(frequencies, wave_spectral_density, n):
 
 
 def compute_mean_wave_period_spectral(zeroth_moment, frequencies, wave_spectral_density):
-    """Compute mean wave period from wave spectral density."""
+    """Compute mean zero-crossing period from wave spectral density."""
     second_moment = compute_nth_moment(frequencies, wave_spectral_density, 2)
     return np.sqrt(zeroth_moment / second_moment)
 
@@ -421,6 +476,23 @@ def compute_dominant_direction(frequencies, direction, energy_density):
     return circular_weighted_average(frequencies, direction, weights=energy_density)
 
 
+def compute_directionality_index(frequencies, spread, spectral_bandwidth, energy_density):
+    """Compute directionality index R.
+
+    Reference:
+
+        Fedele, Francesco, et al. “Large Nearshore Storm Waves off the Irish Coast.”
+        Scientific Reports, vol. 9, no. 1, 1, Nature Publishing Group, Oct. 2019, p. 15406.
+        www.nature.com, doi:10.1038/s41598-019-51706-8.
+
+    """
+    total_squared_spread = (
+        integrate(np.radians(spread) ** 2 * energy_density, frequencies)
+        / integrate(energy_density, frequencies)
+    )
+    return total_squared_spread / (2 * spectral_bandwidth ** 2)
+
+
 # QC
 
 def check_flag_a(zero_crossing_periods, threshold=QC_FLAG_A_THRESHOLD):
@@ -686,9 +758,19 @@ def get_directional_parameters(time, frequencies, directional_spread, mean_direc
     assert len(dominant_directional_spreads) == len(FREQUENCY_INTERVALS)
     assert len(dominant_directions) == len(FREQUENCY_INTERVALS)
 
+    zeroth_moment = compute_nth_moment(frequencies, wave_spectral_density, 0)
+    first_moment = compute_nth_moment(frequencies, wave_spectral_density, 1)
+    spectral_bandwidth = compute_bandwidth_narrowness(
+        zeroth_moment, first_moment, frequencies, wave_spectral_density
+    )
+    directionality_index = compute_directionality_index(
+        frequencies, directional_spread, spectral_bandwidth, wave_spectral_density
+    )
+
     return {
         'sampling_time': time,
         'dominant_spread_in_frequency_interval': dominant_directional_spreads,
         'dominant_direction_in_frequency_interval': dominant_directions,
         'peak_wave_direction': peak_wave_direction,
+        'directionality_index': directionality_index,
     }