scenario comparison/pathways_choice/.ipynb_checkpoints/eclipse_distance-checkpoint.ipynb

+%% Cell type:markdown id:839642a1-ba5e-4fec-beba-565a295541ce tags:
+
+# Pathways and Eclipse
+
+%% Cell type:code id:a6114cac-703c-49e7-8ecd-fec6dbbc95be tags:
+
+``` python
+import xarray as xr
+import numpy as np
+import matplotlib.pyplot as plt
+import pandas as pd
+import glob
+from scipy import stats
+import cftime
+from datetime import datetime
+```
+
+%% Cell type:code id:dc2542b4-60ec-403d-9cf3-34e0582b12e2 tags:
+
+``` python
+specie = 'BC'
+```
+
+%% Cell type:code id:f0285a1d-0acc-4899-ba70-348dc8d7a507 tags:
+
+``` python
+coords_eu = {
+    'lon1': -10,
+    'lon2': 25,
+    'lat1': 35,
+    'lat2': 65
+}
+```
+
+%% Cell type:code id:9fc5f1fb-dbf7-49cb-995c-0207e249dbaa tags:
+
+``` python
+coords_wa = {
+    'lon1': -20,
+    'lon2': 20,
+    'lat1': -5,
+    'lat2': 20
+}
+```
+
+%% Cell type:code id:5be36c60-d718-4e2d-b381-4d38321d7f76 tags:
+
+``` python
+def region_total(emissions, coords):
+    emissions_regional = select_region(emissions, coords)
+    return emissions_regional.sum('lat').sum('lon')
+
+def select_region(data, coords):
+    return data.sel(lon=slice(coords['lon1'], coords['lon2'])).sel(lat=slice(coords['lat1'], coords['lat2']))
+```
+
+%% Cell type:code id:2ae45db9-249c-4004-8108-98361e11f8ba tags:
+
+``` python
+path_ssp126 = '/work/bd1063/b309253/emissions/regional/SSPs_Glb_0.5x0.5_anthro_'+specie+'_SSP126_yearly.nc'
+path_ssp245 = '/work/bd1063/b309253/emissions/regional/SSPs_Glb_0.5x0.5_anthro_'+specie+'_SSP245_yearly.nc'
+path_ssp370 = '/work/bd1063/b309253/emissions/regional/SSPs_Glb_0.5x0.5_anthro_'+specie+'_SSP370_yearly.nc'
+path_ssp460 = '/work/bd1063/b309253/emissions/regional/SSPs_Glb_0.5x0.5_anthro_'+specie+'_SSP460_yearly.nc'
+path_ssp585 = '/work/bd1063/b309253/emissions/regional/SSPs_Glb_0.5x0.5_anthro_'+specie+'_SSP585_yearly.nc'
+```
+
+%% Cell type:code id:c9222316-e2d7-42fa-80b9-e4a42c36e84c tags:
+
+``` python
+ssp126_eu = region_total(xr.open_dataset(path_ssp126)['s3'], coords_eu)
+ssp245_eu = region_total(xr.open_dataset(path_ssp245)['s3'], coords_eu)
+ssp370_eu = region_total(xr.open_dataset(path_ssp370)['s3'], coords_eu)
+ssp460_eu = region_total(xr.open_dataset(path_ssp460)['s3'], coords_eu)
+ssp585_eu = region_total(xr.open_dataset(path_ssp585)['s3'], coords_eu)
+```
+
+%% Cell type:code id:f9614397-43c7-42ce-b8e6-538f6846ce2d tags:
+
+``` python
+ssp126_wa = region_total(xr.open_dataset(path_ssp126)['s3'], coords_wa)
+ssp245_wa = region_total(xr.open_dataset(path_ssp245)['s3'], coords_wa)
+ssp370_wa = region_total(xr.open_dataset(path_ssp370)['s3'], coords_wa)
+ssp460_wa = region_total(xr.open_dataset(path_ssp460)['s3'], coords_wa)
+ssp585_wa = region_total(xr.open_dataset(path_ssp585)['s3'], coords_wa)
+```
+
+%% Cell type:code id:2e6fcbda-57e8-47e1-adf3-10e3e9542441 tags:
+
+``` python
+eclipse6_path = '/work/bd1063/b309253/emissions/regional/ECLIPSE-V6_Glb_0.5x0.5_anthro_BC_CLE_yearly.nc'
+```
+
+%% Cell type:code id:7a2ec56f-a7da-4509-aee9-63c2d1eb14a9 tags:
+
+``` python
+eclipse6_eu = region_total(xr.open_dataset(eclipse6_path)['emis_tra'], coords_eu)
+eclipse6_wa = region_total(xr.open_dataset(eclipse6_path)['emis_tra'], coords_wa)
+```
+
+%% Cell type:code id:2495889a-d1ff-44d4-b66e-3f699fc0e51e tags:
+
+``` python
+eclipse6_eu['time'] = [cftime.DatetimeNoLeap(*date.astype('datetime64[s]').astype(datetime).timetuple()[:5]) for date in eclipse6_eu['time'].values]
+eclipse6_wa['time'] = [cftime.DatetimeNoLeap(*date.astype('datetime64[s]').astype(datetime).timetuple()[:5]) for date in eclipse6_wa['time'].values]
+```
+
+%% Cell type:code id:01252014-a172-4b7b-a3c1-01b73563ce41 tags:
+
+``` python
+eclipse6_eu.plot(label='ECLIPSE v6 EU')
+eclipse6_wa.plot(label='ECLIPSE v6 WA')
+
+plt.title('Yearly '+specie+' emissions in ECLIPSE in Europe and West Africa')
+plt.ylabel('Tg/a')
+plt.xlabel('years')
+plt.legend(bbox_to_anchor=(1.0, 1.0))
+```
+
+%% Output
+
+    <matplotlib.legend.Legend at 0x7fff8dd81600>
+
+
+
+%% Cell type:code id:bfbac69c-d053-4887-a2f1-078ce5ec10f6 tags:
+
+``` python
+ssp126_eu.plot(label='SSP126')
+ssp245_eu.plot(label='SSP245')
+ssp370_eu.plot(label='SSP370')
+ssp460_eu.plot(label='SSP460')
+ssp585_eu.plot(label='SSP585')
+eclipse6_eu.plot(label='ECLIPSE v6 EU')
+
+plt.title('Yearly '+specie+' emissions in available SSPs in Europe')
+plt.ylabel('Tg/a')
+plt.xlabel('years')
+plt.legend(bbox_to_anchor=(1.0, 1.0))
+```
+
+%% Output
+
+    <matplotlib.legend.Legend at 0x7fff8df72e90>
+
+
+
+%% Cell type:code id:c24f5583-d9ed-45a8-8d09-214391bfb2e7 tags:
+
+``` python
+ssp126_wa.plot(label='SSP126')
+ssp245_wa.plot(label='SSP245')
+ssp370_wa.plot(label='SSP370')
+ssp460_wa.plot(label='SSP460')
+ssp585_wa.plot(label='SSP585')
+eclipse6_wa.plot(label='ECLIPSE v6 EU')
+
+plt.title('Yearly '+specie+' emissions in available SSPs in West Africa')
+plt.ylabel('Tg/a')
+plt.xlabel('years')
+plt.legend(bbox_to_anchor=(1.0, 1.0))
+```
+
+%% Output
+
+    <matplotlib.legend.Legend at 0x7fff8dc37130>
+
+
+
+%% Cell type:markdown id:0054ae8c-e9cd-4139-9d8f-7dbe44a06ba7 tags:
+
+## Which SSP scenario is closer to ECLIPSE v6?
+
+%% Cell type:code id:fe185077-b64c-4ab4-8c31-da651ca72e0c tags:
+
+``` python
+common_time_range = slice(max(ssp126_eu['time'].min(), eclipse6_eu['time'].min()),
+                          min(ssp126_eu['time'].max(), eclipse6_eu['time'].max()))
+```
+
+%% Cell type:code id:72b86f44-4fb0-4199-8870-b99162ec1666 tags:
+
+``` python
+def interpolate(dataset):
+    new_time =  xr.cftime_range(start=dataset['time'].values[0], end=dataset['time'].values[-1], freq='1Y')
+    empty_da = xr.DataArray(data=None, dims=['time'], coords={'time': new_time})
+    empty_ds = empty_da.to_dataset(name='emissions')
+
+    empty_ds['emissions'] = dataset.interp(time=empty_ds['time'])
+    return empty_ds.sel(time=common_time_range)['emissions']
+```
+
+%% Cell type:code id:b61aa9d6-fca0-4deb-b4b5-0700a3d9978e tags:
+
+``` python
+eclipse6_eu_interp = interpolate(eclipse6_eu)
+eclipse6_wa_interp = interpolate(eclipse6_wa)
+```
+
+%% Cell type:code id:46a0a0ac-067d-434a-99d8-217c2f35c2c4 tags:
+
+``` python
+ssp126_eu_interp = interpolate(ssp126_eu)
+ssp245_eu_interp = interpolate(ssp245_eu)
+ssp370_eu_interp = interpolate(ssp370_eu)
+ssp460_eu_interp = interpolate(ssp460_eu)
+ssp585_eu_interp = interpolate(ssp585_eu)
+```
+
+%% Cell type:code id:f52c8fd4-71d3-4b1d-86be-7287ac70773f tags:
+
+``` python
+ssp126_wa_interp = interpolate(ssp126_wa)
+ssp245_wa_interp = interpolate(ssp245_wa)
+ssp370_wa_interp = interpolate(ssp370_wa)
+ssp460_wa_interp = interpolate(ssp460_wa)
+ssp585_wa_interp = interpolate(ssp585_wa)
+```
+
+%% Cell type:markdown id:1089cb56-cecf-4a3d-8978-66b50ade309b tags:
+
+## Results
+
+%% Cell type:markdown id:be50d26b-886b-4268-87fc-d697c2ed1a7a tags:
+
+## BC
+
+%% Cell type:code id:ae584401-d35b-4ac3-8f72-d86f7014c3f1 tags:
+
+``` python
+print('ssp126 eu',stats.pearsonr(ssp126_eu_interp, eclipse6_eu_interp)[0])
+print('ssp245 eu',stats.pearsonr(ssp245_eu_interp, eclipse6_eu_interp)[0])
+print('ssp370 eu',stats.pearsonr(ssp370_eu_interp, eclipse6_eu_interp)[0])
+print('ssp460 eu',stats.pearsonr(ssp460_eu_interp, eclipse6_eu_interp)[0])
+print('ssp585 eu',stats.pearsonr(ssp585_eu_interp, eclipse6_eu_interp)[0])
+```
+
+%% Output
+
+    ssp126 eu 0.9838082653687723
+    ssp245 eu 0.9958802454937001
+    ssp370 eu 0.9282812000005504
+    ssp460 eu 0.9097726137748288
+    ssp585 eu 0.9584073024833202
+
+%% Cell type:code id:18ddce31-72d4-4b34-8882-42529c3538e7 tags:
+
+``` python
+print('ssp126 wa',stats.pearsonr(ssp126_wa_interp, eclipse6_wa_interp)[0])
+print('ssp245 wa',stats.pearsonr(ssp245_wa_interp, eclipse6_wa_interp)[0])
+print('ssp370 wa',stats.pearsonr(ssp370_wa_interp, eclipse6_wa_interp)[0])
+print('ssp460 wa',stats.pearsonr(ssp460_wa_interp, eclipse6_wa_interp)[0])
+print('ssp585 wa',stats.pearsonr(ssp585_wa_interp, eclipse6_wa_interp)[0])
+```
+
+%% Output
+
+    ssp126 wa 0.7661857653639685
+    ssp245 wa 0.9904829493369731
+    ssp370 wa 0.9914439475002073
+    ssp460 wa 0.7414116894203907
+    ssp585 wa 0.7924238066121922
+
+%% Cell type:code id:5f258fb4-0e68-4bae-a135-01101b119810 tags:
+
+``` python
+print('ssp126 eu', np.mean(np.abs(eclipse6_eu_interp - ssp126_eu_interp)).values)
+print('ssp245 eu', np.mean(np.abs(eclipse6_eu_interp - ssp245_eu_interp)).values)
+print('ssp370 eu', np.mean(np.abs(eclipse6_eu_interp - ssp370_eu_interp)).values)
+print('ssp460 eu', np.mean(np.abs(eclipse6_eu_interp - ssp460_eu_interp)).values)
+print('ssp585 eu', np.mean(np.abs(eclipse6_eu_interp - ssp585_eu_interp)).values)
+```
+
+%% Output
+
+    ssp126 eu 0.012893578471750778
+    ssp245 eu 0.006442824651262997
+    ssp370 eu 0.022049087697454162
+    ssp460 eu 0.012301613827011825
+    ssp585 eu 0.01274909062937281
+
+%% Cell type:code id:55a0c1d0-9706-47be-942d-cf27e5498c1e tags:
+
+``` python
+print('ssp126 wa', np.mean(np.abs(eclipse6_wa_interp - ssp126_wa_interp)).values)
+print('ssp245 wa', np.mean(np.abs(eclipse6_wa_interp - ssp245_wa_interp)).values)
+print('ssp370 wa', np.mean(np.abs(eclipse6_wa_interp - ssp370_wa_interp)).values)
+print('ssp460 wa', np.mean(np.abs(eclipse6_wa_interp - ssp460_wa_interp)).values)
+print('ssp585 wa', np.mean(np.abs(eclipse6_wa_interp - ssp585_wa_interp)).values)
+```
+
+%% Output
+
+    ssp126 wa 0.01259051129099441
+    ssp245 wa 0.0053062339464891445
+    ssp370 wa 0.0062918563501475604
+    ssp460 wa 0.00735905203764161
+    ssp585 wa 0.01394626363053364
+
+%% Cell type:code id:488e09f8-14da-4696-9421-fd24442a4c94 tags:
+
+``` python
+```