mangrovedigital/tide-engine-api

1"""Wind correction for the corrected-tide product.
2
3Key requirement: a multi-day sustained blow (e.g. N15 for 3 days after a winter
4front) drains the connected system to a new low and HOLDS it. That is a
5cumulative effect, so the features are TRAILING-WINDOW means of wind stress over
6several timescales (hours to days), not just the instantaneous gust.
7
8Decomposition that makes it transferable and bias-free:
9  * high-pass everything at a LONG window (default 14 days) -> removes only the
10    slow seasonal/steric drift, but KEEPS multi-day wind setdowns.
11  * wind enters as eastward/northward stress (~speed^2), trailing-averaged over
12    [3,6,12,24,48,72] h, so both fast response and sustained drain are captured.
13  * pressure as an inverse-barometer term (instant + 24 h mean).
14  * no intercept (an intercept re-introduces period bias).
15
16corrected_level(t) = astronomical_tide(t) + wind_offset(t)
17"""
18import math
19from . import linalg
20
21P_REF = 1013.25
22DEFAULT_WINDOWS = [3, 6, 12, 24, 48, 72]
23DEFAULT_HP = 336  # hours (14 days)
24
25
26def _stress_components(speed, dirf):
27    r = math.radians(dirf)
28    u = -speed * math.sin(r)   # eastward velocity (toward)
29    v = -speed * math.cos(r)   # northward velocity (toward)
30    return speed * u, speed * v  # sx, sy  (~ speed^2, signed by direction)
31
32
33def running_mean(x, win):
34    """Centered running mean (None-aware) — used only for the long high-pass."""
35    n = len(x); half = win // 2; out = [None] * n
36    for i in range(n):
37        lo = max(0, i - half); hi = min(n, i + half + 1)
38        v = [x[j] for j in range(lo, hi) if x[j] is not None]
39        if len(v) >= max(3, win // 3):
40            out[i] = sum(v) / len(v)
41    return out
42
43
44def highpass(x, win):
45    rm = running_mean(x, win)
46    return [None if (x[i] is None or rm[i] is None) else x[i] - rm[i]
47            for i in range(len(x))]
48
49
50def trailing_mean(x, W):
51    """Causal mean over the last W samples (hourly). None if any are missing."""
52    n = len(x); out = [None] * n
53    for i in range(n):
54        if i + 1 < W:
55            continue
56        seg = x[i - W + 1:i + 1]
57        if any(s is None for s in seg):
58            continue
59        out[i] = sum(seg) / W
60    return out
61
62
63class WindCorr:
64    def __init__(self, coeffs, windows, hp_win):
65        self.coeffs = coeffs
66        self.windows = windows
67        self.hp_win = hp_win
68
69    def _feature_arrays(self, speed, dirf, pres):
70        n = len(speed)
71        SX = [0.0] * n; SY = [0.0] * n
72        for i in range(n):
73            SX[i], SY[i] = _stress_components(speed[i], dirf[i])
74        SXh = highpass(SX, self.hp_win); SYh = highpass(SY, self.hp_win)
75        dPh = highpass(pres, self.hp_win)
76        cols = []
77        names = []
78        for W in self.windows:
79            cols.append(trailing_mean(SXh, W)); names.append(f"sx_{W}h")
80            cols.append(trailing_mean(SYh, W)); names.append(f"sy_{W}h")
81        cols.append(dPh); names.append("dP_inst")
82        cols.append(trailing_mean(dPh, 24)); names.append("dP_24h")
83        return cols, names
84
85    def rows(self, speed, dirf, pres):
86        cols, _ = self._feature_arrays(speed, dirf, pres)
87        n = len(speed); out = []
88        for i in range(n):
89            r = [c[i] for c in cols]
90            out.append(None if any(v is None for v in r) else r)
91        return out
92
93    def predict(self, speed, dirf, pres):
94        rows = self.rows(speed, dirf, pres)
95        return [None if r is None else sum(r[k] * self.coeffs[k]
96                                           for k in range(len(self.coeffs)))
97                for r in rows]
98
99    def feature_names(self):
100        _, names = self._feature_arrays([0.0], [0.0], [P_REF])
101        return names
102
103    def sustained_response(self, speed=7.7):
104        """Steady-state offset (m) for a SUSTAINED wind (default ~15 kt) FROM each
105        compass dir — what the bays settle to after a multi-day blow."""
106        out = {}
107        for name, deg in {"N":0,"NE":45,"E":90,"SE":135,"S":180,
108                          "SW":225,"W":270,"NW":315}.items():
109            sx, sy = _stress_components(speed, deg)
110            s = 0.0
111            k = 0
112            for _W in self.windows:
113                s += self.coeffs[k] * sx + self.coeffs[k + 1] * sy
114                k += 2
115            out[name] = s
116        return out
117
118
119def fit(residual, speed, dirf, pres, windows=None, hp_win=DEFAULT_HP, ridge=1.0):
120    # ridge=1.0 (relative to mean diagonal) tuned on held-out winter fronts: tempers
121    # the multi-day gain so the model catches dumps without overshooting (see
122    # scripts/tune_ridge.py -- forecast RMSE 0.140 m, 20% better than persistence).
123    windows = windows or DEFAULT_WINDOWS
124    model = WindCorr(None, windows, hp_win)
125    tgt = highpass(residual, hp_win)
126    rows = model.rows(speed, dirf, pres)
127    X, y = [], []
128    for i in range(len(tgt)):
129        if tgt[i] is None or rows[i] is None:
130            continue
131        X.append(rows[i]); y.append(tgt[i])
132    if len(X) < len(model.feature_names()) * 3:
133        raise RuntimeError("too few samples for windcorr fit (%d)" % len(X))
134    model.coeffs = linalg.solve_lstsq(X, y, ridge=ridge)
135    return model, len(X)
136