# -*- coding: utf-8 -*-

# ------------------------------------------------------------------------------
# # Dünya deniz ticareti hangi deniz seviyesi yükselişini yaşıyor? 2.015 limanda ticaret ağırlıklı göreli deniz seviyesi eğilimi ve gözlem açığı
#
# *What sea-level rise does world seaborne trade experience? Trade-weighted relative sea-level trend at 2,015 ports and the observation gap*
#
# **GOMD Yayınlar · Makale 16** — https://gomd.pirireis.edu.tr/publications.php?id=rp16
#
# Bu not defteri makaledeki **tüm sayıları, tabloları ve şekilleri** üreten kodun Colab sürümüdür. Veri, replikasyon paketinde (12 Kasım 2026 indirmesi) varsa oradan; yoksa özgün sağlayıcılardan (PSMSL, IMF PortWatch, NOAA LSA, NGL, Dünya Bankası) indirilir. Bütün sayılar `out/results.json` dosyasına yazılır.
#
# This notebook reproduces **every number, table and figure** of the paper. Data are taken from the replication package (retrieved 12 November 2026) if present; otherwise they are downloaded from the original providers. All numbers are written to `out/results.json`.
#
# *Makale hazırlık aşamasındadır; dergi adı ve tam metin kabulden sonra eklenecektir. / The paper is in preparation; journal name and full text will be added after acceptance.*
# ------------------------------------------------------------------------------

# ------------------------------------------------------------------------------
# ## 0. Kurulum / Setup
# ------------------------------------------------------------------------------
# [Colab Magic] !pip install -q pandas numpy matplotlib netCDF4 requests

# ------------------------------------------------------------------------------
# ## 1. Veri / Data
# Colab'da birebir yeniden üretim için: GOMD Yayınlar sayfasındaki replikasyon paketini (.zip) yükleyin. / For exact reproduction in Colab: upload the replication package (.zip) from the GOMD Publications page.
# ------------------------------------------------------------------------------
VERI_KAYNAGI = 'kaynak'   # Colab'da 'zip' secenegi icin not defterine bakin / see the notebook for the 'zip' option in Colab
# Veri / Data. Replikasyon paketindeki veri varsa onu kullanır (12 Kasım 2026 indirmesi, birebir yeniden üretim);
# yoksa özgün kaynaktan indirir (güncel sürüm; sağlayıcı güncellemeleri nedeniyle sonuçlar farklılaşabilir).
# Uses the data shipped in the replication package if present (retrieved 12 Nov 2026, exact reproduction);
# otherwise downloads from the original providers (current release; results may differ after provider updates).
import os, json, zipfile
import pandas as pd, numpy as np
def _get(url, path, timeout=600):
    if os.path.exists(path): return
    print('indiriliyor / downloading:', url)
    import urllib.request
    req = urllib.request.Request(url, headers={'User-Agent': 'gomd-replication/1.0'})
    with urllib.request.urlopen(req, timeout=timeout) as r, open(path + '.tmp', 'wb') as f:
        while True:
            b = r.read(1 << 20)
            if not b: break
            f.write(b)
    os.replace(path + '.tmp', path)
for _d in ('data/psmsl', 'data/pw', 'out', 'fig'): os.makedirs(_d, exist_ok=True)
_get('https://psmsl.org/data/obtaining/rlr.monthly.data/rlr_monthly.zip', 'data/psmsl/rlr_monthly.zip')
if not os.path.exists('data/psmsl/rlr_monthly/filelist.txt'): zipfile.ZipFile('data/psmsl/rlr_monthly.zip').extractall('data/psmsl')
_get('https://www.star.nesdis.noaa.gov/socd/lsa/SeaLevelRise/slr/slr_sla_gbl_free_all_66.csv', 'data/gmsl.csv')
_get('https://www.star.nesdis.noaa.gov/socd/lsa/SeaLevelRise/slr/slr_map_ref.nc', 'data/slr_map.nc')
_get('https://geodesy.unr.edu/velocities/midas.IGS14.txt', 'data/midas.txt')
_get('https://coastwatch.pfeg.noaa.gov/erddap/griddap/etopo180.csv?altitude[(-89.75):30:(89.75)][(-179.75):30:(179.75)]', 'data/etopo_30min.csv')
if not os.path.exists('data/wb_country.csv'):
    import urllib.request, csv
    with urllib.request.urlopen('https://api.worldbank.org/v2/country?format=json&per_page=400', timeout=120) as _r: _c = json.load(_r)[1]
    with open('data/wb_country.csv', 'w', newline='') as _f:
        _w = csv.writer(_f); _w.writerow(['iso3', 'name', 'region', 'income'])
        for _x in _c:
            if _x['region']['value'] != 'Aggregates': _w.writerow([_x['id'], _x['name'], _x['region']['value'], _x['incomeLevel']['value']])
