import numpy as np
import pandas as pd
import numba
from scipy.signal import find_peaks
from vtools.functions.filter import cosine_lanczos
__all__ = [
"get_tidal_hl",
"get_tidal_hh_lh",
"get_tidal_ll_hl",
"get_tidal_amplitude",
"get_tidal_hl_zerocrossing",
"get_tidal_phase_diff",
]
[docs]
def get_smoothed_resampled(
df, cutoff_period="2h", resample_period="1min", interpolate_method="pchip"
):
"""Filter then optionally resample a datetime-indexed DataFrame for tidal peak detection.
The cosine_lanczos filter is applied first at the native resolution to remove
high-frequency noise. The series is then resampled to resample_period if it
differs from the input frequency. Filtering before resampling avoids smoothing
interpolated noise.
Args:
df (DataFrame): A single column dataframe indexed by datetime
cutoff_period (str, optional): cutoff period for cosine lanczos filter. Defaults to '2h'.
resample_period (str, optional): Resample to regular period after filtering.
Pass None or the native frequency to skip resampling. Defaults to '1min'.
interpolate_method (str, optional): interpolation method for resampling. Defaults to 'pchip'.
Returns:
DataFrame: filtered (and optionally resampled) dataframe indexed by datetime
"""
# Filter at native resolution first; resampling before filtering would smooth
# interpolated noise rather than the original signal.
df_filtered = cosine_lanczos(df, cutoff_period)
if resample_period is None:
return df_filtered
native_freq = pd.tseries.frequencies.to_offset(df.index.inferred_freq or df.index.freq)
target_freq = pd.tseries.frequencies.to_offset(resample_period)
if native_freq == target_freq:
return df_filtered
dfb = df_filtered.resample(resample_period).bfill()
df_resampled = df_filtered.resample(resample_period).interpolate(method=interpolate_method)
df_resampled[dfb.iloc[:, 0].isna()] = np.nan
return df_resampled
[docs]
@numba.jit(nopython=True)
def lmax(arr):
"""Local maximum: Returns value only when centered on maximum"""
idx = np.argmax(arr)
if idx == len(arr) / 2:
return arr[idx]
else:
return np.nan
[docs]
@numba.jit(nopython=True)
def lmin(arr):
"""Local minimum: Returns value only when centered on minimum"""
idx = np.argmin(arr)
if idx == len(arr) / 2:
return arr[idx]
else:
return np.nan
[docs]
@numba.jit(nopython=True)
def tlmax(arr):
"""return HH(1) or LH (0)"""
idx = np.argmax(arr) # only the first occurence of the maxima is return
# print(arr,idx)
return idx
[docs]
@numba.jit(nopython=True)
def tlmin(arr):
"""return LL(1) or HL(0)"""
idx = np.argmin(arr)
return idx
[docs]
def periods_per_window(moving_window_size: str, period_str: str) -> int:
"""Number of period size in moving window
Args:
moving_window_size (str): moving window size as a string e.g 7H for 7 hour
period_str (str): period as str e.g. 1T for 1 min
Returns:
int: number of periods in the moving window rounded to an integer
"""
return int(
pd.Timedelta(moving_window_size)
/ pd.to_timedelta(pd.tseries.frequencies.to_offset(period_str))
)
[docs]
def tidal_highs_peaks(df, window="25h", prominence=None):
"""Tidal highs using sequential non-overlapping lunar-day windows and scipy find_peaks.
Partitions the series into sequential windows of length *window* (~25h covers
exactly one lunar day, guaranteeing at most two highs per window for a mixed
semidiurnal tide). Within each window ``scipy.signal.find_peaks`` is used to
locate local maxima; these are ranked so the larger is labelled as the higher-high
and the smaller as the lower-high. Returns the same sparse DataFrame format
as :func:`tidal_highs` (column ``"max"``).
This approach is ~400x faster than the rolling-window method for long series
because it visits each sample once rather than W times (W = window width in samples).
Use it after cosine_lanczos filtering when the signal is known to be smooth.
Args:
df (DataFrame): Single-column DataFrame with a regular DatetimeIndex.
window (str, optional): Length of each sequential window. Defaults to ``'25h'``.
prominence (float or None, optional): Minimum peak prominence passed to
``scipy.signal.find_peaks``. ``None`` (default) applies no prominence
filter; the upstream cosine_lanczos filter makes this safe for smooth inputs.
Returns:
DataFrame: Sparse time series of tidal highs with column ``"max"``, indexed
at the times of detected peaks.
"""
col = df.columns[0]
y = df[col].values
window_td = pd.Timedelta(window)
dt = df.index[1] - df.index[0]
win_samples = int(window_td / dt)
# Minimum inter-peak distance within a window: half a semidiurnal period (~6h)
# prevents picking two peaks within the same tidal cycle half.
min_dist = max(1, int(pd.Timedelta("6h") / dt))
kwargs = {"distance": min_dist}
if prominence is not None:
kwargs["prominence"] = prominence
times = []
values = []
n = len(y)
for i0 in range(0, n, win_samples):
i1 = min(i0 + win_samples, n)
seg = y[i0:i1]
# Skip windows that are entirely NaN
if np.all(np.isnan(seg)):
continue
seg_clean = np.where(np.isnan(seg), -np.inf, seg)
peaks, _ = find_peaks(seg_clean, **kwargs)
for p in peaks:
times.append(df.index[i0 + p])
values.append(y[i0 + p])
result = pd.DataFrame({"max": values}, index=pd.DatetimeIndex(times))
result.index.name = df.index.name
return result
[docs]
def tidal_lows_peaks(df, window="25h", prominence=None):
"""Tidal lows using sequential non-overlapping lunar-day windows and scipy find_peaks.
Mirrors :func:`tidal_highs_peaks` for troughs. See that function's docstring for
the rationale. Returns a sparse DataFrame with column ``"min"``.
Args:
df (DataFrame): Single-column DataFrame with a regular DatetimeIndex.
window (str, optional): Length of each sequential window. Defaults to ``'25h'``.
prominence (float or None, optional): Minimum trough prominence. Defaults to ``None``.
Returns:
DataFrame: Sparse time series of tidal lows with column ``"min"``.
"""
col = df.columns[0]
y = df[col].values
window_td = pd.Timedelta(window)
dt = df.index[1] - df.index[0]
win_samples = int(window_td / dt)
min_dist = max(1, int(pd.Timedelta("6h") / dt))
kwargs = {"distance": min_dist}
if prominence is not None:
kwargs["prominence"] = prominence
times = []
values = []
n = len(y)
for i0 in range(0, n, win_samples):
i1 = min(i0 + win_samples, n)
seg = y[i0:i1]
if np.all(np.isnan(seg)):
continue
seg_clean = np.where(np.isnan(seg), np.inf, seg)
peaks, _ = find_peaks(-seg_clean, **kwargs)
for p in peaks:
times.append(df.index[i0 + p])
values.append(y[i0 + p])
result = pd.DataFrame({"min": values}, index=pd.DatetimeIndex(times))
result.index.name = df.index.name
return result
[docs]
def get_tidal_hh_lh(sh):
"""
return HH(1) or LH (0) based from input tide highs (sh) using rolling window of 2 and tlmax function
"""
sth = sh.rolling(2).apply(tlmax, raw=True)
sth.iloc[0] = (
0 if sth.iloc[1, 0] > 0 else 1
) # fill in the first value based on next value
return sth.iloc[:, 0].map({np.nan: "", 0: "LH", 1: "HH"}).astype(str)
[docs]
def get_tidal_ll_hl(sl):
"""return LL(1) or HL (0) based from input tide lows (sl) using rolling window of 2 and tlmin function
"""
stl = sl.rolling(2).apply(tlmin, raw=True)
stl.iloc[0] = (
0 if stl.iloc[1, 0] > 0 else 1
) # fill in the first value based on next value
return stl.iloc[:, 0].map({np.nan: "", 0: "HL", 1: "LL"}).astype(str)
[docs]
def tidal_highs(df, moving_window_size="7h", method="rolling"):
"""Tidal highs (could be upto two highs in a 25 hr period)
Args:
df (DataFrame): a time series with a regular frequency
moving_window_size (str, optional): moving window size to look for highs within. Defaults to '7h'.
method (str, optional): Peak-detection method. ``'rolling'`` (default) uses a
sliding window with numba-accelerated argmax — robust for noisy signals.
``'find_peaks'`` uses sequential 25h windows with scipy find_peaks — ~400x
faster, suitable for pre-filtered smooth signals.
Returns:
DataFrame: an irregular time series with highs at resolution of df.index
"""
if method == "find_peaks":
return tidal_highs_peaks(df)
period_str = df.index.freqstr
periods = periods_per_window(moving_window_size, period_str)
dfmax = df.rolling(moving_window_size, min_periods=periods).apply(lmax, raw=True)
dfmax = dfmax.shift(periods=-(periods // 2 - 1))
dfmax = dfmax.dropna()
dfmax.columns = ["max"]
return dfmax
[docs]
def tidal_lows(df, moving_window_size="7h", method="rolling"):
"""Tidal lows (could be upto two lows in a 25 hr period)
Args:
df (DataFrame): a time series with a regular frequency
moving_window_size (str, optional): moving window size to look for lows within. Defaults to '7h'.
method (str, optional): Peak-detection method. ``'rolling'`` (default) uses a
sliding window with numba-accelerated argmin — robust for noisy signals.
``'find_peaks'`` uses sequential 25h windows with scipy find_peaks — ~400x
faster, suitable for pre-filtered smooth signals.
Returns:
DataFrame: an irregular time series with lows at resolution of df.index
"""
if method == "find_peaks":
return tidal_lows_peaks(df)
period_str = df.index.freqstr
periods = periods_per_window(moving_window_size, period_str)
dfmin = df.rolling(moving_window_size, min_periods=periods).apply(lmin, raw=True)
dfmin = dfmin.shift(periods=-(periods // 2 - 1))
dfmin = dfmin.dropna()
dfmin.columns = ["min"]
return dfmin
[docs]
def get_tidal_hl(
df,
cutoff_period="2h",
resample_period="1min",
interpolate_method="pchip",
moving_window_size="7h",
method="rolling",
):
"""Get Tidal highs and lows
Args:
df (DataFrame): A single column dataframe indexed by datetime
cutoff_period (str, optional): cutoff period for cosine lanczos filter. Defaults to '2h'.
resample_period (str, optional): Resample to regular period after filtering.
Pass None to skip resampling and operate at the native input frequency.
Ignored when method='find_peaks'. Defaults to '1min'.
interpolate_method (str, optional): interpolation for resampling. Defaults to 'pchip'.
moving_window_size (str, optional): moving window size for rolling method. Defaults to '7h'.
method (str, optional): Peak-detection method. ``'rolling'`` (default) uses the
numba-accelerated sliding-window approach \u2014 robust for noisy signals.
``'find_peaks'`` filters at native resolution then uses sequential 25h
windows with scipy find_peaks \u2014 ~400x faster for pre-filtered smooth inputs
such as SF tidal records.
Returns:
tuple of DataFrame: Tidal high and tidal low time series
"""
if method == "find_peaks":
# Filter at native resolution; no resampling needed before find_peaks.
dfs = cosine_lanczos(df, cutoff_period)
return tidal_highs(dfs, method="find_peaks"), tidal_lows(dfs, method="find_peaks")
dfs = get_smoothed_resampled(df, cutoff_period, resample_period, interpolate_method)
return tidal_highs(dfs, method="rolling"), tidal_lows(dfs, method="rolling")
get_tidal_hl_rolling = get_tidal_hl # for older refs. #FIXME
[docs]
def get_tidal_amplitude(dfh, dfl):
"""Tidal amplitude given tidal highs and lows
Args:
dfh (DataFrame): Tidal highs time series
dfl (DataFrame): Tidal lows time series
Returns:
DataFrame: Amplitude timeseries, at the times of the low following the high being used for amplitude calculation
"""
dfamp = pd.concat([dfh, dfl], axis=1)
dfamp = dfamp[["min"]].dropna().join(dfamp[["max"]].ffill())
return pd.DataFrame(dfamp["max"] - dfamp["min"], columns=["amplitude"])
[docs]
def get_tidal_amplitude_diff(dfamp1, dfamp2, percent_diff=False, tolerance="4h"):
"""Get the difference of values within +/- 4H of values in the two amplitude arrays
Args:
dfamp1 (DataFrame): Amplitude time series
dfamp2 (DataFrame): Amplitude time series
percent_diff (bool, optional): If true do percent diff. Defaults to False.
Returns:
DataFrame: Difference dfamp1-dfamp2 or % Difference (dfamp1-dfamp2)/dfamp2*100 for values within +/- 4H of each other
"""
dfamp = pd.merge_asof(
dfamp1,
dfamp2,
left_index=True,
right_index=True,
direction="nearest",
tolerance=pd.Timedelta(tolerance),
)
if percent_diff:
dfdiff = 100.0 * (dfamp.iloc[:, 0] - dfamp.iloc[:, 1]) / dfamp.iloc[:, 1]
else:
dfdiff = dfamp.iloc[:, 0] - dfamp.iloc[:, 1]
return pd.DataFrame(dfdiff, columns=["amplitude_diff"])
[docs]
def get_phase_diff(df1, df2, tolerance="4h"):
df1["time"] = df1.index
df2["time"] = df2.index
df21 = pd.merge_asof(
df2,
df1,
left_index=True,
right_index=True,
direction="nearest",
tolerance=pd.Timedelta(tolerance),
)
return (df21["time_x"] - df21["time_y"]).apply(lambda x: x.total_seconds() / 60)
[docs]
def get_tidal_phase_diff(dfh2, dfl2, dfh1, dfl1, tolerance="4h"):
"""Calculates the phase difference between df2 and df1 tidal highs and lows
Scans +/- 4 hours in df1 to get the highs and lows in that windows for df2 to
get the tidal highs and lows at the times of df1
Args:
dfh2 (DataFrame): Timeseries of tidal highs
dfl2 (DataFrame): Timeseries of tidal lows
dfh1 (DataFrame): Timeseries of tidal highs
dfl1 (DataFRame): Timeseries of tidal lows
Returns:
DataFrame: Phase difference (dfh2-dfh1) and (dfl2-dfl1) in minutes
"""
high_phase_diff = get_phase_diff(dfh2, dfh1, tolerance)
low_phase_diff = get_phase_diff(dfl2, dfl1, tolerance)
merged_diff = pd.merge(
pd.DataFrame(high_phase_diff, index=dfh1.index),
pd.DataFrame(low_phase_diff, index=dfl1.index),
how="outer",
left_index=True,
right_index=True,
)
return merged_diff.iloc[:, 0].fillna(merged_diff.iloc[:, 1])
[docs]
def get_tidal_hl_zerocrossing(df, round_to="1min"):
"""
Finds the tidal high and low times using zero crossings of the first derivative.
This works for all situations but is not robust in the face of noise and perturbations in the signal
"""
zc, zi = zerocross(df)
if round_to:
zc = pd.to_datetime(zc).round(round_to)
return zc
[docs]
def zerocross(df):
"""
Calculates the gradient of the time series and identifies locations where gradient changes sign
Returns the time rounded to nearest minute where the zero crossing happens (based on linear derivative assumption)
"""
diffdfv = pd.Series(np.gradient(df[df.columns[0]].values), index=df.index)
indi = np.where((np.diff(np.sign(diffdfv))) & (diffdfv[1:] != 0))[0]
# Find the zero crossing by linear interpolation
zdb = diffdfv[indi].index
zda = diffdfv[indi + 1].index
x = diffdfv.index
y = diffdfv.values
dx = x[indi + 1] - x[indi]
dy = y[indi + 1] - y[indi]
zc = -y[indi] * (dx / dy) + x[indi]
return zc, indi
##---- FUNCTIONS CACHED BELOW THIS LINE PERHAPS TO USE LATER? ---#
[docs]
def where_changed(df):
""" """
diff = np.diff(df[df.columns[0]].values)
wdiff = np.where(diff != 0)[0]
wdiff = np.insert(wdiff, 0, 0) # insert the first value i.e. zero index
return df.iloc[wdiff + 1, :]
[docs]
def where_same(dfg, df):
"""
return dfg only where its value is the same as df for the same time stamps
i.e. the interesection locations with df
"""
dfall = pd.concat([dfg, df], axis=1)
return dfall[dfall.iloc[:, 0] == dfall.iloc[:, 1]].iloc[:, 0]
[docs]
def limit_to_indices(df, si, ei):
return df[(df.index > si) & (df.index < ei)]
[docs]
def filter_where_na(df, dfb):
"""
remove values in df where dfb has na values
"""
dfx = dfb.loc[df.index]
return df.loc[dfx.dropna().index, :]
