#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Doğu Akdeniz, Ege ve Karadeniz Okyanus-İklim Veri Madenciliği: Günlük Uydu SST, Hobday Deniz Isı Dalgaları, CMEMS Akdeniz Isı İçeriği ve SOCAT fCO₂ (1982–2025)
Eastern Mediterranean, Aegean & Black Sea Ocean-Climate Data Mining: Daily Satellite SST, Hobday Marine Heatwaves, CMEMS Mediterranean Heat Content & SOCAT fCO₂ (1982–2025)
==============================================================================
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):
  - noaa_oisst_area_mean_daily_v1 — NOAA OISST v2.1 günlük alan ortalaması (SST_EMED_BOX, SST_AEGEAN_BOX, SST_BLACKSEA_BOX, SST_WMED_BOX)
  - omi_ohc_anomaly — Copernicus Marine OMI, Akdeniz 0–700 m ısı içeriği anomalisi (SEA_MEDITERRANEAN)
  - socat_fco2_surface, socat_salinity_surface — SOCAT gridli yıllık yüzey fCO₂ ve tuzluluk (SEA_MEDITERRANEAN)

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: Günlük Uydu SST (NOAA OISST v2.1) — Doğu Akdeniz, Ege, Karadeniz, Batı Akdeniz
# STEP 1: Daily satellite SST (NOAA OISST v2.1) for four Mediterranean & Black Sea boxes
# ==============================================================================
from gomd.stats import (
    mann_kendall_sen,
    pettitt_changepoint,
    transfer_entropy_causal,
    multivariate_data_mining_pipeline,
)

BOXES = {
    "SST_EMED_BOX": "Doğu Akdeniz / Eastern Med",
    "SST_AEGEAN_BOX": "Ege / Aegean",
    "SST_BLACKSEA_BOX": "Karadeniz / Black Sea",
    "SST_WMED_BOX": "Batı Akdeniz / Western Med",
}
daily = {}
for code in BOXES:
    d = fetch_gomd_records("noaa_oisst_area_mean_daily_v1", code)
    d["date"] = pd.to_datetime(d["period"].str[:10])
    daily[code] = d.set_index("date")["value"].asfreq("D")

# Yalnız tam yıllar (en az 360 gün) yıllık ortalamaya girer / only complete years (>= 360 days)
annual_sst = {}
for code, s in daily.items():
    g = s.groupby(s.index.year)
    annual_sst[code] = {int(y): float(v) for y, v in g.mean().items() if g.count()[y] >= 360}


def hobday_mhw_days(s: pd.Series, clim_start: int = 1982, clim_end: int = 2011, pct: float = 90.0,
                    window: int = 11, min_len: int = 5, max_gap: int = 2) -> dict:
    """
    Hobday vd. (2016) deniz ısı dalgası tanımı: güne özgü 90. persentil eşiği (1982–2011, 11 günlük pencere),
    en az 5 gün süren ve aralarında en çok 2 gün boşluk olan olaylar. Yıl başına MHW gün sayısını döndürür.
    """
    base = s[(s.index.year >= clim_start) & (s.index.year <= clim_end)].dropna()
    doy = base.index.day_of_year.values
    vals = base.values
    thresh = np.empty(366)
    half = window // 2
    for dd in range(1, 367):
        dist = np.abs(doy - dd)
        dist = np.minimum(dist, 366 - dist)
        thresh[dd - 1] = np.percentile(vals[dist <= half], pct)
    x = s.dropna()
    above = x.values > thresh[x.index.day_of_year.values - 1]
    # olaylari bul / find events
    events, start = [], None
    for i, a in enumerate(above):
        if a and start is None:
            start = i
        if (not a) and start is not None:
            events.append((start, i - 1)); start = None
    if start is not None:
        events.append((start, len(above) - 1))
    events = [e for e in events if e[1] - e[0] + 1 >= min_len]
    merged = []
    for e in events:
        if merged and (x.index[e[0]] - x.index[merged[-1][1]]).days - 1 <= max_gap:
            merged[-1] = (merged[-1][0], e[1])
        else:
            merged.append(e)
    flag = pd.Series(False, index=x.index)
    for a, b in merged:
        flag.iloc[a:b + 1] = True
    per_year = flag.groupby(flag.index.year).sum()
    days_in_year = flag.groupby(flag.index.year).count()
    return {int(y): int(v) for y, v in per_year.items() if days_in_year[y] >= 360}


mhw_days = {code: hobday_mhw_days(daily[code]) for code in ("SST_EMED_BOX", "SST_AEGEAN_BOX")}

