From e180931ce1aeed2df8a340836cda75ff35f00029 Mon Sep 17 00:00:00 2001 From: martavp Date: Thu, 14 Jan 2021 10:06:29 +0100 Subject: [PATCH] Move ploting of carbon budget distribution to plot_summary.py --- scripts/plot_summary.py | 151 ++++++++++++++++++++++++++++-- scripts/prepare_sector_network.py | 151 +++--------------------------- 2 files changed, 159 insertions(+), 143 deletions(-) diff --git a/scripts/plot_summary.py b/scripts/plot_summary.py index 7f37159c..fe28cfed 100644 --- a/scripts/plot_summary.py +++ b/scripts/plot_summary.py @@ -1,6 +1,6 @@ - +import numpy as np import pandas as pd #allow plotting without Xwindows @@ -9,7 +9,7 @@ matplotlib.use('Agg') import matplotlib.pyplot as plt - +from prepare_sector_network import co2_emissions_year #consolidate and rename def rename_techs(label): @@ -237,7 +237,137 @@ def plot_balances(): fig.savefig(snakemake.output.balances[:-10] + k + ".pdf",transparent=True) +def historical_emissions(cts): + """ + read historical emissions to add them to the carbon budget plot + """ + #https://www.eea.europa.eu/data-and-maps/data/national-emissions-reported-to-the-unfccc-and-to-the-eu-greenhouse-gas-monitoring-mechanism-16 + #downloaded 201228 (modified by EEA last on 201221) + fn = "data/eea/UNFCCC_v23.csv" + df = pd.read_csv(fn, encoding="latin-1") + df.loc[df["Year"] == "1985-1987","Year"] = 1986 + df["Year"] = df["Year"].astype(int) + df = df.set_index(['Year', 'Sector_name', 'Country_code', 'Pollutant_name']).sort_index() + e = pd.Series() + e["electricity"] = '1.A.1.a - Public Electricity and Heat Production' + e['residential non-elec'] = '1.A.4.b - Residential' + e['services non-elec'] = '1.A.4.a - Commercial/Institutional' + e['rail non-elec'] = "1.A.3.c - Railways" + e["road non-elec"] = '1.A.3.b - Road Transportation' + e["domestic navigation"] = "1.A.3.d - Domestic Navigation" + e['international navigation'] = '1.D.1.b - International Navigation' + e["domestic aviation"] = '1.A.3.a - Domestic Aviation' + e["international aviation"] = '1.D.1.a - International Aviation' + e['total energy'] = '1 - Energy' + e['industrial processes'] = '2 - Industrial Processes and Product Use' + e['agriculture'] = '3 - Agriculture' + e['LULUCF'] = '4 - Land Use, Land-Use Change and Forestry' + e['waste management'] = '5 - Waste management' + e['other'] = '6 - Other Sector' + e['indirect'] = 'ind_CO2 - Indirect CO2' + e["total wL"] = "Total (with LULUCF)" + e["total woL"] = "Total (without LULUCF)" + + pol = ["CO2"] # ["All greenhouse gases - (CO2 equivalent)"] + cts + if "GB" in cts: + cts.remove("GB") + cts.append("UK") + + year = np.arange(1990,2018).tolist() + + idx = pd.IndexSlice + co2_totals = df.loc[idx[year,e.values,cts,pol],"emissions"].unstack("Year").rename(index=pd.Series(e.index,e.values)) + + co2_totals = (1/1e6)*co2_totals.groupby(level=0, axis=0).sum() #Gton CO2 + + co2_totals.loc['industrial non-elec'] = co2_totals.loc['total energy'] - co2_totals.loc[['electricity', 'services non-elec','residential non-elec', 'road non-elec', + 'rail non-elec', 'domestic aviation', 'international aviation', 'domestic navigation', + 'international navigation']].sum() + + emissions = co2_totals.loc["electricity"] + if "T" in opts: + emissions += co2_totals.loc[[i+ " non-elec" for i in ["rail","road"]]].sum() + if "H" in opts: + emissions += co2_totals.loc[[i+ " non-elec" for i in ["residential","services"]]].sum() + if "I" in opts: + emissions += co2_totals.loc[["industrial non-elec","industrial processes", + "domestic aviation","international aviation", + "domestic navigation","international navigation"]].sum() + return emissions + + + +def plot_carbon_budget_distribution(): + """ + Plot historical carbon emissions in the EU and decarbonization path + """ + + import matplotlib.gridspec as gridspec + import seaborn as sns; sns.set() + sns.set_style('ticks') + plt.style.use('seaborn-ticks') + plt.rcParams['xtick.direction'] = 'in' + plt.rcParams['ytick.direction'] = 'in' + plt.rcParams['xtick.labelsize'] = 20 + plt.rcParams['ytick.labelsize'] = 20 + + plt.figure(figsize=(10, 7)) + gs1 = gridspec.GridSpec(1, 1) + ax1 = plt.subplot(gs1[0,0]) + ax1.set_ylabel('CO$_2$ emissions (Gt per year)',fontsize=22) + ax1.set_ylim([0,5]) + ax1.set_xlim([1990,snakemake.config['scenario']['planning_horizons'][-1]+1]) + + path_cb = snakemake.config['results_dir'] + snakemake.config['run'] + '/csvs/' + countries=pd.read_csv(path_cb + 'countries.csv', index_col=1) + cts=countries.index.to_list() + e_1990 = co2_emissions_year(cts, opts, year=1990) + CO2_CAP=pd.read_csv(path_cb + 'carbon_budget_distribution.csv', + index_col=0) + + + ax1.plot(e_1990*CO2_CAP[o],linewidth=3, + color='dodgerblue', label=None) + + emissions = historical_emissions(cts) + + ax1.plot(emissions, color='black', linewidth=3, label=None) + + #plot commited and uder-discussion targets + #(notice that historical emissions include all countries in the + # network, but targets refer to EU) + ax1.plot([2020],[0.8*emissions[1990]], + marker='*', markersize=12, markerfacecolor='black', + markeredgecolor='black') + + ax1.plot([2030],[0.45*emissions[1990]], + marker='*', markersize=12, markerfacecolor='white', + markeredgecolor='black') + + ax1.plot([2030],[0.6*emissions[1990]], + marker='*', markersize=12, markerfacecolor='black', + markeredgecolor='black') + + ax1.plot([2050, 2050],[x*emissions[1990] for x in [0.2, 0.05]], + color='gray', linewidth=2, marker='_', alpha=0.5) + + ax1.plot([2050],[0.01*emissions[1990]], + marker='*', markersize=12, markerfacecolor='white', + linewidth=0, markeredgecolor='black', + label='EU under-discussion target', zorder=10, + clip_on=False) + + ax1.plot([2050],[0.125*emissions[1990]],'ro', + marker='*', markersize=12, markerfacecolor='black', + markeredgecolor='black', label='EU commited target') + + ax1.legend(fancybox=True, fontsize=18, loc=(0.01,0.01), + facecolor='white', frameon=True) + + path_cb_plot = snakemake.config['results_dir'] + snakemake.config['run'] + '/graphs/' + plt.savefig(path_cb_plot+'carbon_budget_plot.pdf', dpi=300) if __name__ == "__main__": # Detect running outside of snakemake and mock snakemake for testing @@ -249,13 +379,16 @@ if __name__ == "__main__": snakemake.config = yaml.safe_load(f) snakemake.input = Dict() snakemake.output = Dict() - + snakemake.wildcards = Dict() + #snakemake.wildcards['sector_opts']='3H-T-H-B-I-solar3-dist1-cb48be3' + for item in ["costs", "energy"]: snakemake.input[item] = snakemake.config['summary_dir'] + '/{name}/csvs/{item}.csv'.format(name=snakemake.config['run'],item=item) snakemake.output[item] = snakemake.config['summary_dir'] + '/{name}/graphs/{item}.pdf'.format(name=snakemake.config['run'],item=item) - snakemake.input["balances"] = snakemake.config['summary_dir'] + '/test/csvs/supply_energy.csv' - snakemake.output["balances"] = snakemake.config['summary_dir'] + '/test/graphs/balances-energy.csv' - + snakemake.input["balances"] = snakemake.config['summary_dir'] + '/{name}/csvs/supply_energy.csv'.format(name=snakemake.config['run'],item=item) + snakemake.output["balances"] = snakemake.config['summary_dir'] + '/{name}/graphs/balances-energy.csv'.format(name=snakemake.config['run'],item=item) + + n_header = 4 plot_costs() @@ -263,3 +396,9 @@ if __name__ == "__main__": plot_energy() plot_balances() + + for sector_opts in snakemake.config['scenario']['sector_opts']: + opts=sector_opts.split('-') + for o in opts: + if "cb" in o: + plot_carbon_budget_distribution() \ No newline at end of file diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 088b124f..a2c2795e 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -48,7 +48,7 @@ override_component_attrs["Store"].loc["build_year"] = ["integer","year",np.nan," override_component_attrs["Store"].loc["lifetime"] = ["float","years",np.nan,"lifetime","Input (optional)"] -def co2_emissions_year(year): +def co2_emissions_year(cts, opts, year): """ calculate co2 emissions in one specific year (e.g. 1990 or 2018). """ @@ -62,10 +62,6 @@ def co2_emissions_year(year): eurostat_co2 = build_eurostat_co2(year) co2_totals=build_co2_totals(eea_co2, eurostat_co2, year) - - pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0) - pop_layout["ct"] = pop_layout.index.str[:2] - cts = pop_layout.ct.value_counts().index co2_emissions = co2_totals.loc[cts, "electricity"].sum() @@ -80,70 +76,6 @@ def co2_emissions_year(year): co2_emissions *=0.001 #MtCO2 to GtCO2 return co2_emissions -def historical_emissions(): - """ - read historical emissions to add them to the carbon budget plot - """ - #https://www.eea.europa.eu/data-and-maps/data/national-emissions-reported-to-the-unfccc-and-to-the-eu-greenhouse-gas-monitoring-mechanism-16 - #downloaded 201228 (modified by EEA last on 201221) - fn = "data/eea/UNFCCC_v23.csv" - df = pd.read_csv(fn, encoding="latin-1") - df.loc[df["Year"] == "1985-1987","Year"] = 1986 - df["Year"] = df["Year"].astype(int) - df = df.set_index(['Year', 'Sector_name', 'Country_code', 'Pollutant_name']).sort_index() - - e = pd.Series() - e["electricity"] = '1.A.1.a - Public Electricity and Heat Production' - e['residential non-elec'] = '1.A.4.b - Residential' - e['services non-elec'] = '1.A.4.a - Commercial/Institutional' - e['rail non-elec'] = "1.A.3.c - Railways" - e["road non-elec"] = '1.A.3.b - Road Transportation' - e["domestic navigation"] = "1.A.3.d - Domestic Navigation" - e['international navigation'] = '1.D.1.b - International Navigation' - e["domestic aviation"] = '1.A.3.a - Domestic Aviation' - e["international aviation"] = '1.D.1.a - International Aviation' - e['total energy'] = '1 - Energy' - e['industrial processes'] = '2 - Industrial Processes and Product Use' - e['agriculture'] = '3 - Agriculture' - e['LULUCF'] = '4 - Land Use, Land-Use Change and Forestry' - e['waste management'] = '5 - Waste management' - e['other'] = '6 - Other Sector' - e['indirect'] = 'ind_CO2 - Indirect CO2' - e["total wL"] = "Total (with LULUCF)" - e["total woL"] = "Total (without LULUCF)" - - pol = ["CO2"] # ["All greenhouse gases - (CO2 equivalent)"] - - - pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0) - pop_layout["ct"] = pop_layout.index.str[:2] - cts = pop_layout.ct.value_counts().index.to_list() - if "GB" in cts: - cts.remove("GB") - cts.append("UK") - - year = np.arange(1990,2018).tolist() - - idx = pd.IndexSlice - co2_totals = df.loc[idx[year,e.values,cts,pol],"emissions"].unstack("Year").rename(index=pd.Series(e.index,e.values)) - - co2_totals = (1/1e6)*co2_totals.groupby(level=0, axis=0).sum() #Gton CO2 - - co2_totals.loc['industrial non-elec'] = co2_totals.loc['total energy'] - co2_totals.loc[['electricity', 'services non-elec','residential non-elec', 'road non-elec', - 'rail non-elec', 'domestic aviation', 'international aviation', 'domestic navigation', - 'international navigation']].sum() - - emissions = co2_totals.loc["electricity"] - if "T" in opts: - emissions += co2_totals.loc[[i+ " non-elec" for i in ["rail","road"]]].sum() - if "H" in opts: - emissions += co2_totals.loc[[i+ " non-elec" for i in ["residential","services"]]].sum() - if "I" in opts: - emissions += co2_totals.loc[["industrial non-elec","industrial processes", - "domestic aviation","international aviation", - "domestic navigation","international navigation"]].sum() - return emissions - def build_carbon_budget(o): #distribute carbon budget following beta or exponential transition path @@ -155,11 +87,16 @@ def build_carbon_budget(o): #exponential decay carbon_budget = float(o[o.find("cb")+2:o.find("ex")]) r=float(o[o.find("ex")+2:]) - - e_1990 = co2_emissions_year(year=1990) + + + pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0) + pop_layout["ct"] = pop_layout.index.str[:2] + cts = pop_layout.ct.value_counts().index + + e_1990 = co2_emissions_year(cts, opts, year=1990) #emissions at the beginning of the path (last year available 2018) - e_0 = co2_emissions_year(year=2018) + e_0 = co2_emissions_year(cts, opts, year=2018) #emissions in 2019 and 2020 assumed equal to 2018 and substracted carbon_budget -= 2*e_0 planning_horizons = snakemake.config['scenario']['planning_horizons'] @@ -183,69 +120,9 @@ def build_carbon_budget(o): CO2_CAP.to_csv(path_cb + 'carbon_budget_distribution.csv', sep=',', line_terminator='\n', float_format='%.3f') - - """ - Plot historical carbon emissions in the EU and decarbonization path - """ - import matplotlib.pyplot as plt - import matplotlib.gridspec as gridspec - import seaborn as sns; sns.set() - sns.set_style('ticks') - plt.style.use('seaborn-ticks') - plt.rcParams['xtick.direction'] = 'in' - plt.rcParams['ytick.direction'] = 'in' - plt.rcParams['xtick.labelsize'] = 20 - plt.rcParams['ytick.labelsize'] = 20 - - plt.figure(figsize=(10, 7)) - gs1 = gridspec.GridSpec(1, 1) - ax1 = plt.subplot(gs1[0,0]) - ax1.set_ylabel('CO$_2$ emissions (Gt per year)',fontsize=22) - ax1.set_ylim([0,5]) - ax1.set_xlim([1990,planning_horizons[-1]+1]) - ax1.plot(e_1990*CO2_CAP[o],linewidth=3, - color='dodgerblue', label=None) - - emissions = historical_emissions() - - ax1.plot(emissions, color='black', linewidth=3, label=None) - - #plot commited and uder-discussion targets - #(notice that historical emissions include all countries in the - # network, but targets refer to EU) - ax1.plot([2020],[0.8*emissions[1990]], - marker='*', markersize=12, markerfacecolor='black', - markeredgecolor='black') - - ax1.plot([2030],[0.45*emissions[1990]], - marker='*', markersize=12, markerfacecolor='white', - markeredgecolor='black') - - ax1.plot([2030],[0.6*emissions[1990]], - marker='*', markersize=12, markerfacecolor='black', - markeredgecolor='black') - - ax1.plot([2050, 2050],[x*emissions[1990] for x in [0.2, 0.05]], - color='gray', linewidth=2, marker='_', alpha=0.5) - - ax1.plot([2050],[0.01*emissions[1990]], - marker='*', markersize=12, markerfacecolor='white', - linewidth=0, markeredgecolor='black', - label='EU under-discussion target', zorder=10, - clip_on=False) - - ax1.plot([2050],[0.125*emissions[1990]],'ro', - marker='*', markersize=12, markerfacecolor='black', - markeredgecolor='black', label='EU commited target') - - ax1.legend(fancybox=True, fontsize=18, loc=(0.01,0.01), - facecolor='white', frameon=True) - - path_cb_plot = snakemake.config['results_dir'] + snakemake.config['run'] + '/graphs/' - if not os.path.exists(path_cb_plot): - os.makedirs(path_cb_plot) - print('carbon budget distribution saved to ' + path_cb_plot + 'carbon_budget_plot.pdf') - plt.savefig(path_cb_plot+'carbon_budget_plot.pdf', dpi=300) + countries=pd.Series(data=cts) + countries.to_csv(path_cb + 'countries.csv', sep=',', + line_terminator='\n', float_format='%.3f') def add_lifetime_wind_solar(n): """ @@ -1985,7 +1862,7 @@ if __name__ == "__main__": snakemake = MockSnakemake( wildcards=dict(network='elec', simpl='', clusters='37', lv='1.0', opts='', planning_horizons='2020', - sector_opts='120H-T-H-B-I-solar3-dist1-cb40ex0'), + sector_opts='120H-T-H-B-I-solar3-dist1-cb48be3'), input=dict( network='../pypsa-eur/networks/{network}_s{simpl}_{clusters}_ec_lv{lv}_{opts}.nc', energy_totals_name='resources/energy_totals.csv', co2_totals_name='resources/co2_totals.csv', @@ -2024,7 +1901,7 @@ if __name__ == "__main__": retro_cost_energy = "resources/retro_cost_{network}_s{simpl}_{clusters}.csv", floor_area = "resources/floor_area_{network}_s{simpl}_{clusters}.csv" ), - output=['results/version-8/prenetworks/{network}_s{simpl}_{clusters}_lv{lv}__{sector_opts}_{planning_horizons}.nc'] + output=['results/version-cb48be3/prenetworks/{network}_s{simpl}_{clusters}_lv{lv}__{sector_opts}_{planning_horizons}.nc'] ) import yaml with open('config.yaml', encoding='utf8') as f: