#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Küresel Karbon Bütçesi ve Okyanus: GCB 2025 Kaynak–Yutak Hesabı, Okyanus Yutağı Oranı, Mauna Loa CO₂, HOT/ALOHA Okyanus Asitlenmesi, NOAA Okyanus Isı İçeriği ve Arktik Deniz Buzu (1959–2025)
Global Carbon Budget & the Ocean: GCB 2025 Sources and Sinks, Ocean Sink Fraction, Mauna Loa CO₂, HOT/ALOHA Ocean Acidification, NOAA Ocean Heat Content & Arctic Sea Ice (1959–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):
  - gcb2025_fossil, _landuse, _atmosphere, _ocean_sink, _land_sink, _cement_sink, _imbalance — Global Carbon Budget 2025 (GCB_GLOBAL)
  - gml_mauna_loa_co2_ppm — NOAA GML / Scripps, Mauna Loa aylık CO₂
  - hot_seawater_ph — Hawaii Ocean Time-series, ALOHA istasyonu deniz suyu pH'ı
  - ncei_ohc_0_700m — NOAA NCEI okyanus ısı içeriği anomalisi (0–700 m)
  - nsidc_arctic_ice_min — NSIDC Arktik Eylül minimum deniz buzu alanı

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: Küresel Karbon Bütçesi 2025 — kaynaklar, yutaklar ve bütçe kapanışı (1959–2024)
# STEP 1: Global Carbon Budget 2025 — sources, sinks & budget closure (1959–2024)
# ==============================================================================
from gomd.stats import mann_kendall_sen, pettitt_changepoint, cointegration_eg_test

GCB = {k: fetch_gomd_annual_series(f"gcb2025_{k}", "GCB_GLOBAL")
       for k in ("fossil", "landuse", "atmosphere", "ocean_sink", "land_sink", "cement_sink", "imbalance")}
