| """ClearCast Space: EXPERIMENTAL underwater visibility estimate for Finistère (satellite Secchi depth -> diver visibility).""" |
| import joblib, numpy as np, pandas as pd, requests, gradio as gr |
| import plotly.graph_objects as go |
| from features import daily_weather_from_hourly |
| from spots import SPOTS |
|
|
| MODELS = joblib.load("zsd_models.joblib") |
| LATEST = pd.read_csv("latest.csv").set_index("spot") |
| HIND = pd.read_parquet("hindcast.parquet") |
| NAMES = list(SPOTS) |
| K, LINK = 0.76, {"0.1": -0.731, "0.5": 0.008, "0.9": 0.597} |
| SD_LINK = (LINK["0.9"] - LINK["0.1"]) / 2.563 |
| LABEL = ("EXPERIMENTAL estimate. Visibility ≈ 0.76 × satellite Secchi depth, a relation checked against ~6,000 diver surveys " |
| "(Reef Life Survey) worldwide, NOT in Brittany. No validation against Finistère visibility observations exists yet.") |
|
|
|
|
| def live_weather(lat, lon): |
| base = dict(latitude=lat, longitude=lon, past_days=10, forecast_days=1, timezone="UTC") |
| a = requests.get("https://api.open-meteo.com/v1/forecast", params=dict(base, hourly="wind_speed_10m,wind_direction_10m,precipitation"), timeout=30).json() |
| m = requests.get("https://marine-api.open-meteo.com/v1/marine", params=dict(base, hourly="wave_height,wave_period,wave_direction"), timeout=30).json() |
| w = pd.DataFrame(a["hourly"]).merge(pd.DataFrame(m["hourly"]), on="time") |
| w["time"] = pd.to_datetime(w["time"]) |
| return w[w["time"] <= pd.Timestamp.utcnow().tz_localize(None)] |
|
|
|
|
| def predict(spot, horizon, last_zsd): |
| lat, lon, _ = SPOTS[spot] |
| h = int(horizon) |
| d = daily_weather_from_hourly(live_weather(lat, lon)) |
| today = d.iloc[-1:].copy() |
| today["last_kd"] = np.log(float(last_zsd)) |
| today["last_age"] = 0 |
| doy = pd.Timestamp.utcnow().dayofyear |
| today["doy_sin"], today["doy_cos"] = np.sin(2 * np.pi * doy / 365.25), np.cos(2 * np.pi * doy / 365.25) |
| today["lat"], today["lon"] = lat, lon |
| m = MODELS[h] |
| X = today[m["features"]] |
| mg = m["conformal_margin"] |
| l10, l50, l90 = (float(m["q"][k].predict(X)[0]) for k in ("0.1", "0.5", "0.9")) |
| l10, l90 = l10 - mg, l90 + mg |
| sd = np.sqrt(((l90 - l10) / 2.563) ** 2 + SD_LINK ** 2) |
| mu = l50 + np.log(K) + LINK["0.5"] |
| v10, v50, v90 = (float(np.exp(mu + z * sd)) for z in (-1.2816, 0, 1.2816)) |
| c = d.iloc[-1][["hs_max", "tp_mean", "wind", "rain_r3"]] |
| return (f"### {spot}: estimated visibility in +{h} day(s)\n**≈ {v50:.0f} m** (80% range {v10:.0f} – {v90:.0f} m)\n\n" |
| f"Forecast satellite Secchi depth ≈ {np.exp(l50):.1f} m. Today: Hs max {c.hs_max:.1f} m, period {c.tp_mean:.1f} s, " |
| f"wind {c.wind:.1f} km/h, rain 3-day mean {c.rain_r3:.1f} mm.\n\n" |
| f"Starting Secchi depth used: {float(last_zsd):.1f} m. Snapshot: {LATEST.loc[spot,'last_zsd']:.1f} m from **{LATEST.loc[spot,'last_date']}** " |
| f"(stale: enter a recent value if you have one).\n\n> {LABEL}") |
|
|
|
|
| def hindcast_plot(spot): |
| d = HIND[HIND.spot == spot].sort_values("date") |
| f = go.Figure() |
| f.add_traces([go.Scatter(x=d.date, y=d.p90, line=dict(width=0), showlegend=False), |
| go.Scatter(x=d.date, y=d.p10, fill="tonexty", line=dict(width=0), name="80% range", fillcolor="rgba(70,130,180,0.25)"), |
| go.Scatter(x=d.date, y=d.p50, name="Estimated visibility (24 h ahead)", line=dict(color="steelblue")), |
| go.Scatter(x=d.date, y=d.zsd_equiv_vis, mode="markers", name="0.76 × satellite Secchi (not a diver measurement)", marker=dict(color="black", size=4))]) |
| f.update_layout(title=f"2026 test period: {spot}", yaxis_title="Visibility (m)", height=380, margin=dict(l=10, r=10, t=40, b=10)) |
| return f |
|
|
|
|
| with gr.Blocks(title="ClearCast") as demo: |
| gr.Markdown("# ClearCast — Finistère underwater visibility (experimental)\n" + LABEL + |
| "\n\nPipeline: LightGBM quantile forecast of 4 km satellite Secchi depth (Copernicus Marine) from waves, wind and rain (Open-Meteo), " |
| "then a physical scaling to visibility with uncertainty from the diver-survey calibration. Skill vs persistence is small (see model card).") |
| with gr.Row(): |
| spot = gr.Dropdown(NAMES, value="Brest Iroise", label="Spot") |
| hz = gr.Radio([1, 2, 3], value=1, label="Horizon (days)") |
| lk = gr.Number(value=round(float(LATEST.loc["Brest Iroise", "median_zsd"]), 1), label="Latest satellite Secchi depth (m)") |
| btn = gr.Button("Estimate visibility", variant="primary") |
| out = gr.Markdown() |
| plot = gr.Plot() |
| btn.click(predict, [spot, hz, lk], out) |
| spot.change(lambda s: (hindcast_plot(s), round(float(LATEST.loc[s, "median_zsd"]), 1)), spot, [plot, lk]) |
| demo.load(hindcast_plot, spot, plot) |
|
|
| if __name__ == "__main__": |
| demo.launch() |
|
|