# IMF PortWatch: liman veritabanı ve liman–yıl toplamları (2.065 liman; birkaç on dakika sürebilir) / ports database and port–year totals (may take tens of minutes)
if not (os.path.exists('data/pw/ports_global.csv') and os.path.exists('data/pw/port_year_global.csv')):
    import requests, json, time, csv, sys
    from concurrent.futures import ThreadPoolExecutor
    S="https://services9.arcgis.com/weJ1QsnbMYJlCHdG/arcgis/rest/services/"
    # ports db
    feats=[]
    for off in (0,1000,2000):
        r=requests.get(S+"PortWatch_ports_database/FeatureServer/0/query",params={"where":"1=1","outFields":"*","returnGeometry":"false","f":"json","resultOffset":off,"resultRecordCount":1000,"orderByFields":"ObjectId"},timeout=120).json()
        feats+= [f['attributes'] for f in r['features']]
    print(len(feats)); 
    keys=list(feats[0].keys())
    with open('data/pw/ports_global.csv','w',newline='',encoding='utf-8') as f:
        w=csv.DictWriter(f,keys); w.writeheader(); w.writerows(feats)
    F=["portcalls","import","export","import_container","export_container","import_tanker","export_tanker","import_dry_bulk","export_dry_bulk","portcalls_container","portcalls_tanker","portcalls_dry_bulk"]
    st=[{"statisticType":"sum","onStatisticField":f,"outStatisticFieldName":"s%d"%i} for i,f in enumerate(F)]+[{"statisticType":"count","onStatisticField":"portcalls","outStatisticFieldName":"n"}]
    B=S+"Daily_Ports_Data/FeatureServer/0/query"
    def one(pid):
        for k in range(5):
            try:
                d=requests.get(B,params={"where":f"portid='{pid}'","outStatistics":json.dumps(st),"groupByFieldsForStatistics":"year","f":"json"},timeout=120).json()
                if 'features' in d:
                    return [(pid,a['attributes']['year'],a['attributes']['n'])+tuple(a['attributes']['s%d'%i] for i in range(len(F))) for a in d['features']]
            except Exception as e: pass
            time.sleep(2+3*k)
        return [(pid,'ERR')]
    ids=[f['portid'] for f in feats]
    rows=[]
    with ThreadPoolExecutor(8) as ex:
        for i,r in enumerate(ex.map(one,ids)):
            rows+=r
            if i%200==0: print(i,flush=True)
    with open('data/pw/port_year_global.csv','w',newline='') as f:
        w=csv.writer(f); w.writerow(["portid","year","n"]+F); w.writerows(rows)
    print('done',len(rows),sum(1 for r in rows if r[1]=='ERR'))

# ------------------------------------------------------------------------------
# ## 2. Yardımcı işlevler / Helper functions
# ------------------------------------------------------------------------------
# -*- coding: utf-8 -*-
import numpy as np, pandas as pd, os, glob
_DATA='data'
def stations():
    rows=[]
    for l in open(_DATA+'/psmsl/rlr_monthly/filelist.txt',encoding='latin-1'):
        p=[x.strip() for x in l.split(';')]
        rows.append(dict(id=int(p[0]),lat=float(p[1]),lon=float(p[2]),name=p[3],coast=p[4],flag=p[6]))
    return pd.DataFrame(rows).set_index('id')
def series(i):
    a=np.loadtxt(f'{_DATA}/psmsl/rlr_monthly/data/{i}.rlrdata',delimiter=';',usecols=(0,1,3),ndmin=2)
    t,v,fl=a[:,0],a[:,1].astype(float),a[:,2]
    v[v==-99999]=np.nan
    v[(fl.astype(int)%10)==1]=np.nan  # dikkat bayragi (flag for attention) -> haric
    return t,v
def trend(t,v,y0,y1,minfrac=0.7,accel=False):
    """Mevsimselligi arindirilmis aylik seri uzerinde EKK egimi (mm/yil), AR(1) duzeltmeli standart hata."""
    m=(t>=y0)&(t<y1+1); t=t[m]; v=v[m]
    nexp=int(round((y1+1-y0)*12)); ok=~np.isnan(v)
    if ok.sum()<minfrac*nexp: return None
    # kaydin donemin basini ve sonunu kapsamasi (en az 0.85 yayilim)
    if (t[ok].max()-t[ok].min())<0.85*(y1+1-y0): return None
    mon=np.round((t-np.floor(t))*12-0.5).astype(int)
    x=t-(y0+y1+1)/2
    cols=[np.ones_like(x),x]+([x**2] if accel else [])
    for k in range(1,12): cols.append((mon==k).astype(float))
    X=np.column_stack(cols)[ok]; y=v[ok]
    b,*_=np.linalg.lstsq(X,y,rcond=None); r=y-X@b
    n,p=X.shape; s2=r@r/(n-p); XtXi=np.linalg.inv(X.T@X)
    # AR(1): ardisik aylar
    tt=t[ok]; adj=np.isclose(np.diff(tt),1/12,atol=0.01)
    rho=np.corrcoef(r[:-1][adj],r[1:][adj])[0,1] if adj.sum()>20 else 0.0
    rho=min(max(rho,0),0.95); infl=(1+rho)/(1-rho)
    se=np.sqrt(s2*np.diag(XtXi)*infl)
    out=dict(trend=b[1],se=se[1],n=int(n),rho=rho,first=float(tt.min()),last=float(tt.max()))
    if accel: out.update(acc=2*b[2],acc_se=2*se[2])
    return out
def hav(lat1,lon1,lat2,lon2):
    p=np.pi/180; a=np.sin((lat2-lat1)*p/2)**2+np.cos(lat1*p)*np.cos(lat2*p)*np.sin((lon2-lon1)*p/2)**2
    return 2*6371.0*np.arcsin(np.sqrt(a))

# ------------------------------------------------------------------------------
# ## 3. Mareograf eğilimleri / Tide-gauge trends
# ------------------------------------------------------------------------------
st=stations(); rows=[]
for i in st.index:
    t,v=series(i); ok=~np.isnan(v)
    r=dict(id=i,n_all=int(ok.sum()),first=float(t[ok].min()) if ok.any() else np.nan,last=float(t[ok].max()) if ok.any() else np.nan)
    for tag,(y0,y1) in dict(sat=(1993,2024),s2=(2000,2024),long=(1970,2024),early=(1970,1992)).items():
        o=trend(t,v,y0,y1)
        if o: r.update({f'{tag}_tr':o['trend'],f'{tag}_se':o['se'],f'{tag}_n':o['n']})
    o=trend(t,v,1970,2024,accel=True)
    if o: r.update(acc=o['acc'],acc_se=o['acc_se'])
    rows.append(r)
d=st.join(pd.DataFrame(rows).set_index('id')); d.to_csv('out/gauge_trends.csv')
print(len(d), d[['sat_tr','s2_tr','long_tr','early_tr','acc']].describe().round(2).to_string())
print((d['last']>=2020).sum(), (d['last']>=2023).sum(), d.flag.value_counts().to_dict())