# ==============================================================================
# ADIM 2: Akdeniz Isı İçeriği (CMEMS OMI) ve Yüzey fCO₂ (SOCAT) — yıllık
# STEP 2: Mediterranean ocean heat content (CMEMS OMI) & surface fCO₂ (SOCAT), annual
# ==============================================================================
ohc_med = fetch_gomd_annual_series("omi_ohc_anomaly", "SEA_MEDITERRANEAN")            # J m-2 (0–700 m anomaly)
fco2_med = fetch_gomd_annual_series("socat_fco2_surface", "SEA_MEDITERRANEAN")       # µatm
sss_med = fetch_gomd_annual_series("socat_salinity_surface", "SEA_MEDITERRANEAN")    # PSU, shipboard

# ==============================================================================
# ADIM 3: Trend, Kırılma ve Çok Değişkenli Veri Madenciliği
# STEP 3: Trends, change points and multivariate data mining
# ==============================================================================
trend_rows = []
for code, name in BOXES.items():
    yrs = sorted(annual_sst[code]); vals = [annual_sst[code][y] for y in yrs]
    mk = mann_kendall_sen(vals, t=yrs)
    trend_rows.append({"Bölge / Region": name, "Yıllar": f"{yrs[0]}–{yrs[-1]}",
                       "Ortalama SST (°C)": round(float(np.mean(vals)), 2),
                       "Sen eğimi (°C/10 yıl)": round(mk["theil_sen_slope"] * 10, 3),
                       "p": fmt_p(mk['p_value']), "Trend": mk["trend_classification"]})
df_trend = pd.DataFrame(trend_rows)
print("\n🌡️ Yıllık ortalama SST trendleri (NOAA OISST v2.1, tam yıllar):")
print(df_trend.to_string(index=False))

emed_years = sorted(annual_sst["SST_EMED_BOX"])
emed_vals = [annual_sst["SST_EMED_BOX"][y] for y in emed_years]
mk_emed = mann_kendall_sen(emed_vals, t=emed_years)
pt_emed = pettitt_changepoint(emed_vals, labels=emed_years)

for code in mhw_days:
    m = mhw_days[code]
    early = np.mean([m[y] for y in m if 1982 <= y <= 1991])
    late = np.mean([m[y] for y in m if y >= max(m) - 9])
    print(f"🔥 {BOXES[code]:<28}: MHW gün/yıl 1982–1991 ort. = {early:5.1f} | son 10 yıl ort. = {late:5.1f}")

# Ortak yıllar: tüm değişkenlerin gözlendiği yıllar / years where every variable is observed
features = {
    "sst_emed_c": annual_sst["SST_EMED_BOX"],
    "sst_aegean_c": annual_sst["SST_AEGEAN_BOX"],
    "sst_blacksea_c": annual_sst["SST_BLACKSEA_BOX"],
    "mhw_days_emed": mhw_days["SST_EMED_BOX"],
    "ohc_med_jm2": ohc_med,
    "fco2_med_uatm": fco2_med,
}
common = sorted(set.intersection(*[set(v) for v in features.values()]))
df = pd.DataFrame({"year": common, **{k: [v[y] for y in common] for k, v in features.items()}})
print(f"\n📊 Ortak yıllar matrisi (N={len(df)} yıl, {common[0]}–{common[-1]}):")
print(df.tail(8).to_string(index=False))

cols = list(features)
mining = multivariate_data_mining_pipeline(
    feature_matrix=df[cols].values.tolist(), feature_names=cols, period_labels=df["year"].tolist(),
    target_index=0, k_clusters=3, ridge_alphas=(0.01, 0.1, 1.0, 5.0, 10.0),
)
pca_info, km_info, ml_info = mining["pca"], mining["kmeans_regimes"], mining["ridge_ml_forecast"]
df["PC1"] = [s["pc1"] for s in pca_info["scores"]]
df["PC2"] = [s["pc2"] for s in pca_info["scores"]]
df["regime_id"] = km_info["regime_labels"]
te = transfer_entropy_causal(df["sst_emed_c"].tolist(), df["fco2_med_uatm"].tolist(), bins=3)

print(f"\n📈 Doğu Akdeniz SST Sen eğimi      : {mk_emed['theil_sen_slope']*10:+.3f} °C/10 yıl (p = {fmt_p(mk_emed['p_value'])})")
print(f"⚡ Pettitt kırılma yılı            : {pt_emed['changepoint_label']} (Δ = {pt_emed['regime_shift_delta']:+.3f} °C, p ≈ {fmt_p(pt_emed['p_value_approx'])})")
print(f"🔬 PCA açıklanan varyans          : PC1 = %{pca_info['explained_variance_ratio'][0]*100:.1f}, PC2 = %{pca_info['explained_variance_ratio'][1]*100:.1f}")
print(f"🧩 K-Means++ silhouette            : {km_info['silhouette_score']:.3f}")
print(f"🤖 Lag-1 Ridge test RMSE (SST)     : {ml_info['test_rmse']:.3f} °C (naif kalıcılık: {ml_info['naive_baseline_rmse']:.3f} °C)")
print(f"🔗 Transfer entropisi SST → fCO₂   : {te['transfer_entropy_x_to_y_bits']:.4f} bit")
print("ℹ️  N küçük (yıllık seri); kümeleme ve transfer entropisi keşif amaçlıdır, nedensellik kanıtı değildir.")

