mangrovedigital/tide-engine-api
0
1"""The product: a wind-corrected tide FORECAST.
2
3Outputs two curves for a location -- the standard astronomical tide chart and the
4wind-corrected curve -- from recent past through +3 days, riding the real wind
5forecast. Anchored to current observed water level so the slow drift is captured.
6
7 corrected(t) = chart(t) + r0 + [wind_off(t) - wind_off(now)]
8 r0 = latest observed water level - chart (the current real offset)
9
10For an ungauged point, the anchor + wind come from the nearest gauge / the point's
11own wind column; the wind effect is uniform across the connected system (validated).
12Returns HOLD-style status if there's no recent observation or no wind to ride.
13"""
14import os
15import json
16from datetime import datetime, timezone, timedelta
17from . import sources, harmonic, align, windcorr, noaa
18
19HOUR = 3600
20
21# Stage 16b: the ONLY real-time OBSERVED water-level sensors in the coverage
22# water (earned via the NOAA station census, out/evidence/stage16/). Everything
23# else is prediction-only (-> hybrid). These anchor the wind correction to a live
24# observed level, exactly like a USGS gauge. NOTE: the census corrected the
25# amendment — 8725110 "Naples (outer coast)" is PREDICTION-only; the real observed
26# Naples sensor is 8725114 "Naples Bay, north end".
27NOAA_OBS_STATIONS = {"8725520", "8725114"}
28NOAA_OBS_NAMES = {"8725520": "Fort Myers (NOAA CO-OPS 8725520)",
29 "8725114": "Naples Bay (NOAA CO-OPS 8725114)"}
30
31
32def _noaa_observed_series(sid, train):
33 """Live NOAA CO-OPS observed water level (NAVD88) over the training window,
34 chunked to CO-OPS's 31-day request limit, returned as a USGS-compatible
35 Series so build() treats it like any observed gauge."""
36 d0 = datetime.fromisoformat(train[0]); d1 = datetime.fromisoformat(train[1])
37 T, V = [], []
38 cur = d0
39 while cur <= d1:
40 nxt = min(d1, cur + timedelta(days=30))
41 s = noaa.observed_water_level(sid, cur.strftime("%Y%m%d"),
42 nxt.strftime("%Y%m%d"), datum="NAVD")
43 T.extend(s.t); V.extend(s.values)
44 cur = nxt + timedelta(days=1)
45 if not T:
46 raise sources.NoDataError("no NOAA observed water level for %s" % sid)
47 prov = {"source": "NOAA CO-OPS observed water level", "site": sid,
48 "site_name": NOAA_OBS_NAMES.get(sid, sid), "param": "water_level",
49 "param_desc": "NOAA CO-OPS observed water level (NAVD88)",
50 "datum": "NAVD88", "orig_units": "ft"}
51 return sources.Series("noaa_obs:%s" % sid, T, V, "m (NAVD88)", prov)
52
53# H3: empirical confidence bands (out/conf_bands.json), measured per-anchor per
54# lead-bin from the deployed backtest (scripts/build_conf_bands.py). Replaces the
55# old invented linear ramp. Loaded once; regional "_default" is the fallback.
56_CONF_BANDS = None
57_BINS = (("1-24", 1, 24), ("25-48", 25, 48), ("49-72", 49, 72))
58
59
60def _conf_bands():
61 global _CONF_BANDS
62 if _CONF_BANDS is None:
63 p = os.path.join(os.path.dirname(__file__), "..", "out", "conf_bands.json")
64 try:
65 with open(p) as f:
66 _CONF_BANDS = json.load(f)
67 except Exception:
68 _CONF_BANDS = {}
69 return _CONF_BANDS
70
71
72def _band_for_lead(band_rec, lead, M):
73 key = _BINS[-1][0]
74 for k, lo, hi in _BINS:
75 if lo <= lead <= hi:
76 key = k
77 break
78 return round(band_rec[key] * M, 2)
79
80
81def _corr(a, b):
82 """Pearson correlation of two None-aware parallel lists (fit-quality skill)."""
83 pairs = [(a[i], b[i]) for i in range(min(len(a), len(b)))
84 if a[i] is not None and b[i] is not None]
85 if len(pairs) < 30:
86 return None
87 n = len(pairs)
88 ma = sum(x for x, _ in pairs) / n
89 mb = sum(y for _, y in pairs) / n
90 cov = sum((x - ma) * (y - mb) for x, y in pairs)
91 va = sum((x - ma) ** 2 for x, _ in pairs)
92 vb = sum((y - mb) ** 2 for _, y in pairs)
93 if va <= 0 or vb <= 0:
94 return None
95 return cov / (va ** 0.5 * vb ** 0.5)
96
97
98def _iso(e):
99 return datetime.fromtimestamp(e, tz=timezone.utc).strftime("%Y-%m-%d %H:%M")
100
101
102class TideForecaster:
103 def __init__(self, sid, lat, lon, harm, model, prov):
104 self.sid = sid; self.lat = lat; self.lon = lon
105 self.harm = harm; self.model = model; self.prov = prov
106
107 @staticmethod
108 def _fetch_wl(sid, train):
109 """Return the first available water-level param for this gauge.
110 A genuine no-data on a param -> try the next. A network/transport error
111 (UpstreamError) PROPAGATES -> the caller HOLDs; we never silently fall
112 through to a different param (and thus a different datum) on a blip."""
113 for param in ("63160", "72279", "00065", "62620"):
114 try:
115 return sources.usgs_water_level(sid, *train, param=param)
116 except sources.NoDataError:
117 continue
118 # sources.UpstreamError (and any unexpected error) propagates.
119 raise sources.NoDataError("no water-level param for %s" % sid)
120
121 @staticmethod
122 def _observed(sid, window):
123 """Observed water level for an anchor: NOAA CO-OPS for the two census
124 observed sensors, USGS for everything else. Same Series contract."""
125 if sid in NOAA_OBS_STATIONS:
126 return _noaa_observed_series(sid, window)
127 return TideForecaster._fetch_wl(sid, window)
128
129 # H5: rolling training window. 2025-09-01..2026-06-15 is exactly 287 days;
130 # keep that span but roll the END forward to build-time minus the reanalysis
131 # lag, so a spot sold today doesn't silently age. TIDES_TRAIN_END pins the end
132 # (the golden replay sets 2026-06-15 to reproduce the frozen window exactly).
133 TRAIN_SPAN_DAYS = 287
134 TRAIN_LAG_DAYS = int(os.environ.get("TIDES_TRAIN_LAG_DAYS", "10"))
135
136 @classmethod
137 def _default_train(cls):
138 end_env = os.environ.get("TIDES_TRAIN_END")
139 if end_env:
140 end_d = datetime.fromisoformat(end_env).date()
141 else:
142 end_d = (datetime.now(timezone.utc)
143 - timedelta(days=cls.TRAIN_LAG_DAYS)).date()
144 start_d = end_d - timedelta(days=cls.TRAIN_SPAN_DAYS)
145 return (start_d.isoformat(), end_d.isoformat())
146
147 @classmethod
148 def build(cls, sid, lat, lon, train=None):
149 if train is None:
150 train = cls._default_train()
151 wl = cls._observed(sid, train)
152 harm = harmonic.fit(wl.t, wl.values, include_trend=False)
153 met = sources.met_forcing(lat, lon, *train)
154 obs = align.interp(wl.t, wl.values, met.t, 2 * HOUR)
155 astro = harm.predict(met.t)
156 resid = [None if obs[i] is None else obs[i] - astro[i] for i in range(len(met.t))]
157 model, n = windcorr.fit(resid, met.extra["wind_speed"],
158 met.extra["wind_dir_from"], met.values)
159 # fit-quality gate: enough samples AND real training wind-band skill
160 # (corr of high-passed residual vs high-passed model). Below threshold ->
161 # serve the CHART only, wind correction HOLD; never a garbage correction.
162 woff = model.predict(met.extra["wind_speed"], met.extra["wind_dir_from"], met.values)
163 skill = _corr(windcorr.highpass(resid, windcorr.DEFAULT_HP),
164 windcorr.highpass(woff, windcorr.DEFAULT_HP))
165 min_n = int(os.environ.get("TIDES_MIN_NFIT", "500"))
166 min_skill = float(os.environ.get("TIDES_MIN_SKILL", "0.2"))
167 wind_ok = (n >= min_n) and (skill is not None) and (skill >= min_skill)
168 prov = {"gauge": wl.provenance["site_name"],
169 "datum": wl.provenance["datum"], # truthful — from the fetched param
170 "param": wl.provenance["param"],
171 "param_desc": wl.provenance.get("param_desc"),
172 "train": list(train), "n_fit": n,
173 "wind_band_skill": None if skill is None else round(skill, 3),
174 "wind_ok": wind_ok,
175 "chart": "harmonic analysis of the gauge record (NOAA-equivalent)",
176 "wind": met.provenance.get("source", "wind forcing")}
177 # a failed fit-gate keeps the chart but drops the (untrustworthy) model
178 return cls(sid, lat, lon, harm, model if wind_ok else None, prov)
179
180 def forecast(self, now_epoch, history_h=24, horizon_h=72):
181 """Return {status, now, series:[{t, lead_h, tide_ft, corrected_ft, observed_ft,
182 confidence_ft}], provenance}. Times in UTC epoch; levels in feet (anglers)."""
183 now_dt = datetime.fromtimestamp(now_epoch, tz=timezone.utc)
184 # fetch a wide met window so the 14-day high-pass has support
185 s_date = (now_dt - timedelta(days=20)).strftime("%Y-%m-%d")
186 e_date = (now_dt + timedelta(hours=horizon_h) + timedelta(days=1)).strftime("%Y-%m-%d")
187 rec = {"now": now_dt.isoformat(), "provenance": dict(self.prov), "series": []}
188 try:
189 met = sources.met_forcing(self.lat, self.lon, s_date, e_date)
190 except Exception as exc:
191 rec["status"] = f"HOLD: no wind forecast ({exc})"; return rec
192 # provenance.wind reflects the wind ACTUALLY used this call (ERA5 for
193 # historical/replay windows, NWS for live) — not just the training source.
194 rec["provenance"]["wind"] = met.provenance.get("source", rec["provenance"].get("wind"))
195 grid = met.t
196 chart = self.harm.predict(grid)
197 model_ok = self.model is not None # H5: fit-gate may drop the model
198 if model_ok:
199 woff = self.model.predict(met.extra["wind_speed"],
200 met.extra["wind_dir_from"], met.values)
201 else:
202 woff = [None] * len(grid)
203 # recent observed -> anchor r0 at the latest obs at/just before now
204 try:
205 wl = self._observed(self.sid, (s_date, (now_dt + timedelta(days=1)).strftime("%Y-%m-%d")))
206 obs = align.interp(wl.t, wl.values, grid, 2 * HOUR)
207 except Exception:
208 obs = [None] * len(grid)
209 inow = None
210 for i, t in enumerate(grid):
211 if t <= now_epoch + 1800 and obs[i] is not None and (not model_ok or woff[i] is not None):
212 inow = i
213 if inow is None:
214 rec["status"] = "HOLD: no recent observation to anchor to"; return rec
215 # H2: staleness gate. anchor_age_h is ALWAYS reported. A gauge that went
216 # quiet must HOLD — never anchor to days-old water and call it OK.
217 anchor_age_h = (now_epoch - grid[inow]) / HOUR
218 rec["anchor_age_h"] = round(anchor_age_h, 1)
219 stale_h = float(os.environ.get("TIDES_ANCHOR_STALE_H", "3"))
220 if anchor_age_h > stale_h:
221 rec["status"] = "HOLD: anchor observation stale (%.1f h old)" % anchor_age_h
222 return rec
223 r0 = obs[inow] - chart[inow]
224 M = 3.28084 # m -> ft
225 # H3: empirical confidence band, per-anchor if measured else regional default.
226 cb = _conf_bands()
227 band_rec = cb.get(self.sid) or cb.get("_default")
228 default_band = bool(cb) and self.sid not in cb
229 for i, t in enumerate(grid):
230 if t < now_epoch - history_h * HOUR or t > now_epoch + horizon_h * HOUR:
231 continue
232 tide = chart[i]
233 corr = None
234 if woff[i] is not None:
235 corr = chart[i] + r0 + (woff[i] - woff[inow])
236 lead = (t - now_epoch) / HOUR
237 # measured 1-sigma error for this lead bin (not an invented ramp)
238 if lead <= 0:
239 conf = 0.0
240 elif band_rec:
241 conf = _band_for_lead(band_rec, lead, M)
242 else:
243 conf = round((0.03 + 0.012 * (lead / 24.0)) * M, 2) # legacy last resort
244 rec["series"].append({
245 "epoch": int(t), "t": _iso(t), "lead_h": round(lead, 1),
246 "is_now": abs(lead) < 0.5,
247 "tide_ft": round(tide * M, 2),
248 "corrected_ft": None if corr is None else round(corr * M, 2),
249 "observed_ft": round(obs[i] * M, 2) if (obs[i] is not None and lead <= 0) else None,
250 "confidence_ft": conf,
251 "wind_kt": round(met.extra["wind_speed"][i] * 1.94384, 1),
252 "wind_dir_from": round(met.extra["wind_dir_from"][i]),
253 })
254 rec["status"] = "OK"
255 rec["anchor_offset_ft"] = round(r0 * M, 2)
256 if not model_ok:
257 rec["wind_correction"] = (
258 "HOLD: wind-correction fit below quality threshold "
259 "(n_fit=%s, skill=%s) — showing astronomical chart only"
260 % (self.prov.get("n_fit"), self.prov.get("wind_band_skill")))
261 if default_band:
262 # no per-anchor empirical band for this gauge -> using the regional
263 # default; say so honestly (never present an unearned tight band).
264 rec["confidence_default_band"] = True
265 return rec
266