# ------------------------------------------------------------------------------
# ## 4. Eşleştirme, gözlem açığı ve ticaret ağırlıklı eğilim / Matching, observation gap and trade-weighted trend
# ------------------------------------------------------------------------------
# -*- coding: utf-8 -*-
import warnings; warnings.filterwarnings('ignore')
import json, sys
rng=np.random.default_rng(20261112); R={}
G=pd.read_csv('out/gauge_trends.csv',index_col=0)
P=pd.read_csv('data/pw/ports_global.csv'); PY=pd.read_csv('data/pw/port_year_global.csv')
PY=PY[PY.year!='ERR'] if PY.year.dtype==object else PY; PY['year']=PY.year.astype(int)
full=PY[(PY.year>=2019)&(PY.year<=2025)&(PY.n>=360)]
w=full.groupby('portid').agg(W=('import',lambda x:0)).drop(columns='W')
full['trade']=full['import']+full['export']; full['cont']=full.import_container+full.export_container; full['tank']=full.import_tanker+full.export_tanker; full['bulk']=full.import_dry_bulk+full.export_dry_bulk
w=full.groupby('portid')[['trade','cont','tank','bulk','portcalls','import','export']].mean()
P=P.merge(w,left_on='portid',right_index=True,how='left').fillna({'trade':0,'cont':0,'tank':0,'bulk':0,'portcalls':0})
wb=pd.read_csv('data/wb_country.csv'); P=P.merge(wb[['iso3','region','income']],left_on='ISO3',right_on='iso3',how='left')
P['region']=P.region.fillna('Diğer'); P['income']=P.income.fillna('Sınıflanmamış')
R['n_ports']=len(P); R['n_ports_trade']=int((P.trade>0).sum()); R['world_trade_mt']=float(P.trade.sum()/1e6); R['world_calls']=float(P.portcalls.sum())
# GMSL 1993-2024
g=pd.read_csv('data/gmsl.csv',comment='#'); g['v']=g.iloc[:,1:].mean(axis=1); g=g[(g.year>=1993)&(g.year<2025)]
b=np.polyfit(g.year,g.v,1); R['gmsl_9324']=float(b[0]); g2=g[g.year>=2000]; R['gmsl_0024']=float(np.polyfit(g2.year,g2.v,1)[0])
GM=R['gmsl_9324']
# GPS
m=pd.read_csv('data/midas.txt',sep=r'\s+',header=None,usecols=[0,4,10,13,24,25],names=['sta','dur','vu','svu','lat','lon'])
m['lon']=((m.lon+180)%360)-180; m=m[(m.dur>=8)&(m.svu<0.001)]; R['n_gps']=len(m)
def nearest(lat,lon,LAT,LON):
    D=hav(lat[:,None],lon[:,None],LAT[None,:],LON[None,:]); j=D.argmin(1); return j,D[np.arange(len(lat)),j]
j,dd=nearest(G.lat.values,G.lon.values,m.lat.values,m.lon.values)
G['gps_d']=dd; G['vlm']=np.where(dd<=30,m.vu.values[j]*1000,np.nan); G['vlm_se']=np.where(dd<=30,m.svu.values[j]*1000,np.nan)
Gs=G[G.sat_tr.notna()]; R['n_gauges_all']=len(G); R['n_gauges_sat']=len(Gs); R['n_gauges_recent']=int((G['last']>=2020).sum())
R['gauge_mean']=float(Gs.sat_tr.mean()); R['gauge_median']=float(Gs.sat_tr.median()); R['gauge_sd']=float(Gs.sat_tr.std())
v=Gs[Gs.vlm.notna()]; R['gps_pairs']=len(v); R['gps_corr']=float(np.corrcoef(v.sat_tr,-v.vlm)[0,1]); R['gps_slope']=float(np.polyfit(-v.vlm,v.sat_tr,1)[0]); R['gps_mean_geoc']=float((v.sat_tr+v.vlm).mean()); R['gps_sd_geoc']=float((v.sat_tr+v.vlm).std()); R['gps_sd_rel']=float(v.sat_tr.std())
# eslestirme
def match(P,Gx,col,rad):
    j,d=nearest(P.lat.values,P.lon.values,Gx.lat.values,Gx.lon.values)
    ok=d<=rad; return np.where(ok,Gx[col].values[j],np.nan),np.where(ok,Gx.index.values[j],-1),d
for rad in (50,100,200):
    P[f'tr{rad}'],P[f'gid{rad}'],P[f'd{rad}']=match(P,Gs,'sat_tr',rad)
P['se100']=match(P,Gs,'sat_se',100)[0]
P['any100']=match(P,G,'lat',100)[0]; P['rec100']=match(P,G[G['last']>=2020],'lat',100)[0]; P['any25']=match(P,G,'lat',25)[0]; P['rec25']=match(P,G[G['last']>=2020],'lat',25)[0]
P['vlm100']=match(P,Gs,'vlm',100)[0]
# idw ortalamasi (100 km)
def idw(P,Gx,rad=100):
    D=hav(P.lat.values[:,None],P.lon.values[:,None],Gx.lat.values[None,:],Gx.lon.values[None,:]); W=np.where(D<=rad,1/np.maximum(D,5),0); s=W.sum(1)
    return np.where(s>0,(W*Gx.sat_tr.values[None,:]).sum(1)/np.where(s>0,s,1),np.nan)
P['idw100']=idw(P,Gs)
P.to_csv('out/ports_matched.csv',index=False)
T=P[P.trade>0].copy()
def sh(mask,wcol='trade',df=T): return float(df.loc[mask,wcol].sum()/df[wcol].sum())
cov={}
for k,col in dict(any25='any25',rec25='rec25',any100='any100',rec100='rec100',sat50='tr50',sat100='tr100',sat200='tr200').items():
    cov[k]=dict(ports=float(T[col].notna().mean()),trade=sh(T[col].notna()),calls=sh(T[col].notna(),'portcalls'),cont=sh(T[col].notna(),'cont'))
