mangrovedigital/tide-engine-api
0
1"""
2Data sources for the wind-corrected tide engine.
3
4Single responsibility: fetch REAL data from official / managed APIs, attach
5provenance, cache raw responses to disk, and return clean (utc_epoch, value)
6series. No modelling here.
7
8Sources (all commercially clean — public domain or commercial-OK w/ attribution)
9-------
10- USGS NWIS Instantaneous Values (US gov, public domain) -> ground-truth
11 water-surface elevation (NAVD88, param 63160) at gauges inside the target area.
12- Wind + surface pressure forcing at the target lat/lon, by window age:
13 * historical / training -> ERA5 reanalysis (Copernicus CDS). Free for
14 commercial use WITH attribution. Retrieval is async NetCDF and lives in a
15 BUILD-TIME tool (scripts/era5_fetch.py) that writes a processed hourly
16 point-cache; this runtime module reads that cache with stdlib only.
17 * live / forecast -> NWS api.weather.gov (US gov, public domain):
18 gridpoint wind forecast + station observations. NWS does not forecast
19 surface pressure at these grids -> the pressure term is degraded (held at
20 the last observed value); harmless here because windcorr high-passes at
21 14 days and the product is a woff(t)-woff(now) delta (constant cancels).
22
23Open-Meteo (previously used for both roles) is REMOVED: its free tier is
24CC-BY-NC (non-commercial), which is not sale-clean.
25
26Every fetch returns a Series with .provenance describing exactly where the
27numbers came from and when they were pulled.
28"""
29
30import json
31import os
32import time
33import hashlib
34import math
35import random
36import tempfile
37import threading
38import urllib.parse
39from datetime import datetime, timezone, timedelta
40
41import requests
42
43CACHE_DIR = os.path.join(os.path.dirname(__file__), "..", "data")
44USGS_IV = "https://waterservices.usgs.gov/nwis/iv/"
45NWS_API = "https://api.weather.gov"
46NWS_UA = "tides-engine (contact: forprophetoc@gmail.com)"
47ERA5_ATTRIBUTION = ("Generated using Copernicus Climate Change Service "
48 "information 2026 (ERA5 hourly data on single levels).")
49# window older than this many days -> reanalysis (ERA5); else live (NWS).
50LIVE_AGE_DAYS = 2
51# FIX 2 (never-blind 2026-07-18): near-real-time observed-water bodies go stale in
52# 15 minutes. Pinned historical/training windows keep ttl=None (immutable data).
53NRT_TTL_SECONDS = int(os.environ.get("TIDES_NRT_TTL_SECONDS", "900"))
54
55
56def _is_nrt_end(end_date):
57 """True when a fetch window reaches near-now — its body must expire so a
58 recovered gauge is noticed without a restart (Oscar's ruling 2026-07-18)."""
59 try:
60 end_d = datetime.fromisoformat(end_date).date()
61 except (ValueError, TypeError):
62 return False
63 return (_today_utc_date() - end_d).days <= LIVE_AGE_DAYS
64
65
66class UpstreamError(RuntimeError):
67 """Transport/network/5xx failure talking to an upstream — the data MIGHT
68 exist; we just couldn't reach it. Callers must HOLD, never silently fall
69 through to a different param/datum (that would swap the datum on a blip)."""
70
71
72class NoDataError(RuntimeError):
73 """Upstream reached fine but has genuinely no data for this site/param —
74 safe to fall through to the next candidate param."""
75
76
77def _datum_from_desc(desc):
78 """Derive the vertical datum from USGS's own variableDescription text —
79 the truthful source (verified against real fetches, Stage 3):
80 '... above NAVD 1988 ...' -> NAVD88 ; '... NGVD 1929 ...' -> NGVD29 ;
81 'Gage height, feet' -> gage height (local datum, NOT NAVD88)."""
82 d = (desc or "").lower()
83 if "navd" in d:
84 return "NAVD88"
85 if "ngvd" in d:
86 return "NGVD29"
87 if "gage height" in d:
88 return "gage height (local datum)"
89 return "unknown datum"
90
91
92class Series:
93 """A time series aligned to UTC epoch seconds, with provenance."""
94
95 def __init__(self, name, t, values, units, provenance, extra=None):
96 self.name = name
97 self.t = t # list[float] utc epoch seconds, ascending
98 self.values = values # list[float] (None for gaps already removed)
99 self.units = units
100 self.provenance = provenance # dict
101 self.extra = extra or {} # parallel dict of name -> list (e.g. wind dir)
102
103 def __len__(self):
104 return len(self.t)
105
106 def summary(self):
107 if not self.t:
108 return f"<{self.name}: EMPTY>"
109 return (f"<{self.name}: n={len(self.t)} "
110 f"{_iso(self.t[0])}..{_iso(self.t[-1])} "
111 f"units={self.units}>")
112
113
114def _iso(epoch):
115 return datetime.fromtimestamp(epoch, tz=timezone.utc).strftime("%Y-%m-%d %H:%M")
116
117
118# ---- observability (L1): per-request upstream-fetch counter (thread-local) ----
119_fetch_local = threading.local()
120
121
122def reset_fetch_count():
123 _fetch_local.n = 0
124
125
126def get_fetch_count():
127 return getattr(_fetch_local, "n", 0)
128
129
130def _get_with_retry(url, timeout):
131 """GET with bounded, jittered exponential backoff on transient failures
132 (connection/timeout/5xx). 4xx and the final attempt raise. No infinite retry."""
133 retries = int(os.environ.get("TIDES_HTTP_RETRIES", "3"))
134 backoff = float(os.environ.get("TIDES_HTTP_BACKOFF", "0.4"))
135 for attempt in range(retries):
136 try:
137 resp = requests.get(url, headers={"User-Agent": "tide-engine/1.0"},
138 timeout=timeout)
139 except requests.RequestException:
140 if attempt + 1 >= retries:
141 raise
142 time.sleep(backoff * (2 ** attempt) + random.uniform(0, backoff))
143 continue
144 if 500 <= resp.status_code < 600 and attempt + 1 < retries:
145 time.sleep(backoff * (2 ** attempt) + random.uniform(0, backoff))
146 continue
147 resp.raise_for_status() # 4xx, or the final 5xx -> raise
148 return resp
149 resp.raise_for_status()
150 return resp
151
152
153def _fetch_json(url, cache_key, ttl_seconds=None, timeout=90):
154 """GET url, caching the raw body to disk. Returns (parsed_json, cache_path).
155
156 ttl_seconds=None means cache never expires (historical data is immutable).
157 """
158 os.makedirs(CACHE_DIR, exist_ok=True)
159 h = hashlib.sha1(cache_key.encode()).hexdigest()[:16]
160 path = os.path.join(CACHE_DIR, f"{h}.json")
161 meta_path = path + ".meta"
162 if os.path.exists(path):
163 fresh = True
164 if ttl_seconds is not None:
165 if os.path.exists(meta_path):
166 with open(meta_path) as f:
167 fetched_at = float(f.read().strip())
168 fresh = (time.time() - fetched_at) < ttl_seconds
169 else:
170 fresh = False # expiring class without a timestamp = stale
171 if fresh:
172 with open(path, "r", encoding="utf-8") as f:
173 return json.load(f), path
174 _fetch_local.n = getattr(_fetch_local, "n", 0) + 1 # count real upstream hits
175 resp = _get_with_retry(url, timeout)
176 body = resp.text
177 data = json.loads(body)
178 # atomic write: a kill mid-write must leave the OLD file or the NEW complete
179 # file, never a truncated JSON (tempfile in the same dir + os.replace).
180 fd, tmp = tempfile.mkstemp(dir=CACHE_DIR, suffix=".tmp")
181 try:
182 with os.fdopen(fd, "w", encoding="utf-8") as f:
183 f.write(body)
184 os.replace(tmp, path)
185 except BaseException:
186 try:
187 os.unlink(tmp)
188 except OSError:
189 pass
190 raise
191 with open(meta_path, "w") as f:
192 f.write(str(time.time()))
193 return data, path
194
195
196def _parse_iso_to_epoch(s):
197 """USGS dateTime like '2026-06-17T17:15:00.000-04:00' -> utc epoch seconds."""
198 dt = datetime.fromisoformat(s)
199 if dt.tzinfo is None:
200 dt = dt.replace(tzinfo=timezone.utc)
201 return dt.timestamp()
202
203
204# --------------------------------------------------------------------------
205# USGS ground-truth water level
206# --------------------------------------------------------------------------
207def usgs_water_level(site, start_date, end_date, param="63160", ttl=None):
208 """Fetch USGS instantaneous water level.
209
210 param 63160 = water-surface elevation above NAVD88 (ft). Returns Series in
211 METERS (converted from feet) so the whole engine is metric.
212 Missing/flagged values (e.g. -999999) are dropped.
213 ttl=None on a NEAR-REAL-TIME window is upgraded to NRT_TTL_SECONDS (FIX 2) —
214 pinned historical/training windows keep the never-expire default.
215 """
216 if ttl is None and _is_nrt_end(end_date):
217 ttl = NRT_TTL_SECONDS
218 q = {
219 "format": "json",
220 "sites": site,
221 "parameterCd": param,
222 "startDT": start_date,
223 "endDT": end_date,
224 }
225 url = USGS_IV + "?" + urllib.parse.urlencode(q)
226 try:
227 data, path = _fetch_json(url, "usgs|" + url, ttl_seconds=ttl)
228 except requests.RequestException as exc:
229 # network/transport/HTTP error -> upstream unreachable, NOT no-data.
230 raise UpstreamError("USGS unreachable for %s/%s: %s" % (site, param, exc))
231 ts_list = data.get("value", {}).get("timeSeries", [])
232 if not ts_list:
233 raise NoDataError(f"USGS returned no timeSeries for {site}/{param}")
234 ts = ts_list[0]
235 si = ts["sourceInfo"]
236 vals = ts["values"][0]["value"]
237 t, v = [], []
238 for row in vals:
239 raw = row.get("value")
240 if raw is None or raw == "" or float(raw) <= -999998:
241 continue
242 t.append(_parse_iso_to_epoch(row["dateTime"]))
243 v.append(float(raw) * 0.3048) # ft -> m
244 param_desc = ts["variable"]["variableDescription"]
245 datum = _datum_from_desc(param_desc) # truthful — derived from THIS param
246 prov = {
247 "source": "USGS NWIS Instantaneous Values",
248 "site": site,
249 "site_name": si["siteName"],
250 "lat": float(si["geoLocation"]["geogLocation"]["latitude"]),
251 "lon": float(si["geoLocation"]["geogLocation"]["longitude"]),
252 "param": param,
253 "param_desc": param_desc,
254 "datum": datum,
255 "orig_units": "ft",
256 "url": url,
257 "cache": os.path.abspath(path),
258 "fetched_utc": _iso(time.time()),
259 }
260 return Series(f"usgs:{site}", t, v, "m (%s)" % datum, prov)
261
262
263# --------------------------------------------------------------------------
264# Wind + pressure forcing at exact location (ERA5 historical / NWS live)
265# --------------------------------------------------------------------------
266def _today_utc_date():
267 return datetime.now(tz=timezone.utc).date()
268
269
270def _win_bounds(start_date, end_date):
271 """[start 00:00 UTC, end 24:00 UTC) epoch bounds for an inclusive date range."""
272 s = datetime.fromisoformat(start_date).replace(tzinfo=timezone.utc).timestamp()
273 e = (datetime.fromisoformat(end_date).replace(tzinfo=timezone.utc)
274 + timedelta(days=1)).timestamp()
275 return s, e
276
277
278# ---- ERA5 reanalysis: read the processed point-cache written by the build-time
279# tool scripts/era5_fetch.py (NO cdsapi/netCDF at runtime) ---------------
280def era5_cache_path(lat, lon):
281 return os.path.join(CACHE_DIR, "era5_%.3f_%.3f.json" % (round(lat, 3), round(lon, 3)))
282
283
284def _era5_series(lat, lon, start_date, end_date):
285 p = era5_cache_path(lat, lon)
286 if not os.path.exists(p):
287 raise RuntimeError("ERA5 cache miss for %.3f,%.3f — run "
288 "scripts/era5_fetch.py fill" % (lat, lon))
289 with open(p, "r", encoding="utf-8") as f:
290 blob = json.load(f)
291 s, e = _win_bounds(start_date, end_date)
292 t, pres, ws, wd = [], [], [], []
293 for i, ti in enumerate(blob["t"]):
294 if s <= ti < e:
295 t.append(ti); pres.append(blob["pressure_hPa"][i])
296 ws.append(blob["wind_speed"][i]); wd.append(blob["wind_dir_from"][i])
297 if not t:
298 raise RuntimeError("ERA5 cache for %.3f,%.3f lacks %s..%s (have %s..%s)"
299 % (lat, lon, start_date, end_date,
300 _iso(blob["t"][0]), _iso(blob["t"][-1])))
301 prov = {"source": "ERA5 reanalysis (Copernicus CDS)",
302 "role": "historical/training", "lat": lat, "lon": lon,
303 "grid_lat": blob["provenance"].get("grid_lat"),
304 "grid_lon": blob["provenance"].get("grid_lon"),
305 "wind_units": "m/s", "pressure_units": "hPa",
306 "attribution": ERA5_ATTRIBUTION,
307 "licence": "Copernicus Licence (commercial use OK with attribution)",
308 "cache": os.path.abspath(p)}
309 return Series("met", t, pres, "hPa", prov,
310 extra={"wind_speed": ws, "wind_dir_from": wd})
311
312
313# ---- NWS live wind (US gov, public domain) via api.weather.gov -------------
314def _iso_dur_hours(dur):
315 """ISO8601 duration like 'PT1H','PT6H','P1DT2H' -> integer hours (>=1)."""
316 days = hours = 0
317 m = dur.replace("P", "")
318 if "T" in m:
319 d, t = m.split("T")
320 else:
321 d, t = m, ""
322 if d.endswith("D"):
323 days = int(d[:-1] or 0)
324 if t.endswith("H"):
325 hours = int(t[:-1] or 0)
326 return max(1, days * 24 + hours)
327
328
329def _nws_points(lat, lon):
330 url = "%s/points/%.4f,%.4f" % (NWS_API, lat, lon)
331 data, _ = _fetch_json(url, "nws|" + url, ttl_seconds=30 * 86400)
332 p = data["properties"]
333 return p["gridId"], p["gridX"], p["gridY"], p["observationStations"]
334
335
336def _expand_layer(values, factor):
337 """NWS gridpoint layer values (validTime interval + value) -> {epoch_hr: value}."""
338 out = {}
339 for v in values or []:
340 if v.get("value") is None:
341 continue
342 iso, dur = v["validTime"].split("/")
343 t0 = int(datetime.fromisoformat(iso).timestamp())
344 for h in range(_iso_dur_hours(dur)):
345 out[t0 + h * 3600] = float(v["value"]) * factor
346 return out
347
348
349def _nws_forecast_wind(gid, gx, gy):
350 url = "%s/gridpoints/%s/%d,%d" % (NWS_API, gid, gx, gy)
351 data, _ = _fetch_json(url, "nws|" + url, ttl_seconds=3600)
352 pr = data["properties"]
353 ws = _expand_layer(pr.get("windSpeed", {}).get("values"), 1 / 3.6) # km/h -> m/s
354 wd = _expand_layer(pr.get("windDirection", {}).get("values"), 1.0) # deg
355 return ws, wd
356
357
358def _nws_obs(stations_url, start_date, end_date):
359 """Recent hourly station obs -> {epoch_hr: (speed m/s, dir deg, slp hPa)}.
360
361 FIX 2 (2026-07-18): the observations endpoint caps at 500 features,
362 NEWEST-first — a single request silently returns only the last ~1.5 days
363 at a 5-min-cadence station, starving the wind model's 72 h trailing window
364 (the engine went blind with live gauges: 'no recent observation to anchor
365 to' was really 'no wind features at now'). Fetch BACKWARD in 1-day chunks
366 until the window is covered or history runs out."""
367 data, _ = _fetch_json(stations_url, "nws|" + stations_url, ttl_seconds=7 * 86400)
368 feats = data.get("features") or data.get("observationStations")
369 if not feats:
370 return {}
371 sid = feats[0]["properties"]["stationIdentifier"] if isinstance(feats[0], dict) \
372 else feats[0].rsplit("/", 1)[-1]
373 s, e = _win_bounds(start_date, end_date)
374 e = min(e, time.time())
375 out = {}
376
377 def _mps(q):
378 if not q or q.get("value") is None:
379 return None
380 val = float(q["value"])
381 return val / 3.6 if "km_h" in (q.get("unitCode") or "") else val
382
383 chunk = 86400 # 1 day < 500 obs at any reporting cadence
384 e_cur = e
385 while e_cur > s:
386 s_cur = max(s, e_cur - chunk)
387 url = ("%s/stations/%s/observations?start=%s&end=%s" %
388 (NWS_API, sid,
389 datetime.fromtimestamp(s_cur, timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ"),
390 datetime.fromtimestamp(e_cur, timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ")))
391 try:
392 data, _ = _fetch_json(url, "nws|" + url, ttl_seconds=3600)
393 except Exception:
394 break # transport blip: keep what we have
395 feats = data.get("features", [])
396 if not feats:
397 break # history floor reached (NWS keeps ~7 days)
398 for f in feats:
399 p = f["properties"]
400 try:
401 t = int(datetime.fromisoformat(p["timestamp"].replace("Z", "+00:00")).timestamp())
402 except Exception:
403 continue
404 thr = t - (t % 3600)
405 sp = _mps(p.get("windSpeed"))
406 di = p.get("windDirection", {}).get("value")
407 slpq = p.get("seaLevelPressure") or p.get("barometricPressure") or {}
408 slp = None if slpq.get("value") is None else float(slpq["value"]) / 100.0
409 if sp is not None and di is not None:
410 out[thr] = (sp, float(di), slp)
411 e_cur = s_cur
412 return out
413
414
415def _live_series(lat, lon, start_date, end_date):
416 """Wind+pressure for a window reaching near-now/future: NWS obs (recent past)
417 + NWS gridpoint forecast (future), with ERA5 for the older past if cached.
418 NWS has no forecast surface pressure -> pressure is held at the last observed
419 value (documented degradation; cancels in the woff(t)-woff(now) delta)."""
420 s, e = _win_bounds(start_date, end_date)
421 grid = [int(s) + k * 3600 for k in range(int((e - s) // 3600))]
422 gid, gx, gy, stations_url = _nws_points(lat, lon)
423 fc_ws, fc_wd = _nws_forecast_wind(gid, gx, gy)
424 obs = _nws_obs(stations_url, start_date, end_date)
425 # older past from ERA5 cache (best-effort; extends high-pass support)
426 era = {}
427 try:
428 cut = (_today_utc_date() - timedelta(days=6)).isoformat()
429 es = _era5_series(lat, lon, start_date, cut)
430 era = {int(es.t[i]): (es.extra["wind_speed"][i], es.extra["wind_dir_from"][i],
431 es.values[i]) for i in range(len(es.t))}
432 except Exception:
433 pass
434 last_slp = next((v[2] for v in reversed(list(obs.values())) if v[2] is not None), 1013.25)
435 t, pres, ws, wd = [], [], [], []
436 used = set()
437 for g in grid:
438 sp = di = pp = None
439 if g in obs:
440 sp, di, pp = obs[g]; used.add("nws-obs")
441 elif g in era:
442 sp, di, pp = era[g]; used.add("era5")
443 elif g in fc_ws and g in fc_wd:
444 sp, di = fc_ws[g], fc_wd[g]; pp = last_slp; used.add("nws-forecast")
445 if sp is None or di is None:
446 continue
447 t.append(float(g)); ws.append(sp); wd.append(di)
448 pres.append(pp if pp is not None else last_slp)
449 if not t:
450 raise RuntimeError("no live wind for %.3f,%.3f %s..%s" % (lat, lon, start_date, end_date))
451 prov = {"source": "NWS api.weather.gov (gridpoint forecast + station obs)",
452 "role": "live/forecast", "lat": lat, "lon": lon,
453 "components": sorted(used),
454 "wind_units": "m/s", "pressure_units": "hPa",
455 "pressure_note": "NWS provides no forecast surface pressure; held at "
456 "last observed value (delta-cancels in the product)",
457 "licence": "US Government work — public domain",
458 "fetched_utc": _iso(time.time())}
459 if "era5" in used:
460 prov["attribution"] = ERA5_ATTRIBUTION
461 return Series("met", t, pres, "hPa", prov,
462 extra={"wind_speed": ws, "wind_dir_from": wd})
463
464
465def met_forcing(lat, lon, start_date, end_date, ttl=None):
466 """Hourly wind (m/s), wind dir (deg FROM), surface pressure (hPa) at (lat,lon).
467 Routes by window age: windows ending >= LIVE_AGE_DAYS days ago come from ERA5
468 reanalysis (training / backtest / fixed-epoch replay); windows reaching
469 near-now/future come from NWS (live). Series .values is pressure; .extra holds
470 wind. Contract is identical to the previous Open-Meteo implementation."""
471 end_d = datetime.fromisoformat(end_date).date()
472 if (_today_utc_date() - end_d).days >= LIVE_AGE_DAYS:
473 return _era5_series(lat, lon, start_date, end_date)
474 return _live_series(lat, lon, start_date, end_date)
475
476
477if __name__ == "__main__":
478 wl = usgs_water_level("255327081275900", "2026-06-10", "2026-06-17")
479 print(wl.summary())
480 print(" ground-truth provenance:", wl.provenance["site_name"],
481 wl.provenance["param_desc"])
482 met = met_forcing(wl.provenance["lat"], wl.provenance["lon"],
483 "2026-06-10", "2026-06-17")
484 print(met.summary())
485 print(" met source:", met.provenance["source"])
486