# ==============================================================================
# ADIM 4: 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)
colors = {"SST_EMED_BOX": "#f43f5e", "SST_AEGEAN_BOX": "#38bdf8", "SST_BLACKSEA_BOX": "#22c55e", "SST_WMED_BOX": "#f59e0b"}

ax1 = axes[0, 0]
for code, name in BOXES.items():
    yrs = sorted(annual_sst[code]); ax1.plot(yrs, [annual_sst[code][y] for y in yrs], color=colors[code], lw=1.8, label=name)
ax1.axvline(int(pt_emed["changepoint_label"]), color="#e2e8f0", ls=":", lw=1.4, label=f"Pettitt (D. Akdeniz) {pt_emed['changepoint_label']}")
ax1.set_title("A) Yıllık ortalama SST — NOAA OISST v2.1", color="#f8fafc", fontweight="bold")
ax1.set_ylabel("°C", color="#cbd5e1"); ax1.legend(**leg)

ax2 = axes[0, 1]
yrs = sorted(mhw_days["SST_EMED_BOX"])
ax2.bar([y - 0.2 for y in yrs], [mhw_days["SST_EMED_BOX"][y] for y in yrs], width=0.4, color="#f43f5e", label="Doğu Akdeniz")
ax2.bar([y + 0.2 for y in yrs], [mhw_days["SST_AEGEAN_BOX"].get(y, 0) for y in yrs], width=0.4, color="#38bdf8", label="Ege")
ax2.set_title("B) Deniz ısı dalgası günleri / yıl (Hobday 2016, eşik 1982–2011)", color="#f8fafc", fontweight="bold")
ax2.set_ylabel("gün", color="#cbd5e1"); ax2.legend(**leg)

ax3 = axes[1, 0]
oy = sorted(ohc_med); ax3.plot(oy, [ohc_med[y] / 1e9 for y in oy], color="#f59e0b", marker="o", lw=1.8, label="Akdeniz OHC anomalisi 0–700 m (GJ/m²)")
ax3.set_ylabel("GJ/m²", color="#f59e0b")
ax3t = ax3.twinx(); fy = sorted(fco2_med)
ax3t.plot(fy, [fco2_med[y] for y in fy], color="#14b8a6", marker="s", ls="--", lw=1.6, label="Akdeniz yüzey fCO₂ (µatm, SOCAT)")
ax3t.tick_params(colors="#5eead4"); ax3t.set_ylabel("µatm", color="#14b8a6")
ax3.set_title("C) Akdeniz ısı içeriği (CMEMS OMI) ve yüzey fCO₂ (SOCAT)", color="#f8fafc", fontweight="bold")
h1, l1 = ax3.get_legend_handles_labels(); h2, l2 = ax3t.get_legend_handles_labels(); ax3.legend(h1 + h2, l1 + l2, **leg)

ax4 = axes[1, 1]
pal = {1: "#22c55e", 2: "#f59e0b", 3: "#ef4444"}
for rid in sorted(set(df["regime_id"])):
    sub = df[df["regime_id"] == rid]
    ax4.scatter(sub["PC1"], sub["PC2"], s=60, c=pal.get(rid, "#94a3b8"), edgecolors="#ffffff", linewidth=0.5, label=f"Küme {rid}")
    for _, r in sub.iterrows():
        ax4.annotate(str(int(r["year"])), (r["PC1"], r["PC2"]), xytext=(3, 3), textcoords="offset points", color="#e2e8f0", fontsize=7)
ax4.set_title(f"D) PCA + K-Means++ (silhouette = {km_info['silhouette_score']:.2f})", color="#f8fafc", fontweight="bold")
ax4.set_xlabel(f"PC1 (%{pca_info['explained_variance_ratio'][0]*100:.0f})", color="#cbd5e1")
ax4.set_ylabel(f"PC2 (%{pca_info['explained_variance_ratio'][1]*100:.0f})", color="#cbd5e1"); ax4.legend(**leg)

fig.suptitle("GOMD — Vaka Çalışması 1: Doğu Akdeniz, Ege ve Karadeniz okyanus-iklim veri madenciliği (gerçek veri)",
             color="#ffffff", fontsize=12.5, fontweight="bold")
plt.tight_layout()
fig.savefig("cs01_multivariate_ocean_climate_data_mining.png", dpi=140, bbox_inches="tight", facecolor=fig.get_facecolor())
plt.show()