R['coverage']=cov
def covtab(by):
    o={}
    for k,s in T.groupby(by):
        o[k]=dict(n=len(s),trade_share=float(s.trade.sum()/T.trade.sum()),any100=float(s.loc[s.any100.notna(),'trade'].sum()/s.trade.sum()),rec100=float(s.loc[s.rec100.notna(),'trade'].sum()/s.trade.sum()),sat100=float(s.loc[s.tr100.notna(),'trade'].sum()/s.trade.sum()),ports_sat100=float(s.tr100.notna().mean()))
    return o
R['cov_region']=covtab('region'); R['cov_income']=covtab('income')
# en buyuk gozlemsiz limanlar
u=T[T.tr100.isna()].sort_values('trade',ascending=False).head(15); R['top_unobserved']=[dict(port=r.portname,iso=r.ISO3,mt=r.trade/1e6,any=bool(pd.notna(r.any100)),rec=bool(pd.notna(r.rec100))) for r in u.itertuples()]
# endeks
def wmean(x,w): return float((x*w).sum()/w.sum())
def wq(x,w,q):
    o=np.argsort(x); x=x[o]; w=w[o]; c=np.cumsum(w)/w.sum(); return float(x[np.searchsorted(c,q)])
def index(df,col='tr100',wcol='trade',B=2000,secol=None):
    s=df[df[col].notna()&(df[wcol]>0)]; x=s[col].values; ww=s[wcol].values; n=len(s)
    est=wmean(x,ww); bs=[]
    for _ in range(B):
        i=rng.integers(0,n,n); xx=x[i]+(rng.normal(0,s[secol].values[i]) if secol else 0); bs.append(wmean(xx,ww[i]))
    return dict(n=n,est=est,lo=float(np.percentile(bs,2.5)),hi=float(np.percentile(bs,97.5)),eq=float(x.mean()),med_w=wq(x,ww,0.5),q25=wq(x,ww,0.25),q75=wq(x,ww,0.75),q90=wq(x,ww,0.9),
        sh_gt_gm=float(ww[x>GM].sum()/ww.sum()),sh_gt5=float(ww[x>5].sum()/ww.sum()),sh_gt2gm=float(ww[x>2*GM].sum()/ww.sum()),sh_gt10=float(ww[x>10].sum()/ww.sum()),sh_neg=float(ww[x<0].sum()/ww.sum()))
R['tw']=index(T,secol='se100')
R['tw_by_cargo']={k:index(T,wcol=k) for k in ('cont','tank','bulk','portcalls','import','export')}
R['tw_region']={k:index(s,B=500) for k,s in T.groupby('region') if s.tr100.notna().sum()>=8}
R['tw_income']={k:index(s,B=500) for k,s in T.groupby('income') if s.tr100.notna().sum()>=8}
# ayristirma: olcek (gmsl) + konum (liman esit - gmsl) + yogunlasma (tw - esit)
R['decomp']=dict(gmsl=GM,gauge_eq=R['gauge_mean'],port_eq=R['tw']['eq'],tw=R['tw']['est'],location=R['tw']['eq']-GM,concentration=R['tw']['est']-R['tw']['eq'])
# katki: w_i (tr_i - GM)
s=T[T.tr100.notna()].copy(); s['ws']=s.trade/s.trade.sum(); s['contr']=s.ws*(s.tr100-GM); s=s.sort_values('contr',ascending=False)
R['top_contrib']=[dict(port=r.portname,iso=r.ISO3,mt=r.trade/1e6,tr=r.tr100,d=r.d100,gid=int(r.gid100),gname=G.loc[int(r.gid100),'name'],contr=r.contr,share=r.ws) for r in s.head(15).itertuples()]
R['bottom_contrib']=[dict(port=r.portname,iso=r.ISO3,mt=r.trade/1e6,tr=r.tr100,gname=G.loc[int(r.gid100),'name'],contr=r.contr) for r in s.tail(8).itertuples()]
R['excess_total']=float(s.contr.sum()); R['excess_top10_share']=float(s.contr.head(10).sum()/s.contr.sum()); R['pos_contr']=float(s.contr[s.contr>0].sum()); R['neg_contr']=float(s.contr[s.contr<0].sum())
big=T.sort_values('trade',ascending=False).head(30); R['top30']=[dict(port=r.portname,iso=r.ISO3,mt=r.trade/1e6,tr=(None if pd.isna(r.tr100) else r.tr100),se=(None if pd.isna(r.se100) else r.se100),d=r.d100,gname=(G.loc[int(r.gid100),'name'] if r.gid100>0 else None)) for r in big.itertuples()]
R['top30_obs']=int(big.tr100.notna().sum())
# VLM alt kumesi
vv=T[T.vlm100.notna()&T.tr100.notna()]; R['vlm']=dict(n=len(vv),trade_share=sh(T.vlm100.notna()&T.tr100.notna()),tw_rel=wmean(vv.tr100.values,vv.trade.values),tw_vlm=wmean(vv.vlm100.values,vv.trade.values),tw_geoc=wmean((vv.tr100+vv.vlm100).values,vv.trade.values),eq_vlm=float(vv.vlm100.mean()))
# gozlemsiz limanlar icin sinirlar
reg=T[T.tr100.notna()].groupby('region').apply(lambda s:wmean(s.tr100.values,s.trade.values))
allw=T.trade.sum(); obs=T.tr100.notna()
fill_gm=T.tr100.fillna(GM); fill_reg=T.tr100.fillna(T.region.map(reg)).fillna(R['tw']['est'])
R['bounds']=dict(fill_gmsl=wmean(fill_gm.values,T.trade.values),fill_region=wmean(fill_reg.values,T.trade.values),need_unobs=float((GM*allw-(T.tr100[obs]*T.trade[obs]).sum())/T.trade[~obs].sum()))
# saglamlik
rb={}
for rad in (50,200): rb[f'rad{rad}']=index(T,col=f'tr{rad}',B=300)
rb['idw100']=index(T,col='idw100',B=300)
T2=T.copy(); T2['trN']=match(T2,Gs[Gs.flag=='N'],'sat_tr',100)[0]; rb['noflag']=index(T2,col='trN',B=300)
x=T[T.tr100.notna()]; lo,hi=np.percentile(x.tr100,[2.5,97.5]); T2['trW']=T.tr100.clip(lo,hi); rb['winsor']=index(T2,col='trW',B=300)
T2['tr00']=match(T2,G[G.s2_tr.notna()],'s2_tr',100)[0]; rb['p2000_2024']=index(T2,col='tr00',B=300)
# en buyuk katkili ulke haric
for iso in ('CHN','USA'): rb['ex_'+iso]=index(T[T.ISO3!=iso],B=300)
# ulke basina tek agirlik? liman basina ust sinir: en buyuk 10 liman haric
rb['ex_top10ports']=index(T[~T.portid.isin(big.portid.head(10))],B=300)
# sadece <=25 km
T2['tr25']=np.where(T.d100<=25,T.tr100,np.nan); rb['rad25']=index(T2,col='tr25',B=300)
R['robust']={k:{a:v[a] for a in ('n','est','lo','hi','eq')} for k,v in rb.items()}
# ivme: ayni olceklerde erken donem (1970-1992) ve uydu donemi, iki donemi de olan gostergeler
Gb=G[G.early_tr.notna()&G.sat_tr.notna()]; T2['e']=match(T2,Gb,'early_tr',100)[0]; T2['l']=match(T2,Gb,'sat_tr',100)[0]; q=T2[T2.e.notna()&(T2.trade>0)]
R['accel']=dict(n_g=len(Gb),n_p=len(q),trade_share=float(q.trade.sum()/T.trade.sum()),tw_early=wmean(q.e.values,q.trade.values),tw_late=wmean(q.l.values,q.trade.values),g_early=float(Gb.early_tr.mean()),g_late=float(Gb.sat_tr.mean()))
# gosterge aginin yaslanmasi: yillara gore veri gonderen istasyon sayisi
yrs={}
for i in G.index:
    t,vv_=series(i); ok=~np.isnan(vv_)
    for y in np.unique(np.floor(t[ok]).astype(int)): yrs[y]=yrs.get(y,0)+1
R['stations_by_year']={int(k):int(v) for k,v in sorted(yrs.items()) if k>=1950}
# 2050 dogrusal uzatim (aciklama amacli)
R['extrap']=dict(tw_cm=R['tw']['est']*25/10,gm_cm=GM*25/10,q90_cm=R['tw']['q90']*25/10)
json.dump(R,open('out/results.json','w'),ensure_ascii=False,indent=1,default=float)
for k,vv_ in R.items():
    if k in ('stations_by_year',): continue
    print(k, json.dumps(vv_,ensure_ascii=False,default=float)[:1500])

# ------------------------------------------------------------------------------
# ## 5. Altimetre + GNSS, doğrulama ve bileşik tahmin / Altimetry + GNSS, validation and composite estimate
# ------------------------------------------------------------------------------
# -*- coding: utf-8 -*-
import warnings; warnings.filterwarnings('ignore')
import json, sys, netCDF4
R=json.load(open('out/results.json')); GM=R['gmsl_9324']; rng=np.random.default_rng(7)
G=pd.read_csv('out/gauge_trends.csv',index_col=0); P=pd.read_csv('out/ports_matched.csv')
d=netCDF4.Dataset('data/slr_map.nc'); lon=d['lon'][:].data; lat=d['lat'][:].data; Z=np.ma.filled(d['Sea_level_trends'][:].astype(float),np.nan); Z[Z==0]=np.nan
LO,LA=np.meshgrid(lon,lat); ok=~np.isnan(Z); alo,ala,az=LO[ok],LA[ok],Z[ok]
aw=np.cos(np.deg2rad(ala)); R['alt_mean']=float((az*aw).sum()/aw.sum()); R['alt_sd']=float(np.sqrt(((az-R['alt_mean'])**2*aw).sum()/aw.sum()))
def alt_at(la,lo,rad=150):
    out=np.full(len(la),np.nan); dist=np.full(len(la),np.nan)
    for i in range(len(la)):
        m=(np.abs(ala-la[i])<2.5); D=hav(la[i],lo[i],ala[m],alo[m]); 
        if len(D)==0: continue
        k=D<=rad
        if k.any(): w=1/np.maximum(D[k],25); out[i]=(az[m][k]*w).sum()/w.sum(); dist[i]=D.min()
    return out,dist
m=pd.read_csv('data/midas.txt',sep=r'\s+',header=None,usecols=[0,4,10,13,24,25],names=['sta','dur','vu','svu','lat','lon'])
m['lon']=((m.lon+180)%360)-180; m=m[(m.dur>=8)&(m.svu<0.001)]
def gps_at(la,lo,rad=30):
    out=np.full(len(la),np.nan); n=np.zeros(len(la),int)
    for i in range(len(la)):
        k=(np.abs(m.lat.values-la[i])<0.5); 
        if not k.any(): continue
        D=hav(la[i],lo[i],m.lat.values[k],m.lon.values[k]); s=D<=rad
        if s.any(): out[i]=np.median(m.vu.values[k][s])*1000; n[i]=s.sum()
    return out,n
G['alt'],_=alt_at(G.lat.values,G.lon.values); G['vlm_med'],G['ngps']=gps_at(G.lat.values,G.lon.values)
Gs=G[G.sat_tr.notna()&G.alt.notna()]
v=Gs[Gs.vlm_med.notna()]; hyb=v.alt-v.vlm_med
def stats(a,b): e=a-b; return dict(n=len(a),corr=float(np.corrcoef(a,b)[0,1]),bias=float(e.mean()),rmse=float(np.sqrt((e**2).mean())),mad=float(np.median(np.abs(e-np.median(e)))),slope=float(np.polyfit(b,a,1)[0]))
R['valid']=dict(hybrid=stats(v.sat_tr.values,hyb.values),alt_only_same=stats(v.sat_tr.values,v.alt.values),alt_only_all=stats(Gs.sat_tr.values,Gs.alt.values),gmsl_same=dict(rmse=float(np.sqrt(((v.sat_tr-GM)**2).mean()))))
# guvenilir alt kume: |fark|<... degil; dogrudan raporla. Liman duzeyi
P['alt'],P['alt_d']=alt_at(P.lat.values,P.lon.values); P['vlm_p'],P['ngps']=gps_at(P.lat.values,P.lon.values)
P['hyb']=P.alt-P.vlm_p
T=P[P.trade>0].copy(); W=T.trade.sum()
def sh(mask): return float(T.loc[mask,'trade'].sum()/W)
def wm(x,w): return float((x*w).sum()/w.sum())
def boot(x,w,B=2000):
    n=len(x); b=[wm(x[i],w[i]) for i in (rng.integers(0,n,n) for _ in range(B))]; return float(np.percentile(b,2.5)),float(np.percentile(b,97.5))
