# -*- coding: utf-8 -*-
# ------------------------------------------------------------------------------
# Bir tropikal siklon bir limana kaç güne mal olur? 2.065 limanın 1980–2025 fırtına maruziyeti ve 1.321 liman-fırtına karşılaşmasında günlük uğrak kaybı
# How many days does a tropical cyclone cost a port? Storm exposure of 2,065 ports over 1980–2025 and the loss of daily port calls in 1,321 port–storm encounters
#
# GOMD Yayınlar · Makale 26 — https://gomd.pirireis.edu.tr/publications.php?id=rp26
#
# Bu dosya makalenin bütün modüllerini içerir. Replikasyon paketi klasöründe çalıştırıldığında modülleri yerlerine yazar ve sırayla, her birini ayrı süreçte çalıştırır;
# sonuçları `paper/out/` ve şekilleri `paper/fig/` klasörüne yazar.
# This file contains all modules of the paper. Run in the replication-package folder, it writes the modules to their places and runs them in order, each in its own process;
# it writes results to `paper/out/` and figures to `paper/fig/`.
#
# Kurulum / Setup:  pip install pandas numpy matplotlib pyarrow
# Çalıştırma / Run:  python rp26_tropical_cyclones_port_call_losses.py
#
# Storm tracks: International Best Track Archive for Climate Stewardship (IBTrACS) version 4, NOAA NCEI (Knapp et al. 2010), https://www.ncei.noaa.gov/products/international-best-track-archive . Port data: IMF PortWatch / UN Global Platform / University of Oxford (https://portwatch.imf.org), daily port activity derived from AIS.
# GOMD, Piri Reis Üniversitesi ile T.C. Çevre, Şehircilik ve İklim Değişikliği Bakanlığı Çevre Yönetimi Genel Müdürlüğü tarafından geliştirilen açık veri platformudur.
# GOMD is an open data platform developed by Piri Reis University together with the Directorate General of Environmental Management of the Ministry of Environment, Urbanisation and Climate Change of Türkiye.
# ------------------------------------------------------------------------------
import os, subprocess, sys

MODULES = {}
ORDER = ['tc/tc_prep.py', 'tc/tc_analysis.py', 'paper/fig26.py']
SKIP_IF_EXISTS = {'tc/tc_prep.py': 'tc/out/encounters.parquet'}   # hazırlık modülü: çıktısı pakette varsa atlanır / preparation module: skipped if its output is in the package
DIRS = ['tc', 'tc/out', 'paper', 'paper/out', 'paper/fig', 'data']

