# -*- coding: utf-8 -*- # ------------------------------------------------------------------------------ # 2020 küresel kükürt sınırı ve 2025 Akdeniz emisyon kontrol alanından sonra limanların rüzgâr altında kükürt dioksit # Sulphur dioxide downwind of ports after the 2020 global sulphur cap and the 2025 Mediterranean emission control area # # GOMD Yayınlar · Makale 30 — https://gomd.pirireis.edu.tr/publications.php?id=rp30 # # 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 pyfixest scipy shapely # Çalıştırma / Run: python rp30_port_sector_so2_mediterranean_seca.py # # Air quality data: European Environment Agency (EEA), Air Quality e-Reporting, hourly measurements (validated E1a and up-to-date E2a data flows), retrieved from the Air Quality Download Service on 14 November 2026; re-used under CC BY 4.0 with the EEA as the original source. Winds: NOAA National Centers for Environmental Information, Global Historical Climatology Network hourly (GHCNh). Port locations: IMF PortWatch / UN Global Platform (https://portwatch.imf.org). Coastline: Natural Earth. The package contains monthly summaries by sampling point derived from these sources, not the raw hourly data. # 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 = ['run_download.py', 'a6_final.py', 'a10_windreg.py', 'a11_fig.py'] SKIP_IF_EXISTS = {'run_download.py': 'eea/wind_month.parquet'} # hazırlık modülü: çıktısı pakette varsa atlanır / preparation module: skipped if its output is in the package DIRS = ['eea', 'ghcnh', 'data', 'data/pw', 'paper', 'paper/out', 'paper/fig'] # ============================================================================== # run_download.py # ============================================================================== MODULES['run_download.py'] = r''' # Ham veriyi kaynaklardan indirir ve ara tablolari uretir (yaklasik 4 GB indirme, 30-60 dakika). # Paketteki ara tablolar (eea/wind_month.parquet vb.) mevcutsa bu modul atlanir. import subprocess, sys for a in (['a0_meta.py'],['eea_dlx.py','SO2','2'],['eea_dlx.py','SO2','1'],['a1_stations.py'],['a2x_panel.py'],['a7_ghcnh.py'],['a8_hourly.py'],['a9_wind.py']): print('==',a,flush=True); subprocess.run([sys.executable]+a,check=True) ''' # ============================================================================== # a6_final.py # ============================================================================== MODULES['a6_final.py'] = r''' # Akdeniz SOx ECA (1 Mayis 2025) ve IMO 2020: liman istasyonlarinda SO2 -- fark-icinde-fark, olay calismasi, guc import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore') import pyfixest as pf st=pd.read_parquet('eea/stations.parquet'); s=st[st['Air Pollutant']=='SO2'].drop_duplicates('sp') m=pd.read_parquet('eea/panelx_month_SO2.parquet').merge(s,on='sp'); m['t']=m.ym.dt.to_timestamp(); m=m[m.t<'2026-11-01'] MEDC=['Spain','France','Italy','Greece','Croatia','Slovenia','Malta','Cyprus','Türkiye','Albania','Montenegro'] def group(r): if r.d_port<=5: return 'MedPort' if r.basin=='MED' else ('Port_'+r.basin) if r.d_coast>50: return 'Inland_Med' if r.country in MEDC else 'Inland_Oth' if r.basin=='MED' and r.d_port<=25 and r.d_coast<=25: return 'MedNear' return 'other' m['grp']=m.apply(group,axis=1); m['ymc']=m.ym.astype(str); m['cym']=m.country+'_'+m.ymc; m['lv']=np.log(m.v.clip(lower=0.2)) R={} def did(d,tr,ctrl,y='v',fe='sp + ymc',extra=''): d=d[d.grp.isin([tr]+ctrl)].copy(); d['D']=((d.grp==tr)&(d.post==1)).astype(int) mm=pf.feols('%s ~ D %s| %s'%(y,extra,fe),data=d,vcov={'CRV1':'sp'}); t=mm.tidy().loc['D']; pre=d[(d.grp==tr)&(d.post==0)][y].mean() return {'b':float(t['Estimate']),'se':float(t['Std. Error']),'p':float(t['Pr(>|t|)']),'n_sp':int(d.sp.nunique()),'n_tr':int(d[d.grp==tr].sp.nunique()),'n_obs':int(len(d)),'pre_mean':float(pre)} def window(a,c,e,minpre,minpost,src=None): d=(m if src is None else src); d=d[(d.t>=a)&(d.t=c).astype(int) npre=d[d.post==0].groupby('sp').size(); npost=d[d.post==1].groupby('sp').size(); keep=set(npre[npre>=minpre].index)&set(npost[npost>=minpost].index); return d[d.sp.isin(keep)] W={'seca':window('2021-01-01','2025-05-01','2026-11-01',30,9),'imo':window('2016-01-01','2020-01-01','2022-01-01',30,16)} for k,d in W.items(): u=d.drop_duplicates('sp'); R[k+'_sample']={'groups':u.grp.value_counts().to_dict(),'countries_medport':u[u.grp=='MedPort'].country.value_counts().to_dict(),'ports_medport':int(u[u.grp=='MedPort'].portname.nunique()),'months':[str(d.t.min())[:7],str(d.t.max())[:7]], 'desc':{g:{'pre':float(x[x.post==0].v.mean()),'post':float(x[x.post==1].v.mean()),'pre_max':float(x[x.post==0].vmax.mean()),'post_max':float(x[x.post==1].vmax.mean()),'n':int(x.sp.nunique()),'share_below2':float((x.v<2).mean())} for g,x in d.groupby('grp')}} res={} for y in ('v','vmax','hi','lv'): res[y]={'inland_cxm':did(d,'MedPort',['Inland_Med'],y,'sp + cym'),'inland':did(d,'MedPort',['Inland_Med'],y),'nseca_ports':did(d,'MedPort',['Port_NSECA'],y),'atl_ports':did(d,'MedPort',['Port_ATL'],y),'all_inland':did(d,'MedPort',['Inland_Med','Inland_Oth'],y)} R[k]=res; print(k,{y:{c:(round(v['b'],3),round(v['se'],3),round(v['p'],3)) for c,v in res[y].items()} for y in ('v','vmax','hi')},flush=True) # mesafe dd=d.copy(); dd['grp0']=dd.grp out={} for lab,lo,hi in (('0-2',0,2),('2-5',2,5),('5-10',5,10),('10-25',10,25)): x=dd.copy(); x['grp']=np.where((x.basin=='MED')&(x.d_port>lo)&(x.d_port<=hi)&(x.d_coast<=25),'T',np.where(x.grp0=='Inland_Med','Inland_Med','zz')) try: out[lab]=did(x,'T',['Inland_Med'],'v','sp + cym') except Exception as e: out[lab]=None R[k+'_distance']=out; print(k,'dist',{a:(round(b['b'],3),round(b['se'],3),b['n_tr']) for a,b in out.items() if b},flush=True) # heterojenlik het={}; med_calls=d[d.grp=='MedPort'].drop_duplicates('sp').port_calls.median() for lab,mask in (('large_port',(d.port_calls>=med_calls)),('small_port',(d.port_calls=2].index)],'MedPort',['Inland_Med'],'v','sp + cym'), 'background_only':did(d[d.stype=='background'],'MedPort',['Inland_Med'],'v','sp + cym'), 'within3km':did(d[(d.grp!='MedPort')|(d.d_port<=3)],'MedPort',['Inland_Med'],'v','sp + cym')} d=W['seca'] R['seca_rob']['verified_only_2025']=did(d[d.t<'2026-01-01'],'MedPort',['Inland_Med'],'v','sp + cym'); R['seca_rob']['verified_only_2025_nseca']=did(d[d.t<'2026-01-01'],'MedPort',['Port_NSECA'],'v') R['seca_rob']['same_months_May_Oct']=did(d[d.t.dt.month.isin([5,6,7,8,9,10])],'MedPort',['Inland_Med'],'v','sp + cym') R['seca_rob']['summer_JJA']=did(d[d.t.dt.month.isin([6,7,8])],'MedPort',['Inland_Med'],'v','sp + cym') R['seca_rob']['pre_from_2022']=did(d[d.t>='2022-01-01'],'MedPort',['Inland_Med'],'v','sp + cym') # plasebo tarihler (yalniz politika oncesi veri) pl={} for c in ('2023-05-01','2024-05-01'): x=window('2021-01-01',c,'2025-05-01',20,9); pl[c[:7]]=did(x,'MedPort',['Inland_Med'],'v','sp + cym') R['seca_placebo_dates']=pl; print('placebo',{k:(round(v['b'],3),round(v['se'],3)) for k,v in pl.items()}) # olay calismasi def es(d,pol,bins,ref,ctrl=['Inland_Med'],fe='sp + cym'): x=d[d.grp.isin(['MedPort']+ctrl)].copy(); x['k']=bins(x.t); out={} for kk in sorted(x.k.unique()): if kk==ref: continue x['D_%s'%str(kk).replace('-','m')]=((x.grp=='MedPort')&(x.k==kk)).astype(int) cols=[c for c in x.columns if c.startswith('D_')] mm=pf.feols('v ~ %s | %s'%(' + '.join(cols),fe),data=x,vcov={'CRV1':'sp'}); t=mm.tidy() return {c[2:].replace('m','-'):[float(t.loc[c,'Estimate']),float(t.loc[c,'Std. Error'])] for c in cols} half=lambda t: ((t.dt.year-2025)*12+t.dt.month-5)//6 # 6 aylik donemler, 0 = Mayis-Ekim 2025 R['seca_event']=es(W['seca'],'seca',half,-1); R['seca_event_nseca']=es(W['seca'],'seca',half,-1,['Port_NSECA'],'sp + ymc') yr=lambda t: t.dt.year-2020 R['imo_event']=es(W['imo'],'imo',yr,-1); R['imo_event_nseca']=es(W['imo'],'imo',yr,-1,['Port_NSECA'],'sp + ymc') print('seca_event',{k:round(v[0],2) for k,v in R['seca_event'].items()}); print('imo_event',{k:round(v[0],2) for k,v in R['imo_event'].items()}) # ATL limanlari IMO 2020 (ayni ulke ici: Ispanya/Portekiz/Fransa ic istasyonlari) x=W['imo'].copy(); x['grp']=np.where(x.grp=='Port_ATL','T',np.where((x.d_coast>50)&x.country.isin(['Spain','Portugal','France','Ireland']),'C','zz')) R['imo_atl_vs_inland']=did(x,'T',['C'],'v','sp + cym'); print('ATL imo',R['imo_atl_vs_inland']) x=W['seca'].copy(); x['grp']=np.where(x.grp=='Port_ATL','T',np.where((x.d_coast>50)&x.country.isin(['Spain','Portugal','France','Ireland']),'C','zz')) R['seca_atl_vs_inland']=did(x,'T',['C'],'v','sp + cym') x=W['seca'].copy(); x['grp']=np.where(x.grp=='Port_NSECA','T',np.where(x.grp=='Inland_Oth','C','zz')); R['seca_nseca_vs_inland']=did(x,'T',['C'],'v','sp + cym') # guc b=R['seca']['v']['inland_cxm']; R['power']={'mde80':2.8*b['se'],'ci90':[b['b']-1.645*b['se'],b['b']+1.645*b['se']],'ci95':[b['b']-1.96*b['se'],b['b']+1.96*b['se']],'pre_mean':b['pre_mean'], 'imo_effect':R['imo']['v']['inland_cxm']['b'],'expected_seca':0.2*R['imo']['v']['inland_cxm']['b'],'expected_seca_nseca':0.2*R['imo']['v']['nseca_ports']['b']} print(R['power']) json.dump(R,open('eea/res_so2.json','w'),indent=1,default=str) # sekil verisi: grup aylik ortalamalari g=m[m.sp.isin(set(W['seca'].sp)|set(W['imo'].sp))] both=set(W['seca'].sp)&set(W['imo'].sp) ts=m[m.sp.isin(both)&m.grp.isin(['MedPort','Inland_Med','Port_NSECA','Port_ATL'])].groupby(['grp','t']).v.mean().reset_index(); ts.to_csv('eea/ts_groups.csv',index=False); print('both',len(both)) ''' # ============================================================================== # a10_windreg.py # ============================================================================== MODULES['a10_windreg.py'] = r''' import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore') import pyfixest as pf A=pd.read_parquet('eea/wind_month.parquet'); A=A[(A.size_port>=24)&(A.size_away>=24)].copy() A['t']=A.ym.dt.to_timestamp(); A=A[A.t<'2026-11-01']; A['ymc']=A.ym.astype(str); A['moy']=A.t.dt.month.astype(str); A['yr']=A.t.dt.year A['lratio']=np.log(A.mean_port.clip(lower=0.2)/A.mean_away.clip(lower=0.2)) R={'n_sp':int(A.sp.nunique()),'basins':A.drop_duplicates('sp').basin.value_counts().to_dict(),'met_dist_median':float(A.drop_duplicates('sp').met_d0.median()), 'hours_port_share':float(A.size_port.sum()/(A.size_port.sum()+A.size_away.sum()+A.get('size_side',0).sum()))} R['annual']={b:{int(y):{'port':float(g.mean_port.mean()),'away':float(g.mean_away.mean()),'excess':float(g.excess.mean()),'n':int(g.sp.nunique())} for y,g in x.groupby('yr')} for b,x in A.groupby('basin')} def win(a,c,e,minpre,minpost): d=A[(A.t>=a)&(A.t=c).astype(int); npre=d[d.post==0].groupby('sp').size(); npost=d[d.post==1].groupby('sp').size() return d[d.sp.isin(set(npre[npre>=minpre].index)&set(npost[npost>=minpost].index))] def did(d,tr,ctrl,y='excess',fe='sp + ymc'): d=d[d.basin.isin(tr+ctrl)].copy(); d['D']=(d.basin.isin(tr)&(d.post==1)).astype(int) m=pf.feols('%s ~ D | %s'%(y,fe),data=d,vcov={'CRV1':'sp'}); t=m.tidy().loc['D'] return {'b':float(t['Estimate']),'se':float(t['Std. Error']),'p':float(t['Pr(>|t|)']),'n_tr':int(d[d.basin.isin(tr)].sp.nunique()),'n_ct':int(d[d.basin.isin(ctrl)].sp.nunique()),'n_obs':int(len(d)),'pre_tr':float(d[d.basin.isin(tr)&(d.post==0)][y].mean()),'pre_ct':float(d[d.basin.isin(ctrl)&(d.post==0)][y].mean())} def prepost(d,b,y='excess'): d=d[d.basin==b]; m=pf.feols('%s ~ post | sp + moy'%y,data=d,vcov={'CRV1':'sp'}); t=m.tidy().loc['post'] return {'b':float(t['Estimate']),'se':float(t['Std. Error']),'p':float(t['Pr(>|t|)']),'n':int(d.sp.nunique()),'pre':float(d[d.post==0][y].mean())} W={'imo':win('2016-01-01','2020-01-01','2022-01-01',30,16),'seca':win('2021-01-01','2025-05-01','2026-11-01',30,9)} for k,d in W.items(): R[k]={'n':d.drop_duplicates('sp').basin.value_counts().to_dict(),'med_countries':d[d.basin=='MED'].drop_duplicates('sp').country.value_counts().to_dict(),'med_ports':int(d[d.basin=='MED'].portname.nunique())} for y in ('excess','lratio','mean_port','mean_away'): R[k][y]={'med_vs_nseca':did(d,['MED'],['NSECA'],y),'prepost':{b:prepost(d,b,y) for b in ('MED','NSECA','ATL','BLK') if (d.basin==b).any()}} if k=='imo': R[k]['excess']['medatl_vs_nseca']=did(d,['MED','ATL'],['NSECA']); R[k]['excess']['atl_vs_nseca']=did(d,['ATL'],['NSECA']) else: R[k]['excess']['med_vs_atl']=did(d,['MED'],['ATL']); R[k]['excess']['med_vs_nseca_atl']=did(d,['MED'],['NSECA','ATL']); R[k]['excess']['atl_vs_nseca']=did(d,['ATL'],['NSECA']) # heterojenlik ve saglamlik (MED vs NSECA) h={} medc=d[d.basin=='MED'].drop_duplicates('sp').port_calls.median() for lab,mask in (('within2km',d.d_port<=2),('2to5km',d.d_port>2),('large_port',d.port_calls>=medc),('small_port',d.port_calls='2022-01-01'],['MED'],['NSECA']),'summer':did(d[d.t.dt.month.isin([5,6,7,8,9,10])],['MED'],['NSECA'])} pl={} for c in ('2022-05-01','2023-05-01','2024-05-01'): x=A[(A.t>='2021-01-01')&(A.t<'2025-05-01')].copy(); x['post']=(x.t>=c).astype(int); x=x[x.sp.isin(set(d.sp))]; pl[c[:7]]=did(x,['MED'],['NSECA']) R['seca']['placebo']=pl; print('placebo',{k:(round(v['b'],3),round(v['se'],3)) for k,v in pl.items()}); print('rob',{k:(round(v['b'],3),round(v['se'],3)) for k,v in R['seca']['rob'].items()}) # olay calismasi: yillik (MED vs NSECA), ref 2019; ve yari yillik SECA def es(d,key,ref,tr=['MED'],ctrl=['NSECA']): x=d[d.basin.isin(tr+ctrl)].copy(); x['k']=key(x.t); cols=[] for kk in sorted(x.k.unique()): if kk==ref: continue c='D_%s'%str(kk).replace('-','m'); x[c]=(x.basin.isin(tr)&(x.k==kk)).astype(int); cols.append(c) m=pf.feols('excess ~ %s | sp + ymc'%' + '.join(cols),data=x,vcov={'CRV1':'sp'}); t=m.tidy() return {c[2:].replace('m','-'):[float(t.loc[c,'Estimate']),float(t.loc[c,'Std. Error'])] for c in cols} full=A[A.sp.isin(set(W['imo'].sp)&set(W['seca'].sp))]; R['n_full']=full.drop_duplicates('sp').basin.value_counts().to_dict() R['event_year']=es(full[full.t>='2016-01-01'],lambda t:t.dt.year,2019); R['event_year_atl']=es(full,lambda t:t.dt.year,2019,['ATL'],['NSECA']) R['event_half']=es(W['seca'],lambda t:((t.dt.year-2025)*12+t.dt.month-5)//6,-1) print('event_year',{k:(round(v[0],2),round(v[1],2)) for k,v in R['event_year'].items()}); print('event_half',{k:(round(v[0],2),round(v[1],2)) for k,v in R['event_half'].items()}) b=R['seca']['excess']['med_vs_nseca']; i=R['imo']['excess']['med_vs_nseca'] R['power']={'mde80':2.8*b['se'],'ci95':[b['b']-1.96*b['se'],b['b']+1.96*b['se']],'ci90':[b['b']-1.645*b['se'],b['b']+1.645*b['se']],'pre_excess_med':b['pre_tr'],'imo_drop':i['b'],'imo_pre':i['pre_tr']} print(R['power']) json.dump(R,open('eea/res_wind.json','w'),indent=1) ''' # ============================================================================== # a11_fig.py # ============================================================================== MODULES['a11_fig.py'] = r''' from figstyle import * import json, numpy as np, pandas as pd, glob, os W=json.load(open('eea/res_wind.json')); M=json.load(open('eea/res_so2.json')) A=pd.read_parquet('eea/wind_month.parquet'); s=pd.read_parquet('ghcnh/port_stations.parquet') if os.path.isdir('eea/dailyx_SO2_2') and os.path.isfile('eea/panelx_daily_SO2.parquet'): info={'n_points_all':int(pd.read_parquet('eea/panelx_month_SO2.parquet').sp.nunique()),'n_daily':int(len(pd.read_parquet('eea/panelx_daily_SO2.parquet',columns=['sp']))), 'files_verified':len(glob.glob('eea/dailyx_SO2_2/*')),'files_utd':len(glob.glob('eea/dailyx_SO2_1/*')),'port_points':int(len(s)),'port_points_met40':int((s.met_d0<=40).sum()), 'met_stations':int(s[s.met_d0<=40].met0.nunique()),'hours_matched':int(A[['size_port','size_away','size_side']].sum().sum()),'hours_port':int(A.size_port.sum()),'hours_away':int(A.size_away.sum()), 'wind_points':int(A.sp.nunique()),'wind_ports':int(A.portname.nunique()),'wind_countries':int(A.country.nunique()),'wind_basin':A.drop_duplicates('sp').basin.value_counts().to_dict(), 'med_countries':A[A.basin=='MED'].drop_duplicates('sp').country.value_counts().to_dict(),'met_d_median':float(A.drop_duplicates('sp').met_d0.median()),'d_port_median':float(A.drop_duplicates('sp').d_port.median()), 'stype':A.drop_duplicates('sp').stype.value_counts().to_dict()} json.dump({k:info[k] for k in ('n_points_all','n_daily','files_verified','files_utd','port_points','port_points_met40','met_stations')},open('eea/info_static.json','w')) else: info=json.load(open('eea/info_static.json')) # ham indirme klasorleri yoksa sayilar paketteki dosyadan okunur s=s[s.met_d0<=40] info.update({'hours_matched':int(A[['size_port','size_away','size_side']].sum().sum()),'hours_port':int(A.size_port.sum()),'hours_away':int(A.size_away.sum()),'wind_points':int(A.sp.nunique()),'wind_ports':int(A.portname.nunique()),'wind_countries':int(A.country.nunique()),'wind_basin':A.drop_duplicates('sp').basin.value_counts().to_dict(),'med_countries':A[A.basin=='MED'].drop_duplicates('sp').country.value_counts().to_dict(),'met_d_median':float(A.drop_duplicates('sp').met_d0.median()),'d_port_median':float(A.drop_duplicates('sp').d_port.median()),'stype':A.drop_duplicates('sp').stype.value_counts().to_dict()}) try: info['no2']=json.load(open('eea/res_no2.json')) except Exception: pass json.dump({'wind':W,'mean':M,'info':info},open('paper/out/results.json','w'),indent=1); print(info) col={'MED':C['r'],'NSECA':C['b'],'ATL':C['g'],'BLK':C['gr']}; nm={'MED':'Mediterranean','NSECA':'North Sea and Baltic','ATL':'Atlantic','BLK':'Black Sea and Marmara'} fig,ax=plt.subplots(1,2,figsize=(7.4,3.1)) for b in ('MED','ATL','NSECA'): y=sorted(int(k) for k in W['annual'][b]); ax[0].plot(y,[W['annual'][b][str(k)]['excess'] for k in y],color=col[b],marker='o',ms=3,lw=1.4,label=nm[b]) for a in ax: a.axvline(2019.5,color='k',lw=.7,ls=':'); a.axvline(2025.33,color='k',lw=.7,ls='--') ax[0].set_ylabel('Port-sector excess SO$_2$ (µg m$^{-3}$)'); ax[0].set_title('(a) Port-sector excess by sea region'); ax[0].legend(fontsize=7); ax[0].set_ylim(0,2.6) y=sorted(int(k) for k in W['annual']['MED']); ax[1].plot(y,[W['annual']['MED'][str(k)]['port'] for k in y],color=C['r'],marker='o',ms=3,lw=1.4,label='Wind from port sector'); ax[1].plot(y,[W['annual']['MED'][str(k)]['away'] for k in y],color=C['k'],marker='s',ms=3,lw=1.2,label='Wind from opposite sectors') ax[1].set_ylabel('Mean SO$_2$ (µg m$^{-3}$)'); ax[1].set_title('(b) Mediterranean port monitors'); ax[1].legend(fontsize=7); ax[1].set_ylim(0,7) fig.tight_layout(); fig.savefig('paper/fig/Figure_1.png'); plt.close() fig,ax=plt.subplots(1,2,figsize=(7.4,3.0)) for a,(key,tt,c) in zip(ax,(('event_year','(a) Mediterranean relative to North Sea/Baltic',C['r']),('event_year_atl','(b) Atlantic relative to North Sea/Baltic',C['g']))): E=dict(W[key]); E['2019']=[0,0]; ks=sorted(E,key=int); e=np.array([E[k][0] for k in ks]); se=np.array([E[k][1] for k in ks]); x=[int(k) for k in ks] a.errorbar(x,e,yerr=1.96*se,fmt='o-',color=c,ms=3.5,lw=1.2,capsize=2,elinewidth=.7); a.axhline(0,color='k',lw=.6); a.axvline(2019.5,color='k',lw=.7,ls=':'); a.axvline(2025.33,color='k',lw=.7,ls='--'); a.set_title(tt,fontsize=8.5); a.set_ylim(-2.2,1.6) ax[0].set_ylabel('Difference in port-sector excess\nrelative to 2019 (µg m$^{-3}$)') fig.tight_layout(); fig.savefig('paper/fig/Figure_2.png'); plt.close() fig,ax=plt.subplots(1,2,figsize=(7.4,3.0)) E=dict(W['event_half']); E['-1']=[0,0]; ks=sorted([k for k in E if int(k)>=-8],key=int); e=np.array([E[k][0] for k in ks]); se=np.array([E[k][1] for k in ks]); x=[int(k) for k in ks] ax[0].errorbar(x,e,yerr=1.96*se,fmt='o-',color=C['r'],ms=3.5,lw=1.2,capsize=2,elinewidth=.7); ax[0].axhline(0,color='k',lw=.6); ax[0].axvline(-0.5,color='k',lw=.7,ls='--'); ax[0].set_xlabel('Half-years from 1 May 2025'); ax[0].set_ylabel('Difference in port-sector excess (µg m$^{-3}$)'); ax[0].set_title('(a) Around the Mediterranean SECA') rows=[('Main estimate',W['seca']['excess']['med_vs_nseca']),('Control: Atlantic ports',W['seca']['excess']['med_vs_atl']),('Control: both',W['seca']['excess']['med_vs_nseca_atl']),('Verified data only (2025)',W['seca']['rob']['verified_2025_only']),('Pre-period from 2022',W['seca']['rob']['pre_from_2022']),('May–October only',W['seca']['rob']['summer']),('Placebo: May 2022',W['seca']['placebo']['2022-05']),('Placebo: May 2023',W['seca']['placebo']['2023-05']),('Placebo: May 2024',W['seca']['placebo']['2024-05'])] yy=np.arange(len(rows))[::-1]; ax[1].errorbar([r[1]['b'] for r in rows],yy,xerr=[1.96*r[1]['se'] for r in rows],fmt='o',color=C['b'],ms=3.5,capsize=2,elinewidth=.7); ax[1].axvline(0,color='k',lw=.6); ax[1].set_yticks(yy); ax[1].set_yticklabels([r[0] for r in rows],fontsize=7); ax[1].set_xlabel('Estimate (µg m$^{-3}$)'); ax[1].set_title('(b) Alternative controls and placebo dates') fig.tight_layout(); fig.savefig('paper/fig/Figure_3.png'); plt.close() fig,ax=plt.subplots(1,2,figsize=(7.4,3.1),sharey=True) labs=[('Spain','Spain'),('Italy','Italy'),('other_med','Other Mediterranean'),('large_port','Larger ports'),('small_port','Smaller ports'),('within2km','Within 2 km of port'),('2to5km','2–5 km from port'),('industrial','Industrial sites'),('non_industrial','Other sites')] for a,(k,tt) in zip(ax,(('imo','(a) 2020 global sulphur cap'),('seca','(b) 2025 Mediterranean SECA'))): H=W[k]['het']; yy=np.arange(len(labs))[::-1]; a.errorbar([H[x]['b'] for x,_ in labs],yy,xerr=[1.96*H[x]['se'] for x,_ in labs],fmt='o',color=C['r'],ms=3.5,capsize=2,elinewidth=.7); a.axvline(0,color='k',lw=.6); a.set_title(tt,fontsize=8.5); a.set_xlabel('Change in port-sector excess (µg m$^{-3}$)'); a.set_xlim(-2.6,1.6) ax[0].set_yticks(yy); ax[0].set_yticklabels([l for _,l in labs],fontsize=7.5) fig.tight_layout(); fig.savefig('paper/fig/Figure_4.png'); plt.close(); print('ok') ''' # ============================================================================== # 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'} ''' # ============================================================================== # a0_meta.py # ============================================================================== MODULES['a0_meta.py'] = r''' # EEA istasyon metaverisi ve Natural Earth kiyi cizgisini indirir (eea/meta.parquet, ne/ne_10m_coastline.shp) import urllib.request, zipfile, io, os, pandas as pd os.makedirs('eea',exist_ok=True); os.makedirs('ne',exist_ok=True) UA={'User-Agent':'gomd-research/1.0'} b=urllib.request.urlopen(urllib.request.Request('https://discomap.eea.europa.eu/App/AQViewer/download?fqn=Airquality_Dissem.b2g.measurements&f=csv',headers=UA),timeout=600).read() z=zipfile.ZipFile(io.BytesIO(b)); name=[n for n in z.namelist() if n.lower().endswith('.csv')][0] m=pd.read_csv(z.open(name),low_memory=False); m.to_parquet('eea/meta.parquet'); print('metaveri satiri',len(m)) b=urllib.request.urlopen(urllib.request.Request('https://naciscdn.org/naturalearth/10m/physical/ne_10m_coastline.zip',headers=UA),timeout=600).read() zipfile.ZipFile(io.BytesIO(b)).extractall('ne'); print('kiyi cizgisi indirildi') ''' # ============================================================================== # eea_dlx.py # ============================================================================== MODULES['eea_dlx.py'] = r''' # EEA hava kalitesi saatlik verisi (parquet) -> ornekleme noktasi x gun ortalamasi. Kullanim: eea_dl.py POLLUTANT DATASET(1=guncel dogrulanmamis, 2=dogrulanmis) import sys, os, io, json, urllib.request, concurrent.futures as cf, time, tempfile import pyarrow.parquet as pq, pyarrow as pa, pandas as pd, numpy as np POL, DS = sys.argv[1], int(sys.argv[2]); A = 'https://eeadmz1-downloads-api-appservice.azurewebsites.net/'; START = pd.Timestamp('2016-01-01') OUT = 'eea/dailyx_%s_%d' % (POL, DS); os.makedirs(OUT, exist_ok=True) def post(path, body, timeout=300): r = urllib.request.Request(A + path, data=json.dumps(body).encode(), headers={'Content-Type': 'application/json', 'User-Agent': 'gomd-research/1.0'}) return urllib.request.urlopen(r, timeout=timeout).read().decode('utf-8-sig') lf = 'eea/urls_%s_%d.txt' % (POL, DS) if not os.path.isfile(lf): t = post('ParquetFile/urls', dict(countries=[], cities=[], pollutants=[POL], dataset=DS, source='Api')); open(lf, 'w').write(t) urls = [l.strip() for l in open(lf) if l.strip().startswith('http')] def one(u): name = u.split('/')[-2] + '_' + u.split('/')[-1]; f = os.path.join(OUT, name) if os.path.isfile(f): return 0 for k in range(4): try: b = urllib.request.urlopen(urllib.request.Request(u, headers={'User-Agent': 'gomd-research/1.0'}), timeout=300).read() t = pq.read_table(io.BytesIO(b), columns=['Samplingpoint', 'Start', 'Value', 'AggType', 'Validity']).to_pandas() t = t[(t.Start >= START) & (t.Validity > 0)] if len(t) == 0: pd.DataFrame(columns=['sp', 'date', 'value', 'n', 'agg']).to_parquet(f); return 1 t['Value'] = t.Value.astype(float); t['date'] = t.Start.dt.floor('D') g = t.groupby(['Samplingpoint', 'date', 'AggType']).Value.agg(['mean', 'count', 'max']).reset_index(); g.columns = ['sp', 'date', 'agg', 'value', 'n', 'vmax'] g.to_parquet(f + '.part'); os.replace(f + '.part', f); return 1 except Exception as e: err = str(e)[:120]; time.sleep(2 * (k + 1)) open(OUT + '.err', 'a').write(u + ' ' + err + '\n'); return -1 with cf.ThreadPoolExecutor(10) as ex: res = list(ex.map(one, urls)) print(POL, DS, 'dosya', len(urls), 'indirilen', res.count(1), 'vardi', res.count(0), 'hata', res.count(-1), flush=True) ''' # ============================================================================== # a1_stations.py # ============================================================================== MODULES['a1_stations.py'] = r''' import pandas as pd, numpy as np, struct, sys from shapely.geometry import LineString, Point from shapely.strtree import STRtree m=pd.read_parquet('eea/meta.parquet') m=m[m['Air Pollutant'].isin(['SO2','NO2'])] m['sp']=m['Sampling Point Id'] st=m.sort_values('Year').groupby(['sp','Air Pollutant']).agg(country=('Country','last'),eoi=('Air Quality Station EoI Code','last'),name=('Air Quality Station Name','last'),lon=('Longitude','last'),lat=('Latitude','last'),area=('Air Quality Station Area','last'),stype=('Air Quality Station Type','last'),alt=('Altitude','last'),mtype=('Measurement Type','last')).reset_index() st=st[st.lon.notna()&st.lat.notna()] # kiyi cizgisi (shapefile okuyucu: pyshp yoksa elle) def read_shp(fn): b=open(fn,'rb').read(); pos=100; out=[] while posii',b[pos:pos+8]); rec=b[pos+8:pos+8+2*ln]; pos+=8+2*ln t=struct.unpack('=2: out.append(pts[a:bb]) return out lines=read_shp('ne/ne_10m_coastline.shp') # Avrupa-Akdeniz kutusu segs=[] for l in lines: if l[:,0].max()<-35 or l[:,0].min()>50 or l[:,1].max()<25 or l[:,1].min()>75: continue segs.append(l) def tomerc(lon,lat,lat0): return np.column_stack([lon*111.32*np.cos(np.radians(lat0)),lat*110.57]) # istasyon basina mesafe: yerel esit-alan yaklasimi (nokta bazli) from scipy.spatial import cKDTree pts=np.vstack([np.column_stack([np.interp(np.linspace(0,len(l)-1,max(2,int(len(l)*4))),np.arange(len(l)),l[:,0]),np.interp(np.linspace(0,len(l)-1,max(2,int(len(l)*4))),np.arange(len(l)),l[:,1])]) for l in segs]) def xyz(lon,lat): lo,la=np.radians(lon),np.radians(lat); return np.column_stack([np.cos(la)*np.cos(lo),np.cos(la)*np.sin(lo),np.sin(la)])*6371.0 T=cKDTree(xyz(pts[:,0],pts[:,1])) d,_=T.query(xyz(st.lon.values,st.lat.values)); st['d_coast']=d p=pd.read_csv('data/pw/ports_global.csv') TP=cKDTree(xyz(p.lon.values,p.lat.values)); k=min(5,len(p)) d,ix=TP.query(xyz(st.lon.values,st.lat.values),k=1) st['d_port']=d; st['portid']=p.portid.values[ix]; st['portname']=p.portname.values[ix]; st['port_iso']=p.ISO3.values[ix]; st['port_calls']=p.vessel_count_total.values[ix]; st['port_lon']=p.lon.values[ix]; st['port_lat']=p.lat.values[ix] # liman havzasi def basin(lon,lat,iso): # Akdeniz ECA: Cebelitarik (5.6W) dogusu, Canakkale'ye kadar; Karadeniz/Marmara haric if iso in ('ESP',): if lat<44 and lon>-5.6 and lat<42.6 and not (lat>42.9): if lon>-5.6 and (lat<43.0) and not (lon<-1.5 and lat>42): return 'MED' if not (lat>36.6 and lon<-5.3) else 'ATL' return 'ATL' if iso in ('FRA','MCO'): return 'MED' if (lat<44.2 and lon>2.5) else ('NSECA' if (lat>48.3 and lon>-5) else 'ATL') if iso in ('ITA','MLT','SVN','HRV','MNE','ALB','CYP','GIB','BIH'): return 'MED' if iso=='GRC': return 'MED' if iso=='TUR': return 'MED' if (lat<40.0 and lon<36.5) or (lat<38.5) else 'BLK' if iso in ('PRT','IRL','ISL'): return 'ATL' if iso=='GBR': return 'NSECA' if (lon>-5 and lat<58.5 and not (lon<-2.5 and lat>51.5 and lat<55.5)) else 'ATL' if iso in ('BEL','NLD','DEU','DNK','SWE','FIN','EST','LVA','LTU','POL'): return 'NSECA' if iso=='NOR': return 'NSECA' if lat<62 else 'ATL' if iso in ('BGR','ROU','UKR','GEO'): return 'BLK' return 'OTH' p['basin']=[basin(a,b,c) for a,b,c in zip(p.lon,p.lat,p.ISO3)] st=st.merge(p[['portid','basin']],on='portid',how='left') p.to_csv('eea/ports_basin.csv',index=False) st.to_parquet('eea/stations.parquet') print(len(st),st['Air Pollutant'].value_counts().to_dict()) s=st[st['Air Pollutant']=='SO2'] s['near']=pd.cut(s.d_port,[0,5,10,25,50,1e5],labels=['<5','5-10','10-25','25-50','>50']) print(pd.crosstab(s.basin,s.near)); print(pd.crosstab(s.country,s.basin).to_string()) ''' # ============================================================================== # a2x_panel.py # ============================================================================== MODULES['a2x_panel.py'] = r''' import pandas as pd, numpy as np, glob parts=[] for ds in (2,1): for f in glob.glob('eea/dailyx_SO2_%d/*'%ds): try: x=pd.read_parquet(f) except Exception: print('bad',f); continue if len(x)==0: continue x=x[x['agg']=='hour']; x['ds']=ds; parts.append(x) d=pd.concat(parts,ignore_index=True); d['sp']=d.sp.str.split('/').str[-1] d=d[(d.n>=18)&(d.value>-5)&(d.value<2000)] d=d.sort_values('ds').drop_duplicates(['sp','date'],keep='first') d.to_parquet('eea/panelx_daily_SO2.parquet') d['ym']=d.date.dt.to_period('M'); d['hi']=(d.vmax>=20).astype(float) m=d.groupby(['sp','ym']).agg(v=('value','mean'),vmax=('vmax','mean'),p95max=('vmax',lambda s:s.quantile(.95)),hi=('hi','mean'),nd=('value','size'),ds=('ds','max')).reset_index() m['dim']=m.ym.dt.days_in_month; m=m[m.nd>=0.75*m.dim]; m.to_parquet('eea/panelx_month_SO2.parquet') print('nokta',d.sp.nunique(),'gun',len(d),'ay',len(m)); print(m.groupby(m.ym.dt.year).agg(n=('sp','nunique'),v=('v','mean'),vmax=('vmax','mean'),hi=('hi','mean'))) ''' # ============================================================================== # a7_ghcnh.py # ============================================================================== MODULES['a7_ghcnh.py'] = r''' import pandas as pd, numpy as np, os, re, urllib.request, urllib.parse, concurrent.futures as cf, time from scipy.spatial import cKDTree B='https://www.ncei.noaa.gov/oa/global-historical-climatology-network/' os.makedirs('ghcnh/raw',exist_ok=True) def listkeys(prefix): keys=[]; marker='' while True: x=urllib.request.urlopen(B+'?prefix=%s&max-keys=1000&marker=%s'%(urllib.parse.quote(prefix),urllib.parse.quote(marker)),timeout=120).read().decode() k=re.findall(r'([^<]+)[^<]+[^<]+(\d+)',x); keys+=k if 'true' not in x or not k: break marker=k[-1][0] return keys if not os.path.isfile('ghcnh/keys2026.csv'): k=listkeys('hourly/access/by-year/2026/parquet/'); pd.DataFrame(k,columns=['key','size']).to_csv('ghcnh/keys2026.csv',index=False) k=pd.read_csv('ghcnh/keys2026.csv'); k['id']=k.key.str.extract(r'GHCNh_(.+)_2026')[0]; k=k[k['size']>150000] if not os.path.isfile('ghcnh/list.csv'): open('ghcnh/list.csv','wb').write(urllib.request.urlopen(B+'hourly/doc/ghcnh-station-list.csv',timeout=300).read()) L=pd.read_csv('ghcnh/list.csv',dtype=str); L['lat']=L.LATITUDE.astype(float); L['lon']=L.LONGITUDE.astype(float); L=L[L.GHCN_ID.isin(set(k.id))&L.lat.between(25,72)&L.lon.between(-35,50)] print('2026 dosyasi olan istasyon',len(k),'bolgede',len(L)) st=pd.read_parquet('eea/stations.parquet'); s=st[st['Air Pollutant']=='SO2'].drop_duplicates('sp') s=s[s.sp.isin(set(pd.read_parquet('eea/panelx_month_SO2.parquet').sp))&(s.d_port<=5)&s.basin.isin(['MED','NSECA','ATL','BLK'])].copy() dl=np.radians(s.port_lon-s.lon); la1=np.radians(s.lat); la2=np.radians(s.port_lat) s['bearing']=(np.degrees(np.arctan2(np.sin(dl)*np.cos(la2),np.cos(la1)*np.sin(la2)-np.sin(la1)*np.cos(la2)*np.cos(dl)))+360)%360 def xyz(lon,lat): lo,la=np.radians(lon),np.radians(lat); return np.column_stack([np.cos(la)*np.cos(lo),np.cos(la)*np.sin(lo),np.sin(la)])*6371.0 T=cKDTree(xyz(L.lon.values,L.lat.values)); d,ix=T.query(xyz(s.lon.values,s.lat.values),k=2) s['met0']=L.GHCN_ID.values[ix[:,0]]; s['met_d0']=d[:,0]; s['met1']=L.GHCN_ID.values[ix[:,1]]; s['met_d1']=d[:,1] s.to_parquet('ghcnh/port_stations.parquet'); print(s.met_d0.describe(),(s.met_d0<=30).sum(),len(s)) need=sorted(set(s[s.met_d0<=40].met0)); print('met istasyon',len(need)) def one(a): i,y=a; f='ghcnh/raw/%s_%d.parquet'%(i,y) if os.path.isfile(f): return 0 for kk in range(3): try: b=urllib.request.urlopen(B+'hourly/access/by-year/%d/parquet/GHCNh_%s_%d.parquet'%(y,i,y),timeout=120).read(); open(f+'.p','wb').write(b); os.replace(f+'.p',f); return 1 except Exception as e: if '404' in str(e): return -1 time.sleep(2) return -2 with cf.ThreadPoolExecutor(8) as ex: r=list(ex.map(one,[(i,y) for i in need for y in range(2016,2027)])) print('indirilen',r.count(1),'yok',r.count(-1),'hata',r.count(-2)) ''' # ============================================================================== # a8_hourly.py # ============================================================================== MODULES['a8_hourly.py'] = r''' import pandas as pd, numpy as np, io, os, urllib.request, concurrent.futures as cf, time, pyarrow.parquet as pq s=pd.read_parquet('ghcnh/port_stations.parquet'); s=s[s.met_d0<=40]; ids=set(s.sp); ids2={i.replace(':','_') for i in ids} urls=[] for ds in (2,1): for l in open('eea/urls_SO2_%d.txt'%ds): l=l.strip() if l.startswith('http') and l.split('/')[-1].replace('.parquet','') in ids2: urls.append((ds,l)) print('dosya',len(urls),'nokta',len(ids)) def one(a): ds,u=a; f='eea/hourly/%d_%s'%(ds,u.split('/')[-1]) if os.path.isfile(f): return 0 for k in range(4): try: b=urllib.request.urlopen(urllib.request.Request(u,headers={'User-Agent':'gomd-research/1.0'}),timeout=300).read() t=pq.read_table(io.BytesIO(b),columns=['Samplingpoint','Start','Value','AggType','Validity']).to_pandas() t=t[(t.Start>='2016-01-01')&(t.Validity>0)&(t.AggType=='hour')]; t['sp']=t.Samplingpoint.str.split('/').str[-1]; t['v']=t.Value.astype('float32') t[['sp','Start','v']].to_parquet(f+'.p'); os.replace(f+'.p',f); return 1 except Exception as e: err=str(e); time.sleep(2*(k+1)) print('hata',u,err[:80]); return -1 with cf.ThreadPoolExecutor(8) as ex: r=list(ex.map(one,urls)) print('indirilen',r.count(1),'hata',r.count(-1)) ''' # ============================================================================== # a9_wind.py # ============================================================================== MODULES['a9_wind.py'] = r''' # Ruzgar sektoru analizi: liman yonunden esen saatlerde SO2 fazlasi, politika oncesi/sonrasi import pandas as pd, numpy as np, glob, json, os, pyarrow.parquet as pq s=pd.read_parquet('ghcnh/port_stations.parquet'); s=s[s.met_d0<=40].set_index('sp') met={} def loadmet(i): if i in met: return met[i] parts=[] for f in sorted(glob.glob('ghcnh/raw/%s_*.parquet'%i)): try: t=pq.read_table(f,columns=['DATE','wind_direction','wind_speed']).to_pandas() except Exception as e: continue parts.append(t) if not parts: met[i]=None; return None t=pd.concat(parts); t['time']=pd.to_datetime(t.DATE,errors='coerce').dt.floor('h'); t['wd']=pd.to_numeric(t.wind_direction,errors='coerce'); t['ws']=pd.to_numeric(t.wind_speed,errors='coerce') t=t[t.time.notna()&t.wd.between(1,360)&t.ws.between(0.5,60)].drop_duplicates('time').set_index('time')[['wd','ws']].sort_index() met[i]=t; return t rows=[] for sp,g in pd.concat([pd.read_parquet(f).assign(ds=int(os.path.basename(f)[0])) for f in glob.glob('eea/hourly/*.parquet')]).groupby('sp'): if sp not in s.index: continue r=s.loc[sp]; w=loadmet(r.met0) if w is None or len(w)<20000: continue g=g[(g.v>-5)&(g.v<2000)].sort_values('ds').drop_duplicates('Start'); g['time']=pd.to_datetime(g.Start).dt.floor('h') g=g.set_index('time').join(w,how='inner') if len(g)<5000: continue diff=np.abs((g.wd-r.bearing+180)%360-180) g['sector']=np.where(diff<=45,'port',np.where(diff>=90,'away','side')); g['ym']=g.index.to_period('M') a=g.groupby(['ym','sector']).v.agg(['mean','size']).unstack('sector') a.columns=['%s_%s'%(x,y) for x,y in a.columns]; a=a.reset_index(); a['sp']=sp; rows.append(a) A=pd.concat(rows,ignore_index=True) A['excess']=A.mean_port-A.mean_away; A=A.merge(s.reset_index()[['sp','country','basin','d_port','portname','port_calls','stype','area','bearing','met0','met_d0']],on='sp') A.to_parquet('eea/wind_month.parquet'); print('nokta',A.sp.nunique(),A.drop_duplicates('sp').basin.value_counts().to_dict(),'ay',len(A)) print(A.groupby([A.ym.dt.year,'basin'])[['mean_port','mean_away','excess']].mean().unstack().round(2).to_string()) ''' 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, '-W', 'ignore', _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')