T['tier']=np.where(T.tr100.notna(),'gauge',np.where(T.hyb.notna(),'gps',np.where(T.alt.notna(),'alt','none')))
T['best']=np.where(T.tr100.notna(),T.tr100,np.where(T.hyb.notna(),T.hyb,T.alt))
R['tiers']={k:dict(ports=int((T.tier==k).sum()),trade=sh(T.tier==k)) for k in ('gauge','gps','alt','none')}
R['gps_cov']=dict(ports=float(T.vlm_p.notna().mean()),trade=sh(T.vlm_p.notna()),alt_trade=sh(T.alt.notna()))
def est(mask,col):
    s=T[mask&T[col].notna()]; lo,hi=boot(s[col].values,s.trade.values,800); return dict(n=len(s),trade=float(s.trade.sum()/W),est=wm(s[col].values,s.trade.values),lo=lo,hi=hi,eq=float(s[col].mean()))
allm=np.ones(len(T),bool)
R['hyb_est']=dict(alt_all=est(allm,'alt'),hyb_all=est(allm,'hyb'),gauge_or_hyb=est(T.tier.isin(['gauge','gps']),'best'),best_all=est(allm,'best'),
   gauge_ports_alt=est(T.tr100.notna(),'alt'),gauge_ports_gauge=est(T.tr100.notna(),'tr100'),nogauge_alt=est(T.tr100.isna(),'alt'),nogauge_hyb=est(T.tr100.isna(),'hyb'),
   both_gauge=est(T.tr100.notna()&T.hyb.notna(),'tr100'),both_hyb=est(T.tr100.notna()&T.hyb.notna(),'hyb'))
# VLM dagilimi (GPS'li limanlar)
s=T[T.vlm_p.notna()]; R['vlm_ports']=dict(n=len(s),tw=wm(s.vlm_p.values,s.trade.values),eq=float(s.vlm_p.mean()),sh_sub1=float(s.loc[s.vlm_p<-1,'trade'].sum()/s.trade.sum()),sh_sub3=float(s.loc[s.vlm_p<-3,'trade'].sum()/s.trade.sum()),sh_up1=float(s.loc[s.vlm_p>1,'trade'].sum()/s.trade.sum()))
# gozlemsiz ticaretin ne kadar batmasi gerekir ki genel ortalama X olsun: alt-only limanlarda ortalama cokme senaryolari
ao=T[T.tier=='alt']; base=(T.best*T.trade).sum()
R['scen']={str(k):float((base+k*ao.trade.sum())/T.loc[T.best.notna(),'trade'].sum()) for k in (0,1,2,3,5)}
R['alt_only_share']=float(ao.trade.sum()/W)
# bolge / gelir: en iyi tahmin ve katman paylari
def grp(by):
    o={}
    for k,s in T.groupby(by):
        if s.best.notna().sum()<8: continue
        o[k]=dict(n=len(s),tsh=float(s.trade.sum()/W),best=wm(s.best[s.best.notna()].values,s.trade[s.best.notna()].values),alt=wm(s.alt[s.alt.notna()].values,s.trade[s.alt.notna()].values),
           gauge=float(s.loc[s.tier=='gauge','trade'].sum()/s.trade.sum()),gps=float(s.loc[s.tier=='gps','trade'].sum()/s.trade.sum()),altonly=float(s.loc[s.tier=='alt','trade'].sum()/s.trade.sum()))
    return o
R['grp_region']=grp('region'); R['grp_income']=grp('income')
big=T.sort_values('trade',ascending=False).head(25)
R['top25']=[dict(port=r.portname,iso=r.ISO3,mt=r.trade/1e6,tier=r.tier,gauge=None if pd.isna(r.tr100) else r.tr100,alt=None if pd.isna(r.alt) else r.alt,vlm=None if pd.isna(r.vlm_p) else r.vlm_p,best=None if pd.isna(r.best) else r.best) for r in big.itertuples()]
T.to_csv('out/ports_final.csv',index=False); G.to_csv('out/gauges_final.csv')
json.dump(R,open('out/results.json','w'),ensure_ascii=False,indent=1,default=float)
for k in ('alt_mean','alt_sd','valid','tiers','gps_cov','hyb_est','vlm_ports','scen','alt_only_share','grp_region','grp_income','top25'): print(k,json.dumps(R[k],ensure_ascii=False,default=float)[:2600])

# ------------------------------------------------------------------------------
# ## 6. Ülke tablosu ve ek ölçüler / Country table and additional measures
# ------------------------------------------------------------------------------
import json, numpy as np, pandas as pd
R=json.load(open('out/results.json')); T=pd.read_csv('out/ports_final.csv'); GM=R['gmsl_9324']; W=T.trade.sum()
def wm(x,w): return float((x*w).sum()/w.sum())
c=[]
for iso,s in T.groupby('ISO3'):
    o=s[s.tr100.notna()]; b=s[s.best.notna()]
    c.append(dict(iso=iso,n=len(s),tsh=float(s.trade.sum()/W),gauge=float(o.trade.sum()/s.trade.sum()),gps=float(s.loc[s.tier=='gps','trade'].sum()/s.trade.sum()),alt=float(s.loc[s.tier=='alt','trade'].sum()/s.trade.sum()),
       tw=(wm(o.tr100.values,o.trade.values) if len(o)>=3 else None),n_obs=len(o),alt_tw=(wm(s.alt[s.alt.notna()].values,s.trade[s.alt.notna()].values) if s.alt.notna().any() else None)))