yrs = sorted(set.intersection(*[set(s) for s in GCB.values()]))
df = pd.DataFrame({"year": yrs, **{k: [GCB[k][y] for y in yrs] for k in GCB}})
df["emissions"] = df["fossil"] + df["landuse"]
closure = df["emissions"] - (df["atmosphere"] + df["ocean_sink"] + df["land_sink"] + df["cement_sink"] + df["imbalance"])
print(f"\n🧮 Bütçe kapanışı (kaynaklar − yutaklar − dengesizlik): en büyük mutlak fark {closure.abs().max():.3f} GtC/yıl")
df["airborne_frac"] = df["atmosphere"] / df["emissions"]
df["ocean_frac"] = df["ocean_sink"] / df["emissions"]
df["decade"] = (df["year"] // 10) * 10
dec = df[df["year"] >= 1960].groupby("decade")[["fossil", "landuse", "atmosphere", "ocean_sink", "land_sink", "airborne_frac", "ocean_frac"]].mean().round(3)
print("\n🌍 On yıllık ortalamalar, 1960+ (GtC/yıl; oranlar toplam emisyona göre; 2020 = 2020–2024):")
print(dec.to_string())
mk_oc = mann_kendall_sen(df["ocean_frac"].tolist(), t=yrs)
mk_af = mann_kendall_sen(df["airborne_frac"].tolist(), t=yrs)
print(f"🌊 Okyanus yutağı oranı trendi: {mk_oc['theil_sen_slope']*10:+.4f} /10 yıl (p = {fmt_p(mk_oc['p_value'])}); "
      f"havada kalan oran trendi: {mk_af['theil_sen_slope']*10:+.4f} /10 yıl (p = {fmt_p(mk_af['p_value'])})")
eg = cointegration_eg_test(df["emissions"].tolist(), df["ocean_sink"].tolist())
print(f"🔗 Engle-Granger (toplam emisyon → okyanus yutağı): β = {eg['long_run_slope_beta']:.3f}, ADF t = {eg['engle_granger_adf_t_stat']:.2f} "
      f"(%5 kritik {eg['critical_value_5pct']}; eşbütünleşik: {'evet' if eg['cointegrated_at_5pct'] else 'hayır'})")

# ==============================================================================
# ADIM 2: Mauna Loa CO₂ — yıllık ortalama, büyüme hızı ve GCB atmosfer artışıyla karşılaştırma
# STEP 2: Mauna Loa CO₂ — annual mean, growth rate & comparison with the GCB atmospheric growth
# ==============================================================================
mlo = fetch_gomd_records("gml_mauna_loa_co2_ppm", "WLD")
mlo["date"] = pd.to_datetime(mlo["period"].str[:7])
mlo_a = mlo.groupby(mlo["date"].dt.year)["value"].agg(["mean", "count"])
mlo_a = mlo_a[mlo_a["count"] == 12]["mean"]                   # yalnız 12 ayı tam yıllar
growth = mlo_a.diff().dropna()                                 # ppm / yıl
PPM_TO_GTC = 2.124                                             # 1 ppm CO₂ = 2,124 GtC (GCB'nin kullandığı dönüşüm)
cmp_y = [y for y in growth.index if y in GCB["atmosphere"]]
r = np.corrcoef([growth[y] * PPM_TO_GTC for y in cmp_y], [GCB["atmosphere"][y] for y in cmp_y])[0, 1]
print(f"\n🏔️ Mauna Loa yıllık ortalama CO₂: {mlo_a.iloc[0]:.2f} ppm ({mlo_a.index[0]}) → {mlo_a.iloc[-1]:.2f} ppm ({mlo_a.index[-1]})")
print(f"   Son 10 yılın ortalama artışı {growth.tail(10).mean():.2f} ppm/yıl; en hızlı artış {growth.max():.2f} ppm ({growth.idxmax()})")
print(f"   Mauna Loa artışı × {PPM_TO_GTC} ile GCB atmosfer artışının korelasyonu: r = {r:.3f} ({len(cmp_y)} yıl)")
print("   Not: GCB atmosfer artışı küresel deniz yüzeyi ağı ortalamasına dayanır; Mauna Loa tek istasyondur, küçük farklar beklenir.")

# ==============================================================================
# ADIM 3: Okyanus asitlenmesi — HOT/ALOHA istasyonu deniz suyu pH'ı (1988–2024)
# STEP 3: Ocean acidification — HOT/ALOHA station seawater pH (1988–2024)
# ==============================================================================
ph = fetch_gomd_records("hot_seawater_ph", "HOT_ALOHA")
ph["date"] = pd.to_datetime(ph["period"].str[:10])
ph_a = ph.groupby(ph["date"].dt.year)["value"].agg(["mean", "count"])
ph_a = ph_a[ph_a["count"] >= 4]["mean"]
mk_ph = mann_kendall_sen(ph_a.values.tolist(), t=ph_a.index.tolist())
both = [y for y in ph_a.index if y in mlo_a.index]
r_ph = np.corrcoef([ph_a[y] for y in both], [mlo_a[y] for y in both])[0, 1]
print(f"\n🧪 HOT/ALOHA pH: {len(ph_a)} yıl (yılda ≥4 ölçüm); Sen eğimi {mk_ph['theil_sen_slope']*10:+.4f} pH birimi /10 yıl (p = {fmt_p(mk_ph['p_value'])})")
print(f"   Yıllık pH ile Mauna Loa CO₂ korelasyonu: r = {r_ph:.3f}. Tek istasyondur; küresel okyanusu temsil etmez.")

# ==============================================================================
# ADIM 4: Okyanus ısı içeriği (NOAA NCEI, 0–700 m) ve Arktik Eylül deniz buzu (NSIDC)
# STEP 4: Ocean heat content (NOAA NCEI, 0–700 m) & Arctic September sea ice (NSIDC)
# ==============================================================================
ohc = fetch_gomd_annual_series("ncei_ohc_0_700m", "WLD")                  # 10²² J anomali
ice = fetch_gomd_annual_series("nsidc_arctic_ice_min", "ARCTIC")          # milyon km², Eylül minimumu
mk_ohc = mann_kendall_sen([ohc[y] for y in ohc if y >= 1979], t=[y for y in ohc if y >= 1979])
mk_ice = mann_kendall_sen(list(ice.values()), t=list(ice))
pt_ice = pettitt_changepoint(list(ice.values()), labels=list(ice))
print(f"\n🔥 OHC 0–700 m (1979+): Sen eğimi {mk_ohc['theil_sen_slope']*10:+.2f} ×10²² J /10 yıl (p = {fmt_p(mk_ohc['p_value'])})")
print(f"🧊 Arktik Eylül buz alanı: {ice[min(ice)]:.2f} → {ice[max(ice)]:.2f} milyon km² ({min(ice)}–{max(ice)}); "
      f"Sen eğimi {mk_ice['theil_sen_slope']*10:+.3f} milyon km² /10 yıl (p = {fmt_p(mk_ice['p_value'])}); Pettitt kırılma yılı {pt_ice['changepoint_label']}")

# ==============================================================================
# ADIM 5: 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]
ax1.plot(df["year"], df["fossil"], color="#f43f5e", lw=2, label="Fosil emisyon")
ax1.plot(df["year"], df["landuse"], color="#f59e0b", lw=1.5, label="Arazi kullanımı değişikliği")
ax1.plot(df["year"], df["atmosphere"], color="#e2e8f0", lw=1.3, label="Atmosferde birikim")
ax1.plot(df["year"], df["ocean_sink"], color="#38bdf8", lw=2, label="Okyanus yutağı")
ax1.plot(df["year"], df["land_sink"], color="#22c55e", lw=1.5, label="Kara yutağı")
ax1.set_ylabel("GtC / yıl", color="#cbd5e1")
ax1.set_title("A) Küresel Karbon Bütçesi 2025: kaynaklar ve yutaklar", color="#f8fafc", fontweight="bold"); ax1.legend(**leg)

ax2 = axes[0, 1]
ax2.plot(df["year"], df["airborne_frac"].rolling(5, center=True).mean(), color="#e2e8f0", lw=2, label="Havada kalan oran (5 yıllık ort.)")
ax2.plot(df["year"], df["ocean_frac"].rolling(5, center=True).mean(), color="#38bdf8", lw=2, label="Okyanusun aldığı oran (5 yıllık ort.)")
ax2.scatter(df["year"], df["ocean_frac"], color="#38bdf8", s=8, alpha=0.5)
ax2.set_ylabel("toplam emisyona oran", color="#cbd5e1")
ax2.set_title("B) Emisyonun nereye gittiği", color="#f8fafc", fontweight="bold"); ax2.legend(**leg)

ax3 = axes[1, 0]
ax3.plot(mlo["date"], mlo["value"], color="#f43f5e", lw=1.0, label="Mauna Loa CO₂ (aylık, ppm)")
ax3.set_ylabel("ppm", color="#f43f5e")
ax3t = ax3.twinx()
ax3t.plot(pd.to_datetime(ph_a.index.astype(str) + "-07-01"), ph_a.values, color="#38bdf8", marker="o", ms=3, lw=1.4, label="HOT/ALOHA yıllık pH")
ax3t.tick_params(colors="#7dd3fc"); ax3t.set_ylabel("pH", color="#38bdf8")
ax3.set_title("C) Atmosferik CO₂ ve okyanus asitlenmesi", color="#f8fafc", fontweight="bold")
h1, l1 = ax3.get_legend_handles_labels(); h2, l2 = ax3t.get_legend_handles_labels(); ax3.legend(h1 + h2, l1 + l2, loc="center left", **leg)

ax4 = axes[1, 1]
oy = [y for y in sorted(ohc) if y >= 1979]
ax4.plot(oy, [ohc[y] for y in oy], color="#f59e0b", marker="o", ms=3, lw=1.6, label="OHC 0–700 m (10²² J)")
ax4.set_ylabel("10²² J", color="#f59e0b")
ax4t = ax4.twinx(); iy = sorted(ice)
ax4t.plot(iy, [ice[y] for y in iy], color="#38bdf8", marker="s", ms=3, ls="--", lw=1.4, label="Arktik Eylül buz alanı (milyon km²)")
ax4t.tick_params(colors="#7dd3fc"); ax4t.set_ylabel("milyon km²", color="#38bdf8")
ax4.set_title("D) Okyanus ısı içeriği ve Arktik deniz buzu", color="#f8fafc", fontweight="bold")
h1, l1 = ax4.get_legend_handles_labels(); h2, l2 = ax4t.get_legend_handles_labels(); ax4.legend(h1 + h2, l1 + l2, loc="center left", **leg)

fig.suptitle("GOMD — Vaka Çalışması 7: Küresel karbon bütçesi, okyanus yutağı, asitlenme ve okyanus ısınması (gerçek veri)", color="#ffffff", fontsize=12.5, fontweight="bold")
plt.tight_layout()
fig.savefig("cs07_global_carbon_budget_mauna_loa_ocean_sink_cointegration.png", dpi=140, bbox_inches="tight", facecolor=fig.get_facecolor())
plt.show()