#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Dalga Koşulları ve Darboğazlarda Rüzgâr Ekstremleri: 10 Havzada Günlük Dalga Yüksekliği (Open-Meteo Marine / CMEMS) ve NASA MERRA-2 ile 30 Yıllık Gumbel Dönüş Seviyeleri, Fırtına Günleri (1996–2026)
Wave Conditions & Wind Extremes at Chokepoints: Daily Wave Height in 10 Basins (Open-Meteo Marine / CMEMS) and 30-Year Gumbel Return Levels & Gale Days from NASA MERRA-2 (1996–2026)
==============================================================================
GOMD Python SDK — Gerçek veriyle yeniden üretilebilir vaka çalışması / Reproducible real-data case study
Google Colab Install: !pip install -q --upgrade gomd

Veri kaynakları / Data sources (hepsi GOMD sunucusundan canlı çekilir / all fetched live from the GOMD server):
  - metocean_global_basin_wave_height_max_m, _swell_height_max_m, _wind_wave_height_max_m, _sst_mean_c — Open-Meteo Marine API / CMEMS küresel analiz (son 93 gün)
  - nasa_merra2_daily_wind_10m_max_ms — NASA POWER / MERRA-2 günlük en yüksek 10 m rüzgâr (1996–2025)

Not / Note: Bu betik sentetik veri üretmez. Bir seri sunucuda yoksa betik açık bir hata mesajıyla durur.
            This script never generates synthetic data. If a series is missing on the server, it stops with an explicit error.