c=sorted(c,key=lambda r:-r['tsh']); R['country']=c[:20]; R['n_countries']=len(c)
R['n_countries_zero_gauge']=int(sum(1 for r in c if r['gauge']==0)); R['zero_gauge_trade']=float(sum(r['tsh'] for r in c if r['gauge']==0))
s=T.sort_values('trade',ascending=False); cs=s.trade.cumsum()/W
R['conc']=dict(top10=float(cs.iloc[9]),top50=float(cs.iloc[49]),top100=float(cs.iloc[99]),top200=float(cs.iloc[199]),n_half=int((cs<0.5).sum()+1))
o=T[T.tr100.notna()]; z=(o.tr100-GM)/o.se100
R['signif']=dict(above=float(o.loc[z>1.96,'trade'].sum()/o.trade.sum()),below=float(o.loc[z<-1.96,'trade'].sum()/o.trade.sum()),med_se=float(o.se100.median()))
# mareograf-liman uzakligi
R['dist']=dict(med=float(o.d100.median()),w_med=float(o.sort_values('d100').pipe(lambda d:d.d100.values[np.searchsorted(d.trade.cumsum().values/d.trade.sum(),0.5)])),sh25=float(o.loc[o.d100<=25,'trade'].sum()/o.trade.sum()))
# ayni mareografi paylasan liman sayisi
R['n_gauges_used']=int(o.gid100.nunique())
# gozlem aciginin yogunlasmasi: mareograf katmani disindaki en buyuk limanlarin hacim payi
u=T[T.tr100.isna()].sort_values('trade',ascending=False); a_=T[T.tier=='alt'].sort_values('trade',ascending=False)
R['close']=dict(g10=float(u.trade.head(10).sum()/W),g20=float(u.trade.head(20).sum()/W),g50=float(u.trade.head(50).sum()/W),g100=float(u.trade.head(100).sum()/W),a20=float(a_.trade.head(20).sum()/W),a50=float(a_.trade.head(50).sum()/W),a100=float(a_.trade.head(100).sum()/W),n_unobs=len(u),n_alt=len(a_))
json.dump(R,open('out/results.json','w'),ensure_ascii=False,indent=1,default=float)
for k in ('country','n_countries','n_countries_zero_gauge','zero_gauge_trade','conc','signif','dist','n_gauges_used','extrap'): print(k,json.dumps(R[k],ensure_ascii=False,default=lambda x:round(float(x),3))[:3000])

# ------------------------------------------------------------------------------
# ## 7. Şekiller / Figures
# ------------------------------------------------------------------------------
# -*- coding: utf-8 -*-
import warnings; warnings.filterwarnings('ignore')
import json, sys, matplotlib; matplotlib.use('Agg'); import matplotlib.pyplot as plt
import numpy as np, pandas as pd
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':9,'axes.spines.top':False,'axes.spines.right':False,'savefig.dpi':300})
C=dict(gauge='#1b4f72',gps='#1e8449',alt='#b03a2e',none='#999999')
R=json.load(open('out/results.json')); T=pd.read_csv('out/ports_final.csv'); G=pd.read_csv('out/gauges_final.csv',index_col=0); GM=R['gmsl_9324']
e=pd.read_csv('data/etopo_30min.csv',skiprows=[1]); e=e[(e.altitude>=0)]; e=e.iloc[::2]
# S1 harita
fig,ax=plt.subplots(figsize=(7.4,3.9))
ax.scatter(e.longitude,e.latitude,s=.25,color='#dcdcdc',marker='s',lw=0,rasterized=True)
lab=dict(gauge='Mareograf (≤100 km, 1993–2024 kaydı)',gps='Uydu altimetresi + GNSS',alt='Yalnız uydu altimetresi (kara hareketi bilinmiyor)',none='Veri yok')
for k in ('alt','gps','gauge','none'):
    s=T[T.tier==k]; ax.scatter(s.lon,s.lat,s=3+s.trade/2.5e6,color=C[k],alpha=.6,lw=0,label=f"{lab[k]}: ticaretin %{100*R['tiers'][k]['trade']:.1f}".replace('.',','))
ax.set_xlim(-180,180); ax.set_ylim(-80,80); ax.set_xticks([]); ax.set_yticks([])
for s in ax.spines.values(): s.set_visible(False)
lg=ax.legend(frameon=False,fontsize=7.5,loc='lower left',bbox_to_anchor=(0,-0.02),markerscale=1)
for h in lg.legend_handles: h.set_sizes([25])
fig.tight_layout(); fig.savefig('fig/Figure_1.png'); plt.close()
# S2 katman paylari gelir + bolge
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1),gridspec_kw={'width_ratios':[1,1.6]})
inc=[('High income','Yüksek gelir'),('Upper middle income','Üst-orta gelir'),('Lower middle income','Alt-orta gelir'),('Low income','Düşük gelir')]
reg=[('North America','Kuzey Amerika'),('Europe & Central Asia','Avrupa ve Orta Asya'),('East Asia & Pacific','Doğu Asya ve Pasifik'),('South Asia','Güney Asya'),('Latin America & Caribbean ','Latin Amerika ve Karayipler'),('Sub-Saharan Africa ','Sahra Altı Afrika'),('Middle East, North Africa, Afghanistan & Pakistan','Orta Doğu ve Kuzey Afrika')]
for a,(lst,key,ttl) in zip(ax,[(inc,'grp_income','(a) Gelir grubu'),(reg,'grp_region','(b) Bölge')]):
    y=np.arange(len(lst))[::-1]; left=np.zeros(len(lst))
    for k,l in (('gauge','Mareograf'),('gps','Altimetre + GNSS'),('altonly','Yalnız altimetre')):
        v=np.array([100*R[key][c][k] for c,_ in lst]); a.barh(y,v,left=left,color=C['alt' if k=='altonly' else k],label=l,height=.65); left+=v
    a.set_yticks(y); a.set_yticklabels([f"{n}\n(%{100*R[key][c]['tsh']:.1f})".replace('.',',') for c,n in lst],fontsize=7.5); a.set_xlim(0,100); a.set_xlabel('Liman ticaretindeki pay (%)'); a.set_title(ttl,fontsize=9,loc='left')