# ==============================================================================
# tc/tc_prep.py
# ==============================================================================
MODULES['tc/tc_prep.py'] = r'''
# P26 hazirlik: IBTrACS izleri (saatlik ara degerleme) x PortWatch limanlari -> liman-firtina karsilasmalari
import pandas as pd, numpy as np, json
use = ['SID', 'SEASON', 'BASIN', 'NAME', 'ISO_TIME', 'LAT', 'LON', 'WMO_WIND', 'USA_WIND', 'TRACK_TYPE', 'NATURE']
t = pd.read_csv('tc/ibtracs.since1980.csv', usecols=use, skiprows=[1], low_memory=False, na_values=[' ', ''], keep_default_na=False)
t = t[t.TRACK_TYPE.str.strip().str.lower() != 'spur'].copy()
t['time'] = pd.to_datetime(t.ISO_TIME); t['LAT'] = pd.to_numeric(t.LAT, errors='coerce'); t['LON'] = pd.to_numeric(t.LON, errors='coerce')
uw = pd.to_numeric(t.USA_WIND, errors='coerce'); ww = pd.to_numeric(t.WMO_WIND, errors='coerce')
t['wind'] = uw.where(uw.notna(), ww / 0.88)       # 1 dakikalik esdeger (kt); ABD kaydi yoksa WMO 10 dk / 0.88
t['wsrc'] = np.where(uw.notna(), 'USA', np.where(ww.notna(), 'WMO', 'none'))
t = t.dropna(subset=['LAT', 'LON']).sort_values(['SID', 'time'])
t['SEASON'] = pd.to_numeric(t.SEASON, errors='coerce').astype(int)
print('iz noktasi', len(t), 'firtina', t.SID.nunique(), 'sezon', t.SEASON.min(), t.SEASON.max(), 'ruzgar kaynagi', t.wsrc.value_counts().to_dict())
# saatlik ara degerleme
rows = []
for sid, g in t.groupby('SID', sort=False):
    g = g.drop_duplicates('time')
    if len(g) < 2: continue
    tt = g.time.values.astype('datetime64[h]').astype('int64'); lon = np.unwrap(np.radians(g.LON.values)); lat = g.LAT.values; w = g.wind.values
    hh = np.arange(tt[0], tt[-1] + 1)
    wi = np.interp(hh, tt[~np.isnan(w)], w[~np.isnan(w)]) if (~np.isnan(w)).sum() >= 2 else np.full(len(hh), np.nan)
    if (~np.isnan(w)).sum() >= 2:
        first, last = tt[~np.isnan(w)][0], tt[~np.isnan(w)][-1]; wi[(hh < first) | (hh > last)] = np.nan
    lo = np.degrees(np.interp(hh, tt, lon)); lo = (lo + 180) % 360 - 180
    rows.append(pd.DataFrame({'SID': sid, 'h': hh, 'lat': np.interp(hh, tt, lat), 'lon': lo, 'wind': wi, 'season': g.SEASON.iloc[0], 'basin': g.BASIN.iloc[0], 'name': g.NAME.iloc[0]}))
H = pd.concat(rows, ignore_index=True); H = H[H.wind >= 34].reset_index(drop=True)
print('saatlik nokta (>=34 kt)', len(H), 'firtina', H.SID.nunique())
P = pd.read_csv('data/pw/ports_global.csv')[['portid', 'portname', 'country', 'ISO3', 'continent', 'lat', 'lon']]
pl = np.radians(P.lat.values); pn = np.radians(P.lon.values)
out = []
R = 6371.0
for i0 in range(0, len(H), 20000):
    h = H.iloc[i0:i0 + 20000]; la = np.radians(h.lat.values)[:, None]; lo = np.radians(h.lon.values)[:, None]
    box = (np.abs(h.lat.values[:, None] - P.lat.values[None, :]) < 2.0)
    ii, jj = np.where(box)
    if len(ii) == 0: continue
    a = np.sin((pl[jj] - la[ii, 0]) / 2) ** 2 + np.cos(la[ii, 0]) * np.cos(pl[jj]) * np.sin((pn[jj] - lo[ii, 0]) / 2) ** 2
    d = 2 * R * np.arcsin(np.sqrt(a)); k = d <= 200
    out.append(pd.DataFrame({'SID': h.SID.values[ii[k]], 'h': h.h.values[ii[k]], 'portid': P.portid.values[jj[k]], 'dist': d[k], 'wind': h.wind.values[ii[k]]}))
E_ = pd.concat(out, ignore_index=True)
print('saatlik liman-firtina eslesmesi (<=200 km)', len(E_))
def enc(radius):
    e = E_[E_.dist <= radius]
    g = e.groupby(['SID', 'portid'])
    r = g.agg(first_h=('h', 'min'), last_h=('h', 'max'), hours=('h', 'count'), min_dist=('dist', 'min'), max_wind=('wind', 'max')).reset_index()
    r['radius'] = radius; return r
ENC = pd.concat([enc(r) for r in (50, 100, 200)], ignore_index=True)
meta = H.groupby('SID').agg(season=('season', 'first'), basin=('basin', 'first'), name=('name', 'first'), life_max=('wind', 'max')).reset_index()
ENC = ENC.merge(meta, on='SID'); ENC['t0'] = pd.to_datetime(ENC.first_h.astype('int64') * 3600, unit='s').dt.floor('D')
ENC['year'] = ENC.t0.dt.year
ENC.to_parquet('tc/out/encounters.parquet'); P.to_csv('tc/out/ports.csv', index=False)
print(ENC.groupby('radius').agg(n=('SID', 'size'), storms=('SID', 'nunique'), ports=('portid', 'nunique')))
json.dump({'track_points': int(len(t)), 'storms_all': int(t.SID.nunique()), 'storms_ts': int(H.SID.nunique()), 'season_min': int(t.SEASON.min()), 'season_max': int(t.SEASON.max()),
           'wind_src': {k: int(v) for k, v in t.wsrc.value_counts().items()}, 'ports': int(len(P))}, open('tc/out/prep.json', 'w'))
'''

