scenario comparison/.ipynb_checkpoints/pathways-checkpoint.ipynb

+%% Cell type:markdown id:6a46baa6-8954-43ad-9528-a6561ae50130 tags:
+
+# Pathway choice
+
+%% Cell type:markdown id:86d94dbb-8b96-41b4-8d35-82d7829ab214 tags:
+
+Here we explore how to select adequate scenarios for our simulations
+
+%% Cell type:code id:c084c3d5-bd57-42bb-b8dd-c79d528fd21f 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
+```
+
+%% Cell type:code id:568ff7d2-1e7f-4999-a943-77eeef4d9f12 tags:
+
+``` python
+specie = 'SO2'
+```
+
+%% Cell type:code id:60940ffb-c3a9-47f3-99e2-a06809a68698 tags:
+
+``` python
+coords_eu = {
+    'lon1': -10,
+    'lon2': 25,
+    'lat1': 35,
+    'lat2': 65
+}
+```
+
+%% Cell type:code id:db61d6bd-ce0a-4acd-ae92-46271ea08017 tags:
+
+``` python
+coords_wa = {
+    'lon1': -20,
+    'lon2': 20,
+    'lat1': -5,
+    'lat2': 20
+}
+```
+
+%% Cell type:code id:effb9f30-db8a-4201-80ae-304ad08c5bef 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:6d069dd3-3429-4acd-8e21-83e7d9dfca23 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:e93feca7-4814-42e0-90d8-6254a04c44d1 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:f5b20a04-b58f-47ab-884d-2298bf8f8671 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:3b54b088-8b9c-4677-ae5b-074615703166 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')
+
+plt.title('Yearly '+specie+' emissions in available SSPs in ' )
+plt.ylabel('Tg/a')
+plt.xlabel('years')
+plt.legend(bbox_to_anchor=(1.0, 1.0))
+```
+
+%% Output
+
+    <matplotlib.legend.Legend at 0x7fff9717bbb0>
+
+
+
+%% Cell type:code id:386da6af-40cb-4c38-8f78-3e63ec4560d8 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')
+
+plt.title('Yearly '+specie+' emissions in available SSPs in ' )
+plt.ylabel('Tg/a')
+plt.xlabel('years')
+plt.legend(bbox_to_anchor=(1.0, 1.0))
+```
+
+%% Output
+
+    <matplotlib.legend.Legend at 0x7fff972165f0>
+
+
+
+%% Cell type:code id:839c7257-a2c8-4fb8-9a29-169e1f7b045d tags:
+
+``` python
+path_ceds = '/work/bd1063/b309253/emissions/regional/CEDS_Glb_0.5x0.5_anthro_'+specie+'__yearly.nc'
+```
+
+%% Cell type:code id:86b45d42-9d01-46f9-aba1-f52f9ecb221e tags:
+
+``` python
+ceds_eu = region_total(xr.open_dataset(path_ceds)['transportation'], coords_eu)
+ceds_wa = region_total(xr.open_dataset(path_ceds)['transportation'], coords_wa)
+```
+
+%% Cell type:code id:58a6f6f5-7818-406c-a040-e0e2ca60c368 tags:
+
+``` python
+# Select the years that overlap with ssp projections
+ceds_recentpast_eu = ceds_eu.values[-5:]
+ceds_recentpast_wa = ceds_wa.values[-5:]
+```
+
+%% Cell type:markdown id:1084e1aa-a5ec-44b4-8660-a0799498edd9 tags:
+
+### We want to linearly 'fill in' the values of the projections between 2015 and 2019 to compare with the data from CEDS
+
+%% Cell type:code id:9d901dfa-419a-4bec-96d4-5082db215d6d tags:
+
+``` python
+def fill_values(ssp):
+    step = (ssp[1] - ssp[0])/5
+    return np.arange(ssp[0], ssp[1], step)
+```
+
+%% Cell type:code id:f78b0e7e-c273-428c-b2a5-8bb9fbcc6056 tags:
+
+``` python
+ssp126_recentpast_eu = fill_values(ssp126_eu)
+ssp245_recentpast_eu = fill_values(ssp245_eu)
+ssp370_recentpast_eu = fill_values(ssp370_eu)
+ssp460_recentpast_eu = fill_values(ssp460_eu)
+ssp585_recentpast_eu = fill_values(ssp585_eu)
+```
+
+%% Cell type:code id:5972a211-2a8f-45c3-a2e3-e0e0f16be001 tags:
+
+``` python
+ssp126_recentpast_wa = fill_values(ssp126_wa)
+ssp245_recentpast_wa = fill_values(ssp245_wa)
+ssp370_recentpast_wa = fill_values(ssp370_wa)
+ssp460_recentpast_wa = fill_values(ssp460_wa)
+ssp585_recentpast_wa = fill_values(ssp585_wa)
+```
+
+%% Cell type:code id:2f7d5474-6000-4c31-860c-9493a62bff05 tags:
+
+``` python
+print('ssp126 eu',stats.pearsonr(ceds_recentpast_eu, ssp126_recentpast_eu)[0])
+print('ssp245 eu',stats.pearsonr(ceds_recentpast_eu, ssp245_recentpast_eu)[0])
+print('ssp370 eu',stats.pearsonr(ceds_recentpast_eu, ssp370_recentpast_eu)[0])
+print('ssp460 eu',stats.pearsonr(ceds_recentpast_eu, ssp460_recentpast_eu)[0])
+print('ssp585 eu',stats.pearsonr(ceds_recentpast_eu, ssp585_recentpast_eu)[0])
+```
+
+%% Output
+
+    ssp126 eu -0.8880639532330912
+    ssp245 eu -0.8880639532330912
+    ssp370 eu 0.8880639532330912
+    ssp460 eu -0.8880639532330914
+    ssp585 eu -0.8880639532330914
+
+%% Cell type:code id:a0d11060-8e66-4e2e-828b-0c72b923c83e tags:
+
+``` python
+print('ssp126 wa',stats.pearsonr(ceds_recentpast_wa, ssp126_recentpast_wa)[0])
+print('ssp245 wa',stats.pearsonr(ceds_recentpast_wa, ssp245_recentpast_wa)[0])
+print('ssp370 wa',stats.pearsonr(ceds_recentpast_wa, ssp370_recentpast_wa)[0])
+print('ssp460 wa',stats.pearsonr(ceds_recentpast_wa, ssp460_recentpast_wa)[0])
+print('ssp585 wa',stats.pearsonr(ceds_recentpast_wa, ssp585_recentpast_wa)[0])
+```
+
+%% Output
+
+    ssp126 wa -0.9900271062887032
+    ssp245 wa -0.9900271062887032
+    ssp370 wa 0.990027106288703
+    ssp460 wa 0.9900271062887032
+    ssp585 wa 0.9900271062887032
+
+%% Cell type:markdown id:4acd08c1-01cc-43ad-b1a1-218358d7c757 tags:
+
+Since we set a linear trend for the projections, this tells us only if the trend is in the right direction
+
+%% Cell type:markdown id:ba60d5ca-12e5-4536-84b2-71d171c0502e tags:
+
+## MAE
+
+%% Cell type:code id:2f78a66e-f06e-4ebc-94b2-5295104f6f1a tags:
+
+``` python
+print('ssp126 eu', np.mean(np.abs(ceds_recentpast_eu - ssp126_recentpast_eu)))
+print('ssp245 eu', np.mean(np.abs(ceds_recentpast_eu - ssp245_recentpast_eu)))
+print('ssp370 eu', np.mean(np.abs(ceds_recentpast_eu - ssp370_recentpast_eu)))
+print('ssp460 eu', np.mean(np.abs(ceds_recentpast_eu - ssp460_recentpast_eu)))
+print('ssp585 eu', np.mean(np.abs(ceds_recentpast_eu - ssp585_recentpast_eu)))
+```
+
+%% Output
+
+    ssp126 wa 0.01895322799682616
+    ssp245 wa 0.01734040677547454
+    ssp370 wa 0.009080237150192261
+    ssp460 wa 0.013976210355758656
+    ssp585 wa 0.012740185856819153
+
+%% Cell type:code id:b2476d75-dd00-4847-8139-320e670c560e tags:
+
+``` python
+print('ssp126 wa', np.mean(np.abs(ceds_recentpast_wa - ssp126_recentpast_wa)))
+print('ssp245 wa', np.mean(np.abs(ceds_recentpast_wa - ssp245_recentpast_wa)))
+print('ssp370 wa', np.mean(np.abs(ceds_recentpast_wa - ssp370_recentpast_wa)))
+print('ssp460 wa', np.mean(np.abs(ceds_recentpast_wa - ssp460_recentpast_wa)))
+print('ssp585 wa', np.mean(np.abs(ceds_recentpast_wa - ssp585_recentpast_wa)))
+```
+
+%% Output
+
+    ssp126 wa 0.026111802458763135
+    ssp245 wa 0.047006168961524965
+    ssp370 wa 0.021110220253467572
+    ssp460 wa 0.019767563045024868
+    ssp585 wa 0.02138192802667619

scenario comparison/catalogues_comparisons/.ipynb_checkpoints/pathways-checkpoint.ipynb

-{
- "cells": [],
- "metadata": {},
- "nbformat": 4,
- "nbformat_minor": 5
-}