ax[0].legend(frameon=False,fontsize=7.5,loc='upper center',bbox_to_anchor=(0.3,-0.2),ncol=1)
fig.tight_layout(); fig.savefig('fig/Figure_2.png'); plt.close()
# S3 dogrulama
v=G[G.sat_tr.notna()&G.alt.notna()&G.vlm_med.notna()]
fig,ax=plt.subplots(1,2,figsize=(7.4,3.3))
for a,x,t,k in ((ax[0],v.alt,'(a) Yalnız uydu altimetresi','alt_only_same'),(ax[1],v.alt-v.vlm_med,'(b) Altimetre − GNSS düşey hızı','hybrid')):
    a.scatter(x,v.sat_tr,s=6,color='#1b4f72',alpha=.5,lw=0); a.plot([-12,16],[-12,16],color='#999',lw=.8,ls='--'); a.set_xlim(-12,16); a.set_ylim(-12,16)
    a.set_xlabel('Tahmin edilen göreli eğilim (mm/yıl)'); a.set_ylabel('Mareograf eğilimi, 1993–2024 (mm/yıl)'); a.set_title(t,fontsize=9,loc='left')
    a.text(0.04,0.93,f"r = {R['valid'][k]['corr']:.2f}; KOKH = {R['valid'][k]['rmse']:.2f} mm/yıl".replace('.',','),transform=a.transAxes,fontsize=8)
fig.tight_layout(); fig.savefig('fig/Figure_3.png'); plt.close()
# S4 agirlikli dagilim
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1))
s=T[T.tr100.notna()].sort_values('tr100'); ax[0].plot(s.tr100,100*s.trade.cumsum()/s.trade.sum(),color=C['gauge'],lw=1.5,label='Liman ticareti (mareograflı limanlar)')
g=G[G.sat_tr.notna()].sort_values('sat_tr'); ax[0].plot(g.sat_tr,100*np.arange(1,len(g)+1)/len(g),color='#b9770e',lw=1.2,label='Mareograflar (eşit ağırlık)')
ax[0].axvline(GM,color='#999',ls='--',lw=.8); ax[0].text(GM+.2,4,'Küresel ortalama',fontsize=7.5,color='#666'); ax[0].set_xlim(-6,14); ax[0].set_xlabel('Göreli deniz seviyesi eğilimi (mm/yıl)'); ax[0].set_ylabel('Birikimli pay (%)'); ax[0].legend(frameon=False,fontsize=7.5,loc='center right'); ax[0].set_title('(a) Birikimli dağılım',fontsize=9,loc='left')
steps=[('Küresel\nokyanus ort.',R['alt_mean']),('Limanlarda\naltimetre',R['hyb_est']['alt_all']['est']),('Mareograflı\nlimanlar:\naltimetre',R['hyb_est']['gauge_ports_alt']['est']),('Mareograflı\nlimanlar:\nmareograf',R['hyb_est']['gauge_ports_gauge']['est']),('Bileşik\ntahmin',R['hyb_est']['best_all']['est'])]
ci=[None,R['hyb_est']['alt_all'],R['hyb_est']['gauge_ports_alt'],R['hyb_est']['gauge_ports_gauge'],R['hyb_est']['best_all']]
x=np.arange(len(steps)); ax[1].bar(x,[s[1] for s in steps],color=['#999',C['alt'],C['alt'],C['gauge'],'#6c3483'],width=.62)
for i,c in enumerate(ci):
    if c: ax[1].plot([i,i],[c['lo'],c['hi']],color='k',lw=1)
    ax[1].text(i,steps[i][1]+(0.5 if c else .08),f"{steps[i][1]:.2f}".replace('.',','),ha='center',fontsize=8)
ax[1].set_xticks(x); ax[1].set_xticklabels([s[0] for s in steps],fontsize=7); ax[1].set_ylim(2.5,4.9); ax[1].set_ylabel('Ticaret ağırlıklı eğilim (mm/yıl)'); ax[1].set_title('(b) Ortalamanın katmanları',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('fig/Figure_4.png'); plt.close()
# S5 istasyon sayisi + senaryo
fig,ax=plt.subplots(1,2,figsize=(7.4,3.0))
sy=pd.Series({int(k):v for k,v in R['stations_by_year'].items()}).loc[1950:2024]; ax[0].plot(sy.index,sy.values,color=C['gauge'],lw=1.4); ax[0].set_ylabel('Veri gönderen mareograf sayısı'); ax[0].set_title('(a) PSMSL RLR istasyonları',fontsize=9,loc='left'); ax[0].set_ylim(0,None)
k=np.array([0,1,2,3,5]); ax[1].plot(k,[R['scen'][str(i)] for i in k],marker='o',color='#6c3483',lw=1.4); ax[1].axhline(GM,color='#999',ls='--',lw=.8); ax[1].text(0.1,GM+.07,'Küresel ortalama',fontsize=7.5,color='#666')
ax[1].set_xlabel('Ölçülmeyen limanlarda varsayılan ortalama çökme (mm/yıl)'); ax[1].set_ylabel('Ticaret ağırlıklı göreli eğilim (mm/yıl)'); ax[1].set_ylim(3,6.4); ax[1].set_title('(b) Bilinmeyen kara hareketine duyarlılık',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('fig/Figure_5.png'); plt.close()
print(sy.loc[[1970,1980,1990,2000,2010,2015,2020,2023,2024]].to_dict(), sy.idxmax(), sy.max())

# ------------------------------------------------------------------------------
# ## Şekiller / Figures
# ------------------------------------------------------------------------------
import glob
print('Sekiller / figures:', sorted(glob.glob('fig/*.png')))