"""

# ==============================================================================
# ADIM 0: GOMD Python SDK Kurulumu ve Canlı Veri Yardımcıları
# STEP 0: Install the GOMD Python SDK & Live Data Helpers
# ==============================================================================
import importlib
import os
import subprocess
import sys
import math
import json
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

if hasattr(sys.stdout, "reconfigure"):
    try:
        sys.stdout.reconfigure(encoding="utf-8", errors="replace")
    except Exception:
        pass

_IN_COLAB = "google.colab" in sys.modules
_GOMD_WHEEL_URL = "https://gomd.pirireis.edu.tr/sdk/dist/gomd-2.0.1-py3-none-any.whl"
try:
    import gomd  # noqa: F401
    if _IN_COLAB:
        subprocess.check_call([sys.executable, "-m", "pip", "install", "-q", "--force-reinstall", "--no-cache-dir", _GOMD_WHEEL_URL])
except ImportError:
    subprocess.check_call([sys.executable, "-m", "pip", "install", "-q", "--force-reinstall", "--no-cache-dir", _GOMD_WHEEL_URL])

import gomd
for _mod_name in list(sys.modules.keys()):
    if _mod_name == "gomd" or _mod_name.startswith("gomd."):
        importlib.reload(sys.modules[_mod_name])
import gomd
from gomd import GOMDClient

client = GOMDClient("https://gomd.pirireis.edu.tr")


def fetch_gomd_records(indicator: str, entity: str, period_from=None, period_to=None, limit: int = 50000) -> pd.DataFrame:
    """
    GOMD sunucusundan bir göstergenin TEK bir varlık (ülke/havza) için ham kayıtlarını çeker.
    Yalnız istenen varlığın kayıtları tutulur; sunucu başka bir varlık döndürürse kullanılmaz.
    Returns a DataFrame with columns: period (str), value (float).
    """
    rows = client.get_series(indicator=indicator, entity=entity, period_from=period_from,
                             period_to=period_to, limit=limit, as_df=False)
    keep = [r for r in rows if str(r.get("entity_code") or "") == entity and r.get("value") is not None]
    if not keep:
        raise RuntimeError(f"GOMD sunucusunda veri yok: indicator={indicator}, entity={entity}")
    df = pd.DataFrame({"period": [str(r["period"]) for r in keep], "value": [float(r["value"]) for r in keep]})
    df = df.drop_duplicates("period").sort_values("period").reset_index(drop=True)
    print(f"📡 [GOMD] {indicator:<44} | {entity:<18} | N={len(df):<6} ({df['period'].iloc[0][:10]} – {df['period'].iloc[-1][:10]})")
    return df


def fetch_gomd_annual_series(indicator: str, entity: str, period_from: int = 1900, period_to: int = 2100,
                             how: str = "mean", min_count: int = 1) -> dict:
    """
    Yıllık seri döndürür: {yıl: değer}. Kayıtlar yıllık değilse (aylık/günlük) yıl içi ortalama alınır
    ve yalnız en az `min_count` kaydı olan yıllar tutulur. Returns {year: value}.
    """
    df = fetch_gomd_records(indicator, entity)
    df["year"] = df["period"].str[:4].astype(int)
    df = df[(df["year"] >= period_from) & (df["year"] <= period_to)]
    g = df.groupby("year")["value"]
    agg = g.mean() if how == "mean" else g.last()
    cnt = g.count()
    out = {int(y): float(v) for y, v in agg.items() if cnt[y] >= min_count}
    if not out:
        raise RuntimeError(f"Seçilen dönemde veri yok: indicator={indicator}, entity={entity}")
    return dict(sorted(out.items()))


def fmt_p(p: float) -> str:
    """p değerini okunur biçimde yazar / readable p-value."""
    return "<1e-6" if p < 1e-6 else f"{p:.4f}"


print(f"✅ GOMD Python SDK Sürümü / Version : v{gomd.__version__}")
print(f"🌐 Canlı Veri Sunucusu / Server     : {client.base_url}")
print(f"🏛️ Kurum / Institution              : Pîrî Reis Üniversitesi & T.C. Çevre, Şehircilik ve İklim Değişikliği Bakanlığı Çevre Yönetimi Genel Müdürlüğü — GOMD (COP31 Türkiye)")

# ==============================================================================
# ADIM 1: Son 93 günün dalga koşulları — 10 havza ve darboğaz (Open-Meteo Marine / CMEMS modeli, günlük)
# STEP 1: Wave conditions over the last 93 days — 10 basins & chokepoints (Open-Meteo Marine / CMEMS model, daily)
# ==============================================================================
from gomd.stats import extreme_return_levels, mann_kendall_sen, panel_fixed_effects

BASINS = {
    "OCN_MED_LEVANTINE": "Levant (Doğu Akdeniz)", "OCN_MED_AEGEAN": "Ege", "OCN_MARMARA_STRAITS": "Marmara ve Boğazlar",
    "OCN_BLACKSEA_CENTRAL": "Orta Karadeniz", "OCN_MED_IONIAN": "İyon", "OCN_MED_WEST": "Batı Akdeniz",
    "OCN_CHOKE_SUEZ_REDSEA": "Süveyş / Kızıldeniz", "OCN_CHOKE_GIBRALTAR": "Cebelitarık", "OCN_CHOKE_HORMUZ": "Hürmüz",
    "OCN_CHOKE_BABELMANDEB": "Babülmendep",
}
VARS = {"hs": "metocean_global_basin_wave_height_max_m", "swell": "metocean_global_basin_swell_height_max_m",
        "ww": "metocean_global_basin_wind_wave_height_max_m", "sst": "metocean_global_basin_sst_mean_c"}
daily = {}
for code in BASINS:
    parts = []
    for v, ind in VARS.items():
        d = fetch_gomd_records(ind, code)
        parts.append(pd.Series(d["value"].values, index=pd.to_datetime(d["period"].str[:10]), name=v))
    daily[code] = pd.concat(parts, axis=1).dropna()

rows, fe_rows = [], []
for code, name in BASINS.items():
    d = daily[code]
    rows.append({"Havza": name, "Gün": len(d), "Dönem": f"{d.index.min():%Y-%m-%d} – {d.index.max():%Y-%m-%d}",
                 "Ort. SST (°C)": round(d["sst"].mean(), 1), "Hs medyan (m)": round(d["hs"].median(), 2),
                 "Hs %95 (m)": round(d["hs"].quantile(0.95), 2), "Hs en yüksek (m)": round(d["hs"].max(), 2),
                 "Hs ≥ 2,5 m gün": int((d["hs"] >= 2.5).sum()), "Soluğan payı (medyan)": round((d["swell"] / d["hs"]).median(), 2)})
    fe_rows += [{"basin": code, "hs": r.hs, "swell": r.swell, "ww": r.ww} for r in d.itertuples()]
df_w = pd.DataFrame(rows)
print("\n🌊 Son 93 günün dalga koşulları (günlük en yüksek Hs; model analizi):")
print(df_w.to_string(index=False))
fe_s = panel_fixed_effects(fe_rows, "basin", "hs", "swell")
fe_w = panel_fixed_effects(fe_rows, "basin", "hs", "ww")
print(f"📐 Havza içi ilişki (N = {fe_s['n_obs']} gün): Hs ~ soluğan β = {fe_s['beta_fe']:+.2f} (R² = {fe_s['r_squared_within']:.2f}); "
      f"Hs ~ rüzgâr dalgası β = {fe_w['beta_fe']:+.2f} (R² = {fe_w['r_squared_within']:.2f})")
print("ℹ️  93 günlük tek mevsim verisi, dönüş seviyesi (ör. 100 yıllık dalga) hesaplamaya yetmez; ekstremler ADIM 2'de 30 yıllık seriyle hesaplanır.")

# ==============================================================================
# ADIM 2: 30 yıllık rüzgâr ekstremleri — NASA MERRA-2 günlük en yüksek 10 m rüzgârı (1996–2025)
# STEP 2: 30-year wind extremes — NASA MERRA-2 daily maximum 10 m wind (1996–2025)
# ==============================================================================
SITES = {
    "STRAIT_BOSPHORUS": "İstanbul Boğazı", "STRAIT_DARDANELLES": "Çanakkale Boğazı", "BASIN_MARMARA_SEA": "Marmara Denizi",
    "BASIN_AEGEAN_SEA": "Ege Denizi", "BASIN_BLACK_SEA": "Karadeniz", "BASIN_MED_LEVANTINE": "Levant", "PORT_MERSIN": "Mersin Limanı",
    "PORT_PIRAEUS": "Pire Limanı", "CANAL_SUEZ": "Süveyş Kanalı", "STRAIT_BAB_EL_MANDEB": "Babülmendep", "STRAIT_GIBRALTAR": "Cebelitarık",
    "STRAIT_HORMUZ": "Hürmüz",
}
GALE = 17.2   # m/s, Beaufort 8 (fırtına) alt sınırı
wind, ext_rows, gale = {}, [], {}
for code, name in SITES.items():
    d = fetch_gomd_records("nasa_merra2_daily_wind_10m_max_ms", code)
    s = pd.Series(d["value"].values, index=pd.to_datetime(d["period"].str[:10]))
    cnt = s.groupby(s.index.year).count()
    full = cnt[cnt >= 360].index                                  # yalnız tam yıllar
    amax = s.groupby(s.index.year).max().loc[full]
    gale[code] = (s >= GALE).groupby(s.index.year).sum().loc[full]
    wind[code] = amax
    ev = extreme_return_levels(amax.values.tolist(), return_periods=(10, 50, 100))
    mk = mann_kendall_sen(amax.values.tolist(), t=amax.index.tolist())
    ext_rows.append({"Yer": name, "Yıl": len(amax), "Yıllık maks. ort. (m/s)": round(amax.mean(), 1),
                     "Gözlenen en yüksek (m/s)": round(amax.max(), 1), "10 yıl (m/s)": round(ev["return_levels"]["return_level_10yr"], 1),
                     "50 yıl (m/s)": round(ev["return_levels"]["return_level_50yr"], 1), "100 yıl (m/s)": round(ev["return_levels"]["return_level_100yr"], 1),
                     "Fırtına günü / yıl": round(gale[code].mean(), 1), "Yıllık maks. eğimi (m/s/10 yıl)": round(mk["theil_sen_slope"] * 10, 2),
                     "p": fmt_p(mk["p_value"])})
df_ext = pd.DataFrame(ext_rows).sort_values("100 yıl (m/s)", ascending=False)
print("\n💨 Rüzgâr ekstremleri, Gumbel (yıllık maksimumlar, momentler yöntemi):")
print(df_ext.to_string(index=False))
print("ℹ️  MERRA-2 yaklaşık 50 km çözünürlüklü yeniden analizdir; dar boğazlarda yerel rüzgârı düşük tahmin edebilir. Tasarım için yerel ölçümle doğrulanmalıdır.")

# ==============================================================================
# ADIM 3: 4 Panelli Görselleştirme / 4-panel figure
# ==============================================================================
fig, axes = plt.subplots(2, 2, figsize=(14, 9.5), dpi=120)
fig.patch.set_facecolor("#071324")
for ax in axes.flat:
    ax.set_facecolor("#0b1d36"); ax.tick_params(colors="#cbd5e1")
    for sp in ax.spines.values():
        sp.set_color("#1e3a5f")
    ax.grid(alpha=0.18, color="#64748b")
leg = dict(facecolor="#0f172a", edgecolor="#334155", labelcolor="#f8fafc", fontsize=8)

ax1 = axes[0, 0]
lev = daily["OCN_MED_AEGEAN"]
ax1.plot(lev.index, lev["hs"], color="#38bdf8", lw=1.8, label="Toplam Hs")
ax1.plot(lev.index, lev["swell"], color="#22c55e", lw=1.3, ls="--", label="Soluğan")
ax1.plot(lev.index, lev["ww"], color="#f59e0b", lw=1.3, ls=":", label="Rüzgâr dalgası")
ax1.set_ylabel("m (günlük en yüksek)", color="#cbd5e1")
ax1.set_title("A) Ege: son 93 günün dalga yüksekliği (model)", color="#f8fafc", fontweight="bold"); ax1.legend(**leg)
ax1.tick_params(axis="x", labelrotation=20)

ax2 = axes[0, 1]
order = df_w.sort_values("Hs %95 (m)")["Havza"].tolist()
name_to_code = {v: k for k, v in BASINS.items()}
_bx = dict(patch_artist=True, boxprops=dict(facecolor="#1e3a5f", color="#38bdf8"), medianprops=dict(color="#f59e0b"),
           whiskerprops=dict(color="#94a3b8"), capprops=dict(color="#94a3b8"), flierprops=dict(markeredgecolor="#94a3b8", markersize=3))
_data = [daily[name_to_code[n]]["hs"].values for n in order]
try:
    ax2.boxplot(_data, orientation="horizontal", **_bx)      # matplotlib ≥ 3.10
except TypeError:
    ax2.boxplot(_data, vert=False, **_bx)                    # eski matplotlib sürümleri
ax2.set_yticks(range(1, len(order) + 1)); ax2.set_yticklabels(order, fontsize=8)   # matplotlib sürümünden bağımsız etiketleme
ax2.set_xlabel("Hs (m)", color="#cbd5e1")
ax2.set_title("B) Havzalara göre günlük en yüksek Hs dağılımı (son 93 gün)", color="#f8fafc", fontweight="bold")

ax3 = axes[1, 0]
de = df_ext.sort_values("100 yıl (m/s)")
y = np.arange(len(de))
ax3.barh(y - 0.2, de["Gözlenen en yüksek (m/s)"], height=0.38, color="#38bdf8", label="1996–2025 gözlenen en yüksek")
ax3.barh(y + 0.2, de["100 yıl (m/s)"], height=0.38, color="#f43f5e", label="100 yıllık dönüş seviyesi (Gumbel)")
ax3.axvline(GALE, color="#e2e8f0", ls=":", lw=1.2, label="Beaufort 8 (17,2 m/s)")
ax3.set_yticks(y); ax3.set_yticklabels(de["Yer"], fontsize=8)
ax3.set_xlabel("günlük en yüksek 10 m rüzgâr (m/s)", color="#cbd5e1")
ax3.set_title("C) Rüzgâr ekstremleri — MERRA-2, 30 yıl", color="#f8fafc", fontweight="bold"); ax3.legend(loc="lower right", **leg)

ax4 = axes[1, 1]
for code, c in [("BASIN_AEGEAN_SEA", "#38bdf8"), ("BASIN_BLACK_SEA", "#22c55e"), ("STRAIT_DARDANELLES", "#f59e0b"), ("BASIN_MARMARA_SEA", "#f43f5e"), ("PORT_MERSIN", "#a78bfa")]:
    g = gale[code]; ax4.plot(g.index, g.values, color=c, lw=1.6, marker="o", ms=3, label=SITES[code])
ax4.set_ylabel("gün / yıl (rüzgâr ≥ 17,2 m/s)", color="#cbd5e1")
ax4.set_title("D) Yıllık fırtına günü sayısı — MERRA-2", color="#f8fafc", fontweight="bold"); ax4.legend(**leg)

fig.suptitle("GOMD — Vaka Çalışması 8: Dalga koşulları ve darboğazlarda rüzgâr ekstremleri (gerçek veri)", color="#ffffff", fontsize=12.5, fontweight="bold")
plt.tight_layout()
fig.savefig("cs08_high_frequency_metocean_wave_extremes_chokepoints.png", dpi=140, bbox_inches="tight", facecolor=fig.get_facecolor())
plt.show()