# ==============================================================================
# tc/tc_analysis.py
# ==============================================================================
MODULES['tc/tc_analysis.py'] = r'''
# P26 analiz: liman tropikal siklon maruziyeti (1980-2025) ve gunluk ugrak olay calismasi (2019-2026)
import pandas as pd, numpy as np, json
rng = np.random.default_rng(26)
ENC = pd.read_parquet('tc/out/encounters.parquet'); P = pd.read_csv('tc/out/ports.csv')
D = pd.read_parquet('data/pwd/daily.parquet', columns=['date', 'portid', 'portcalls', 'import', 'export'])
D['cargo'] = D['import'] + D['export']
calls = D.pivot(index='date', columns='portid', values='portcalls').sort_index(); cargo = D.pivot(index='date', columns='portid', values='cargo').sort_index()
full = pd.date_range(calls.index.min(), calls.index.max()); calls = calls.reindex(full); cargo = cargo.reindex(full)
yr = D[(D.date >= '2019-01-01') & (D.date <= '2025-12-31')]
w = (yr.groupby('portid').cargo.sum() / 7.0).rename('cargo_yr'); cy = (yr.groupby('portid').portcalls.sum() / 7.0).rename('calls_yr')
P = P.merge(w, on='portid', how='left').merge(cy, on='portid', how='left').fillna({'cargo_yr': 0, 'calls_yr': 0})
res = {}
# ---------------- 1) maruziyet 1980-2025
def cls(w): return np.where(w >= 96, 'major', np.where(w >= 64, 'hurricane', 'storm'))
E1 = ENC[(ENC.radius == 100) & (ENC.year >= 1980) & (ENC.year <= 2025)].copy(); E1['cls'] = cls(E1.max_wind.values)
res['n_encounters_100'] = int(len(E1)); res['n_storms_100'] = int(E1.SID.nunique()); res['n_ports_hit_100'] = int(E1.portid.nunique()); res['n_ports'] = int(len(P))
ex = E1.groupby('portid').agg(ts_events=('SID', 'size'), hu_events=('max_wind', lambda x: int((x >= 64).sum())), major_events=('max_wind', lambda x: int((x >= 96).sum())), max_wind=('max_wind', 'max')).reset_index()
for a, b, nm in ((1980, 2002, 'early'), (2003, 2025, 'late')):
    s = E1[(E1.year >= a) & (E1.year <= b)]
    ex = ex.merge(s.groupby('portid').size().rename('ts_' + nm).reset_index(), on='portid', how='left').merge(s[s.max_wind >= 64].groupby('portid').size().rename('hu_' + nm).reset_index(), on='portid', how='left')
ex = P.merge(ex, on='portid', how='left').fillna({c: 0 for c in ['ts_events', 'hu_events', 'major_events', 'ts_early', 'ts_late', 'hu_early', 'hu_late']})
ex['ts_per_decade'] = ex.ts_events / 4.6; ex['hu_per_decade'] = ex.hu_events / 4.6
ex.to_csv('tc/out/exposure_port.csv', index=False)
tot_c = ex.cargo_yr.sum()
res['share_ports_exposed'] = float((ex.ts_events > 0).mean()); res['share_cargo_exposed'] = float(ex.cargo_yr[ex.ts_events > 0].sum() / tot_c)
res['share_ports_hu'] = float((ex.hu_events > 0).mean()); res['share_cargo_hu'] = float(ex.cargo_yr[ex.hu_events > 0].sum() / tot_c)
res['share_ports_ts_decadal'] = float((ex.ts_per_decade >= 1).mean()); res['share_cargo_ts_decadal'] = float(ex.cargo_yr[ex.ts_per_decade >= 1].sum() / tot_c)
res['mean_ts_per_decade_port'] = float(ex.ts_per_decade.mean()); res['tw_ts_per_decade'] = float((ex.ts_per_decade * ex.cargo_yr).sum() / tot_c)
res['mean_hu_per_decade_port'] = float(ex.hu_per_decade.mean()); res['tw_hu_per_decade'] = float((ex.hu_per_decade * ex.cargo_yr).sum() / tot_c)
top = ex.sort_values('ts_events', ascending=False).head(15)[['portname', 'country', 'ts_events', 'hu_events', 'major_events', 'cargo_yr']]
res['top_ports'] = top.to_dict('records')
ctry = ex.groupby('country').apply(lambda g: pd.Series({'ports': len(g), 'ts_pd': g.ts_per_decade.mean(), 'tw_ts_pd': (g.ts_per_decade * g.cargo_yr).sum() / max(g.cargo_yr.sum(), 1), 'cargo': g.cargo_yr.sum()})).reset_index()
ctry.to_csv('tc/out/exposure_country.csv', index=False)
res['top_countries_tw'] = ctry[ctry.ports >= 5].sort_values('tw_ts_pd', ascending=False).head(12).to_dict('records')
# yillik seriler: karsilasma sayisi, ticaret agirlikli, enlem
E1 = E1.merge(P[['portid', 'lat', 'cargo_yr']], on='portid')
ann = E1.groupby('year').apply(lambda g: pd.Series({'enc': len(g), 'enc_hu': int((g.max_wind >= 64).sum()), 'enc_major': int((g.max_wind >= 96).sum()), 'tw': g.cargo_yr.sum() / tot_c, 'tw_hu': g.cargo_yr[g.max_wind >= 64].sum() / tot_c,
                                               'abs_lat': g.lat.abs().mean(), 'storms': g.SID.nunique(), 'mean_wind': g.max_wind.mean(), 'share_hu': (g.max_wind >= 64).mean()})).reset_index()
ann.to_csv('tc/out/annual.csv', index=False)
def trend(y, x):
    x = np.asarray(x, float); y = np.asarray(y, float); n = len(x); xm = x - x.mean(); b = (xm * (y - y.mean())).sum() / (xm ** 2).sum(); e = y - y.mean() - b * xm
    r1 = np.corrcoef(e[:-1], e[1:])[0, 1]; neff = n * (1 - r1) / (1 + r1) if r1 > 0 else n
    se = np.sqrt((e ** 2).sum() / max(neff - 2, 1) / (xm ** 2).sum()); return {'slope': float(b), 'se': float(se), 't': float(b / se), 'mean': float(y.mean()), 'pct_per_decade': float(1000 * b / y.mean())}
for c in ('enc', 'enc_hu', 'enc_major', 'tw', 'tw_hu', 'abs_lat', 'mean_wind', 'share_hu'): res['trend_' + c] = trend(ann[c], ann.year)
res['early_late'] = {k: [float(ann[ann.year <= 2002][k].mean()), float(ann[ann.year >= 2003][k].mean())] for k in ('enc', 'enc_hu', 'enc_major', 'tw', 'tw_hu', 'abs_lat', 'share_hu')}
# havza bazinda
bas = E1.groupby(['basin', 'year']).size().unstack(0).reindex(range(1980, 2026)).fillna(0)
res['basin_trend'] = {b: trend(bas[b].values, bas.index.values) | {'total': int(bas[b].sum())} for b in bas.columns}
for r in (50, 200):
    e = ENC[(ENC.radius == r) & (ENC.year >= 1980) & (ENC.year <= 2025)]; a = e.groupby('year').size().reindex(range(1980, 2026)).fillna(0)
    res['trend_enc_r%d' % r] = trend(a.values, a.index.values) | {'n': int(len(e))}
# ---------------- 2) olay calismasi 2019-2026
W0, W1 = -7, 28
def events(radius=100, minwind=34):
    e = ENC[(ENC.radius == radius) & (ENC.max_wind >= minwind) & (ENC.t0 >= '2019-03-01') & (ENC.t0 <= calls.index.max() - pd.Timedelta(days=W1))].sort_values(['portid', 't0']).copy()
    keep = []; last = {}
    for r in e.itertuples():
        if r.portid in last and (r.t0 - last[r.portid]).days < 36: continue
        last[r.portid] = r.t0; keep.append(r.Index)
    return e.loc[keep]
def window(e, mat, base_min=0.5, placebo_shift=0):
    X = []; meta = []
    for r in e.itertuples():
        if r.portid not in mat.columns: continue
        t0 = r.t0 + pd.Timedelta(days=placebo_shift); s = mat[r.portid]
        if t0 + pd.Timedelta(days=W1) > s.index[-1] or t0 - pd.Timedelta(days=56) < s.index[0]: continue
        b = s.loc[t0 - pd.Timedelta(days=56): t0 - pd.Timedelta(days=8)]
        if b.isna().any() or b.mean() < base_min: continue
        x = s.loc[t0 + pd.Timedelta(days=W0): t0 + pd.Timedelta(days=W1)].values
        if np.isnan(x).any(): continue
        X.append(x / b.mean()); meta.append((r.SID, r.portid, r.max_wind, r.min_dist, b.mean(), r.basin, t0))
    return np.array(X), pd.DataFrame(meta, columns=['SID', 'portid', 'max_wind', 'min_dist', 'base', 'basin', 't0'])
days = np.arange(W0, W1 + 1)
def summarize(X, M, nb=1000):
    wts = M.base.values; pooled = (X * wts[:, None]).sum(0) / wts.sum()
    loss_i = (1 - X[:, (days >= -2) & (days <= 14)]).sum(1); reb_i = (X[:, (days >= 15) & (days <= 28)] - 1).sum(1)
    sid = M.SID.values; u = np.unique(sid); idx = {s: np.where(sid == s)[0] for s in u}
    bl = []; bp = []; br = []; bm = []
    for _ in range(nb):
        pick = np.concatenate([idx[s] for s in rng.choice(u, len(u))]); ww = wts[pick]
        pp = (X[pick] * ww[:, None]).sum(0) / ww.sum(); bp.append(pp)
        bl.append((1 - pp[(days >= -2) & (days <= 14)]).sum()); br.append((pp[(days >= 15) & (days <= 28)] - 1).sum()); bm.append(pp.min())
    bp = np.array(bp)
    return {'n': int(len(X)), 'storms': int(len(u)), 'ports': int(M.portid.nunique()), 'pooled': pooled.tolist(), 'lo': np.percentile(bp, 2.5, 0).tolist(), 'hi': np.percentile(bp, 97.5, 0).tolist(),
            'min_ratio': float(pooled.min()), 'min_day': int(days[pooled.argmin()]), 'min_ci': [float(np.percentile(bm, 2.5)), float(np.percentile(bm, 97.5))],
            'loss_days': float((1 - pooled[(days >= -2) & (days <= 14)]).sum()), 'loss_ci': [float(np.percentile(bl, 2.5)), float(np.percentile(bl, 97.5))],
            'rebound_days': float((pooled[(days >= 15) & (days <= 28)] - 1).sum()), 'rebound_ci': [float(np.percentile(br, 2.5)), float(np.percentile(br, 97.5))],
            'loss_median_event': float(np.median(loss_i)), 'loss_mean_event': float(loss_i.mean()), 'share_events_loss_gt1': float((loss_i > 1).mean()),
            'recover_day': int(next((d for d in days[days >= 0] if pooled[days == d][0] >= 0.95 and pooled[(days >= d) & (days <= d + 2)].mean() >= 0.95), 99))}
ev = events(100, 34); X, M = window(ev, calls)
res['es_days'] = days.tolist(); res['es_all'] = summarize(X, M)
M['cls'] = cls(M.max_wind.values)
for c in ('storm', 'hurricane', 'major'):
    k = (M.cls == c).values; res['es_' + c] = summarize(X[k], M[k])
k = (M.max_wind >= 64).values; res['es_hu_all'] = summarize(X[k], M[k])
# plasebo: ayni liman, ayni takvim gunu, olay olmayan yillar
def placebo():
    out = []; evset = {(r.portid, r.t0.year) for r in ENC[ENC.radius == 200].itertuples()}
    for sh in (-365 * 2, -365, 365, 365 * 2):
        e2 = ev.copy(); Xp, Mp = window(e2, calls, placebo_shift=sh)
        if len(Xp) == 0: continue
        t = Mp.t0
        ok = np.array([not any((p, y) in evset for y in (tt.year,)) or True for p, tt in zip(Mp.portid, t)])
        # olay penceresine baska firtina dusuyorsa at
        near = ENC[ENC.radius == 200][['portid', 't0']].rename(columns={'t0': 'te'})
        mm = Mp.reset_index().merge(near, on='portid', how='left'); mm['gap'] = (mm.te - mm.t0).dt.days.abs()
        bad = set(mm[mm.gap <= 35]['index']); keep = np.array([i not in bad for i in range(len(Mp))])
        out.append((Xp[keep], Mp[keep]))
    Xp = np.vstack([o[0] for o in out]); Mp = pd.concat([o[1] for o in out], ignore_index=True); return Xp, Mp
Xp, Mp = placebo(); res['es_placebo'] = summarize(Xp, Mp)
# heterojenlik (kasirga ve ustu)
k64 = (M.max_wind >= 64).values
het = {}
het['dist_le50'] = summarize(X[k64 & (M.min_dist <= 50).values], M[k64 & (M.min_dist <= 50).values]); het['dist_50_100'] = summarize(X[k64 & (M.min_dist > 50).values], M[k64 & (M.min_dist > 50).values])
q = M.base.quantile([1 / 3, 2 / 3]).values
for nm, kk in (('small', M.base <= q[0]), ('mid', (M.base > q[0]) & (M.base <= q[1])), ('large', M.base > q[1])): het['size_' + nm] = summarize(X[k64 & kk.values], M[k64 & kk.values])
for b in M[k64].basin.value_counts().index[:4]: het['basin_' + b] = summarize(X[k64 & (M.basin == b).values], M[k64 & (M.basin == b).values])
res['het'] = het; res['size_terciles'] = [float(q[0]), float(q[1])]
# yuk (tonaj)
Xc, Mc = window(ev, cargo, base_min=1.0); kc = (Mc.max_wind >= 64).values
res['es_cargo_all'] = summarize(Xc, Mc); res['es_cargo_hu'] = summarize(Xc[kc], Mc[kc])
# yaricap duyarliligi
for r in (50, 200):
    e2 = events(r, 64); X2, M2 = window(e2, calls); res['es_hu_r%d' % r] = summarize(X2, M2)
# doz-tepki: olay bazinda kayip ~ ruzgar
loss_i = (1 - X[:, (days >= -2) & (days <= 14)]).sum(1); M['loss'] = loss_i; M.to_csv('tc/out/events.csv', index=False)
xw = M.max_wind.values; A = np.c_[np.ones(len(xw)), (xw - 34) / 10, M.min_dist.values / 50]; wt = M.base.values
beta = np.linalg.lstsq(A * np.sqrt(wt)[:, None], loss_i * np.sqrt(wt), rcond=None)[0]
bb = []
u = M.SID.unique(); idx = {s: np.where(M.SID.values == s)[0] for s in u}
for _ in range(1000):
    pick = np.concatenate([idx[s] for s in rng.choice(u, len(u))]); bb.append(np.linalg.lstsq(A[pick] * np.sqrt(wt[pick])[:, None], loss_i[pick] * np.sqrt(wt[pick]), rcond=None)[0])
bb = np.array(bb); res['dose'] = {'const': float(beta[0]), 'per10kt': float(beta[1]), 'per50km': float(beta[2]), 'per10kt_ci': [float(np.percentile(bb[:, 1], 2.5)), float(np.percentile(bb[:, 1], 97.5))], 'per50km_ci': [float(np.percentile(bb[:, 2], 2.5)), float(np.percentile(bb[:, 2], 97.5))]}
# ---------------- 3) beklenen yillik kayip: maruziyet x olay basina kayip
Lts, Lhu, Lmaj = res['es_storm']['loss_days'], res['es_hurricane']['loss_days'], res['es_major']['loss_days']
ex['exp_loss_days_yr'] = ((ex.ts_events - ex.hu_events) * Lts + (ex.hu_events - ex.major_events) * Lhu + ex.major_events * Lmaj) / 46.0
ex.to_csv('tc/out/exposure_port.csv', index=False)
res['exp_loss'] = {'mean_port': float(ex.exp_loss_days_yr.mean()), 'tw': float((ex.exp_loss_days_yr * ex.cargo_yr).sum() / tot_c), 'p95_port': float(ex.exp_loss_days_yr.quantile(.95)), 'max_port': float(ex.exp_loss_days_yr.max()),
                   'ports_gt1': int((ex.exp_loss_days_yr >= 1).sum()), 'cargo_share_gt1': float(ex.cargo_yr[ex.exp_loss_days_yr >= 1].sum() / tot_c),
                   'top': ex.sort_values('exp_loss_days_yr', ascending=False).head(12)[['portname', 'country', 'exp_loss_days_yr', 'ts_events', 'hu_events', 'major_events']].to_dict('records')}
res['events_2019_2026'] = {'n': int(len(M)), 'storms': int(M.SID.nunique()), 'ports': int(M.portid.nunique()), 'by_cls': M.cls.value_counts().to_dict(), 'first': str(M.t0.min().date()), 'last': str(M.t0.max().date())}
json.dump(res, open('tc/out/results.json', 'w'), indent=1, default=float, allow_nan=False)
for k in ('n_encounters_100', 'n_storms_100', 'n_ports_hit_100', 'share_ports_exposed', 'share_cargo_exposed', 'share_ports_hu', 'share_cargo_hu', 'mean_ts_per_decade_port', 'tw_ts_per_decade', 'early_late', 'trend_enc', 'trend_enc_hu', 'trend_tw', 'trend_abs_lat', 'trend_share_hu', 'events_2019_2026', 'dose', 'exp_loss'): print(k, json.dumps(res[k], default=float)[:420])
for k in ('es_all', 'es_storm', 'es_hurricane', 'es_major', 'es_hu_all', 'es_placebo', 'es_cargo_hu', 'es_hu_r50', 'es_hu_r200'): print(k, {a: (round(b, 3) if isinstance(b, float) else b) for a, b in res[k].items() if a not in ('pooled', 'lo', 'hi')})
for k, v in res['het'].items(): print('het', k, v['n'], round(v['loss_days'], 2), [round(x, 2) for x in v['loss_ci']], round(v['min_ratio'], 2))
'''

# ==============================================================================
# paper/fig26.py
# ==============================================================================
MODULES['paper/fig26.py'] = r'''
import sys; sys.path.insert(0, '.')
from figstyle import *
import json, numpy as np, pandas as pd
R = json.load(open('tc/out/results.json')); P = pd.read_csv('tc/out/exposure_port.csv'); A = pd.read_csv('tc/out/annual.csv'); EV = pd.read_csv('tc/out/events.csv')
# ek kuyruk istatistikleri
tail = {'n': int(len(EV)), 'share_gt3': float((EV.loss > 3).mean()), 'share_gt5': float((EV.loss > 5).mean()), 'share_le0': float((EV.loss <= 0).mean()), 'p90': float(EV.loss.quantile(.9)), 'p99': float(EV.loss.quantile(.99)),
        'major_share_gt5': float((EV[EV.cls == 'major'].loss > 5).mean()), 'top10_share_of_total_loss': float(EV.loss.clip(lower=0).nlargest(int(len(EV) * .1)).sum() / EV.loss.clip(lower=0).sum())}
json.dump(tail, open('tc/out/tail.json', 'w'), indent=1); print(tail)
# Sekil 1: harita
fig, ax = plt.subplots(figsize=(7.4, 3.7))
z = P[P.ts_events == 0]; ax.scatter(z.lon, z.lat, s=2, c='#c8ccd0', lw=0)
h = P[P.ts_events > 0].sort_values('ts_per_decade'); sc = ax.scatter(h.lon, h.lat, s=4 + 2.2 * h.ts_per_decade, c=h.ts_per_decade, cmap='YlOrRd', vmin=0, vmax=10, lw=.2, edgecolors='#444')
cb = fig.colorbar(sc, ax=ax, shrink=.75, pad=.01); cb.set_label('On yılda fırtına karşılaşması (≥34 kt, ≤100 km)')
ax.set_xlim(-180, 180); ax.set_ylim(-60, 75); ax.set_xlabel('Boylam'); ax.set_ylabel('Enlem'); ax.grid(False)
fig.savefig('paper/fig/Figure_1.png'); plt.close()
# Sekil 2: yillik seriler
fig, ax = plt.subplots(1, 2, figsize=(7.4, 3))
ax[0].plot(A.year, A.enc, color=C['b'], label='Tümü (≥34 kt)'); ax[0].plot(A.year, A.enc_hu, color=C['o'], label='Kasırga gücü (≥64 kt)'); ax[0].plot(A.year, A.enc_major, color=C['r'], label='Majör (≥96 kt)')
for col, c in (('enc', 'b'), ('enc_hu', 'o'), ('enc_major', 'r')):
    b = np.polyfit(A.year, A[col], 1); ax[0].plot(A.year, np.polyval(b, A.year), color=C[c], ls='--', lw=.8)
ax[0].set_title('(a) Liman–fırtına karşılaşmaları'); ax[0].set_ylabel('Yıllık sayı'); ax[0].legend(loc='upper left', ncol=1)
ax[1].plot(A.year, 10 * A.tw, color=C['b'], label='Tümü'); ax[1].plot(A.year, 10 * A.tw_hu, color=C['o'], label='Kasırga gücü')
ax[1].set_title('(b) Yük ağırlıklı maruziyet'); ax[1].set_ylabel('Liman başına karşılaşma (on yıllık eşdeğer)'); ax[1].legend()
fig.tight_layout(); fig.savefig('paper/fig/Figure_2.png'); plt.close()
# Sekil 3: olay calismasi
d = np.array(R['es_days']); fig, ax = plt.subplots(1, 2, figsize=(7.4, 3.1))
for k, c, lab in (('es_storm', 'b', 'Tropikal fırtına (34–63 kt)'), ('es_hurricane', 'o', 'Kasırga (64–95 kt)'), ('es_major', 'r', 'Majör kasırga (≥96 kt)')):
    e = R[k]; m = d <= 14; ax[0].plot(d[m], 100 * np.array(e['pooled'])[m], color=C[c], label='%s, n = %d' % (lab, e['n']), marker='o', ms=2.5); ax[0].fill_between(d[m], 100 * np.array(e['lo'])[m], 100 * np.array(e['hi'])[m], color=C[c], alpha=.13, lw=0)
ax[0].axhline(100, color='k', lw=.6); ax[0].axvline(0, color='k', lw=.6, ls=':'); ax[0].set_xlabel('En yakın geçişe göre gün'); ax[0].set_ylabel('Uğrak sayısı (olağan düzey = 100)'); ax[0].set_title('(a) Fırtına şiddetine göre'); ax[0].legend(loc='lower right', fontsize=7)
for k, c, lab in (('es_all', 'b', 'Gerçek olaylar'), ('es_placebo', 'gr', 'Plasebo tarihler')):
    e = R[k]; ax[1].plot(d, 100 * np.array(e['pooled']), color=C[c], label='%s, n = %d' % (lab, e['n'])); ax[1].fill_between(d, 100 * np.array(e['lo']), 100 * np.array(e['hi']), color=C[c], alpha=.15, lw=0)
ax[1].axhline(100, color='k', lw=.6); ax[1].axvline(0, color='k', lw=.6, ls=':'); ax[1].set_xlabel('En yakın geçişe göre gün'); ax[1].set_title('(b) Tüm olaylar ve plasebo'); ax[1].legend(loc='lower right')
fig.tight_layout(); fig.savefig('paper/fig/Figure_3.png'); plt.close()
# Sekil 4: doz-tepki ve beklenen yillik kayip
fig, ax = plt.subplots(1, 2, figsize=(7.4, 3.2))
bins = [34, 50, 64, 80, 96, 113, 200]; lab = ['34–49', '50–63', '64–79', '80–95', '96–112', '≥113']; mm = []; lo = []; hi = []; nn = []; rng = np.random.default_rng(1)
for a, b in zip(bins[:-1], bins[1:]):
    q = EV[(EV.max_wind >= a) & (EV.max_wind < b)]; st = q.SID.unique(); g = {s: x.loss.values for s, x in q.groupby('SID')}
    bs = [np.concatenate([g[s] for s in rng.choice(st, len(st))]).mean() for _ in range(1000)]; mm.append(q.loss.mean()); lo.append(np.percentile(bs, 2.5)); hi.append(np.percentile(bs, 97.5)); nn.append(len(q))
x = np.arange(len(lab)); ax[0].bar(x, mm, color=C['lb'], edgecolor=C['b'], lw=.6); ax[0].errorbar(x, mm, yerr=[np.array(mm) - lo, np.array(hi) - mm], fmt='none', ecolor=C['k'], lw=.8, capsize=2)
for i, n_ in enumerate(nn): ax[0].text(i, -0.42, 'n=%d' % n_, ha='center', fontsize=6.5, color='#555')
ax[0].set_xticks(x); ax[0].set_xticklabels(lab, fontsize=7); ax[0].set_xlabel('Limana 100 km içindeki en yüksek rüzgâr (kt)'); ax[0].set_ylabel('Kayıp (liman-günü)'); ax[0].set_title('(a) Şiddete göre uğrak kaybı'); ax[0].axhline(0, color='k', lw=.6); ax[0].set_ylim(-0.6, None)
json.dump({'bins': lab, 'mean': mm, 'lo': lo, 'hi': hi, 'n': nn}, open('tc/out/dose_bins.json', 'w'))
t = P.sort_values('exp_loss_days_yr').tail(15); ax[1].barh(t.portname + ' (' + t.ISO3 + ')', t.exp_loss_days_yr, color=C['o']); ax[1].set_xlabel('Beklenen yıllık kayıp (liman-günü)'); ax[1].set_title('(b) En yüksek beklenen kayıp'); ax[1].tick_params(axis='y', labelsize=7)
fig.tight_layout(); fig.savefig('paper/fig/Figure_4.png'); plt.close()
print(json.load(open('tc/out/dose_bins.json')))
# makalenin okudugu birlesik sonuc dosyasi
import os, shutil
os.makedirs('paper/out', exist_ok=True)
R['tail'] = tail; R['dose_bins'] = json.load(open('tc/out/dose_bins.json')); R['prep'] = json.load(open('tc/out/prep.json'))
json.dump(R, open('paper/out/results.json', 'w'))
for f in ('exposure_port.csv', 'exposure_country.csv', 'annual.csv', 'events.csv'): shutil.copy('tc/out/' + f, 'paper/out/' + f)
'''

# ==============================================================================
# figstyle.py
# ==============================================================================
MODULES['figstyle.py'] = r'''
import matplotlib; matplotlib.use('Agg')
import matplotlib.pyplot as plt
plt.rcParams.update({'font.family': 'DejaVu Sans', 'font.size': 9, 'axes.spines.top': False, 'axes.spines.right': False, 'axes.grid': True, 'grid.alpha': .25, 'grid.linewidth': .5,
                     'axes.titlesize': 9.5, 'axes.titleweight': 'bold', 'axes.titlelocation': 'left', 'legend.frameon': False, 'legend.fontsize': 8, 'savefig.dpi': 300, 'savefig.bbox': 'tight', 'figure.dpi': 110})
C = {'b': '#1f5a94', 'r': '#c0392b', 'g': '#2e7d5b', 'o': '#d98200', 'k': '#333333', 'p': '#6a4c93', 'lb': '#8fb8de', 'gr': '#9aa0a6'}
'''

# ==============================================================================
# dl_all.py
# ==============================================================================
MODULES['dl_all.py'] = r'''
# n5: bes makale icin ham veri indirme. Kullanim: dl_all.py tc|ice|mp|ices|obis
import sys, os, json, time, urllib.request, urllib.parse, concurrent.futures as cf
UA = {'User-Agent': 'gomd-research/1.0 (gomd.pirireis.edu.tr)'}
def get(url, timeout=120, tries=4):
    last = None
    for k in range(tries):
        try:
            with urllib.request.urlopen(urllib.request.Request(url, headers=UA), timeout=timeout) as r: return r.read()
        except Exception as e:
            last = e; time.sleep(2 + 3 * k)
    raise last
what = sys.argv[1]
if what == 'tc':
    u = 'https://www.ncei.noaa.gov/data/international-best-track-archive-for-climate-stewardship-ibtracs/v04r01/access/csv/ibtracs.since1980.list.v04r01.csv'
    b = get(u, 900); open('tc/ibtracs.since1980.csv.part', 'wb').write(b); os.replace('tc/ibtracs.since1980.csv.part', 'tc/ibtracs.since1980.csv'); print('tc', len(b))
elif what == 'ice':
    for h in 'NS':
        f = '%s_Sea_Ice_Index_Regional_Daily_Data_G02135_v4.0.xlsx' % h
        b = get('https://noaadata.apps.nsidc.org/NOAA/G02135/seaice_analysis/' + f, 600); open('ice/' + f, 'wb').write(b); print(f, len(b))
    for h, d in (('N', 'north'), ('S', 'south')):
        f = '%s_seaice_extent_daily_v4.0.csv' % h
        b = get('https://noaadata.apps.nsidc.org/NOAA/G02135/%s/daily/data/%s' % (d, f), 600); open('ice/' + f, 'wb').write(b); print(f, len(b))
elif what == 'mp':
    base = 'https://services2.arcgis.com/C8EMgrsFcRFL6LrL/arcgis/rest/services/Marine_Microplastics_WGS84/FeatureServer/0/query'
    out = open('mp/ncei.jsonl.part', 'w', encoding='utf-8'); n = 0; off = 0
    while True:
        q = urllib.parse.urlencode(dict(where='1=1', outFields='*', returnGeometry='false', orderByFields='OBJECTID', resultOffset=off, resultRecordCount=2000, f='json'))
        d = json.loads(get(base + '?' + q)); fs = d.get('features', [])
        for x in fs: out.write(json.dumps(x['attributes'], ensure_ascii=False) + '\n')
        n += len(fs); off += len(fs)
        if len(fs) < 2000 and not d.get('exceededTransferLimit'): break
        if not fs: break
    out.close(); os.replace('mp/ncei.jsonl.part', 'mp/ncei.jsonl'); print('mp', n)
elif what == 'ices':
    A = 'https://sag.ices.dk/SAG_API/api/'
    lists = {}
    for y in range(2014, 2027):
        try: lists[y] = json.loads(get(A + 'StockList?year=%d' % y))
        except Exception as e: print('liste', y, e); lists[y] = []
        print(y, len(lists[y]), flush=True)
    json.dump(lists, open('ices/stocklists.json', 'w'))
    keys = sorted({s['AssessmentKey'] for y in lists for s in lists[y]})
    def one(k):
        f = 'ices/a/%d.json' % k
        if os.path.isfile(f): return 0
        try:
            s = json.loads(get(A + 'SummaryTable?assessmentKey=%d' % k, 90)); r = json.loads(get(A + 'FishStockReferencePoints?assessmentKey=%d' % k, 90))
            json.dump({'s': s, 'r': r}, open(f + '.part', 'w')); os.replace(f + '.part', f); return 1
        except Exception as e: return -1
    os.makedirs('ices/a', exist_ok=True)
    with cf.ThreadPoolExecutor(6) as ex: res = list(ex.map(one, keys))
    print('ices anahtar', len(keys), 'indirilen', res.count(1), 'hata', res.count(-1))
elif what == 'obis':
    A = 'https://api.obis.org/v3/'
    areas = [a for a in json.load(open('raw/obis_areas.json'))['results'] if a['type'] in ('lme', 'iho')]
    os.makedirs('obis/c', exist_ok=True)
    def one(arg):
        a, tag, extra = arg; f = 'obis/c/%s_%s.json' % (a['id'], tag)
        if os.path.isfile(f): return 0
        try:
            rows = []; skip = 0
            while True:
                d = json.loads(get(A + 'checklist?areaid=%s&size=10000&skip=%d%s' % (a['id'], skip, extra), 300))
                rs = d.get('results', [])
                rows += [(r.get('taxonID'), r.get('taxonRank'), r.get('records'), r.get('phylum'), r.get('class'), r.get('kingdom')) for r in rs]
                skip += len(rs)
                if len(rs) < 10000 or skip >= d.get('total', 0): break
            st = json.loads(get(A + 'statistics?areaid=%s%s' % (a['id'], extra), 120))
            json.dump({'area': a, 'tag': tag, 'stats': st, 'rows': rows}, open(f + '.part', 'w')); os.replace(f + '.part', f); return 1
        except Exception as e: return -1
    jobs = [(a, tag, extra) for a in areas for tag, extra in (('dated', '&enddate=2026-12-31'), ('all', ''), ('to1999', '&enddate=1999-12-31'), ('to2009', '&enddate=2009-12-31'), ('to2019', '&enddate=2019-12-31'))]
    with cf.ThreadPoolExecutor(4) as ex: res = list(ex.map(one, jobs))
    print('obis is', len(jobs), 'indirilen', res.count(1), 'hata', res.count(-1))
'''

if __name__ == "__main__":
    for _d in DIRS: os.makedirs(_d, exist_ok=True)
    for _m in MODULES:
        with open(_m, 'w', encoding='utf-8') as _f: _f.write(MODULES[_m].lstrip('\n'))
    for _m in ORDER:
        if _m in SKIP_IF_EXISTS and os.path.exists(SKIP_IF_EXISTS[_m]):
            print('== atlandi / skipped (cikti mevcut / output present):', _m); continue
        print('==', _m, flush=True); subprocess.run([sys.executable, _m], check=True, env=dict(os.environ, MPLBACKEND='Agg'))   # her modul ayri surecte calisir / each module runs in its own process
    print('Bitti / done. Sonuclar / results: paper/out/results.json')