diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs index 3f1edbd81..c85a3de7d 100644 --- a/.git-blame-ignore-revs +++ b/.git-blame-ignore-revs @@ -6,4 +6,5 @@ 5d1ef8a64055a039aa4a0834d2d26fe7752fe9a0 92080b1cd2ca5f123158571481722767b99c2b27 13769f90af4500948b0376d57df4cceaa13e78b5 +cc299f969e1535082bcf668e4afe206c57c8a920 9865a970893d9e515786f33c629b14f71645bf1e diff --git a/config/config.default.yaml b/config/config.default.yaml index 1033d49d7..89b650f19 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -44,7 +44,7 @@ scenario: opts: - '' sector_opts: - - Co2L0-3H-T-H-B-I-A-dist1 + - Co2L0-3h-T-H-B-I-A-dist1 planning_horizons: # - 2020 # - 2030 @@ -314,9 +314,9 @@ pypsa_eur: # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#energy energy: - energy_totals_year: 2011 + energy_totals_year: 2019 base_emissions_year: 1990 - eurostat_report_year: 2016 + eurostat_report_year: 2023 emissions: CO2 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#biomass @@ -501,6 +501,7 @@ sector: SMR_cc: true regional_methanol_demand: false regional_oil_demand: false + regional_coal_demand: false regional_co2_sequestration_potential: enable: false attribute: 'conservative estimate Mt' @@ -636,6 +637,15 @@ industry: 2040: 0.12 2045: 0.16 2050: 0.20 + sector_ratios_fraction_future: + 2020: 0.0 + 2025: 0.1 + 2030: 0.3 + 2035: 0.5 + 2040: 0.7 + 2045: 0.9 + 2050: 1.0 + basic_chemicals_without_NH3_production_today: 69. #Mt/a, = 86 Mtethylene-equiv - 17 MtNH3 HVC_production_today: 52. MWh_elec_per_tHVC_mechanical_recycling: 0.547 MWh_elec_per_tHVC_chemical_recycling: 6.9 @@ -815,6 +825,13 @@ plotting: energy_max: 20000 energy_min: -20000 energy_threshold: 50. + countries: + - all + carriers: + - electricity + - heat + carrier_groups: + electricity: [AC, low_voltage] nice_names: OCGT: "Open-Cycle Gas" diff --git a/doc/configtables/plotting.csv b/doc/configtables/plotting.csv index 82fc203c0..3124ce306 100644 --- a/doc/configtables/plotting.csv +++ b/doc/configtables/plotting.csv @@ -1,13 +1,16 @@ ,Unit,Values,Description map,,, -- boundaries,°,"[x1,x2,y1,y2]",Boundaries of the map plots in degrees latitude (y) and longitude (x) -projection,,,, --- name,--,"Valid Cartopy projection name","See https://scitools.org.uk/cartopy/docs/latest/reference/projections.html for list of available projections." +projection,,, +-- name,--,Valid Cartopy projection name,See https://scitools.org.uk/cartopy/docs/latest/reference/projections.html for list of available projections. -- args,--,--,"Other entries under 'projection' are passed as keyword arguments to the projection constructor, e.g. ``central_longitude: 10.``." costs_max,bn Euro,float,Upper y-axis limit in cost bar plots. costs_threshold,bn Euro,float,Threshold below which technologies will not be shown in cost bar plots. energy_max,TWh,float,Upper y-axis limit in energy bar plots. energy_min,TWh,float,Lower y-axis limit in energy bar plots. energy_threshold,TWh,float,Threshold below which technologies will not be shown in energy bar plots. -tech_colors,--,carrier -> HEX colour code,Mapping from network ``carrier`` to a colour (`HEX colour code `_). +countries,--,[str],List of countries you want to plot in your statistics figures. Default is all coutries but also single countries can be defined like ``DE``. +carriers,--,[str],Carrier you want to plot from in your statistics figures. This carrier must be an available bus carrier in your network or substring of a carrier or a defined carrier group. +carrier_groups,--,str -> carrier,Mapping from ``carrier_groups`` to individual bus carriers in the network. nice_names,--,str -> str,Mapping from network ``carrier`` to a more readable name. +tech_colors,--,carrier -> HEX colour code,Mapping from network ``carrier`` to a colour (`HEX colour code `_). diff --git a/doc/release_notes.rst b/doc/release_notes.rst index fce4ae1bd..baceed0ba 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -10,6 +10,8 @@ Release Notes Upcoming Release ================ +* Add plotting routine with statistics where we allow for plotting of individual energy carriers and countries. Besides the plots, we create all necessary csv files for the plotting routine. + * Merged two OPSD time series data versions into such that the option ``load: power_statistics:`` becomes superfluous and was hence removed. diff --git a/rules/build_sector.smk b/rules/build_sector.smk index c25c86739..bec9aa7aa 100644 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -433,6 +433,30 @@ rule build_industry_sector_ratios: "../scripts/build_industry_sector_ratios.py" +rule build_industry_sector_ratios_intermediate: + params: + industry=config["industry"], + input: + industry_sector_ratios=RESOURCES + "industry_sector_ratios.csv", + industrial_energy_demand_per_country_today=RESOURCES + + "industrial_energy_demand_per_country_today.csv", + industrial_production_per_country=RESOURCES + + "industrial_production_per_country.csv", + output: + industry_sector_ratios=RESOURCES + "industry_sector_ratios_{planning_horizons}.csv", + threads: 1 + resources: + mem_mb=1000, + log: + LOGS + "build_industry_sector_ratios_{planning_horizons}.log", + benchmark: + BENCHMARKS + "build_industry_sector_ratios_{planning_horizons}" + conda: + "../envs/environment.yaml" + script: + "../scripts/build_industry_sector_ratios_intermediate.py" + + rule build_industrial_production_per_country: params: industry=config["industry"], @@ -535,7 +559,7 @@ rule build_industrial_production_per_node: rule build_industrial_energy_demand_per_node: input: - industry_sector_ratios=RESOURCES + "industry_sector_ratios.csv", + industry_sector_ratios=RESOURCES + "industry_sector_ratios_{planning_horizons}.csv", industrial_production_per_node=RESOURCES + "industrial_production_elec_s{simpl}_{clusters}_{planning_horizons}.csv", industrial_energy_demand_per_node_today=RESOURCES @@ -566,7 +590,6 @@ rule build_industrial_energy_demand_per_country_today: industry=config["industry"], input: jrc="data/bundle-sector/jrc-idees-2015", - ammonia_production=RESOURCES + "ammonia_production.csv", industrial_production_per_country=RESOURCES + "industrial_production_per_country.csv", output: diff --git a/rules/collect.smk b/rules/collect.smk index dc0a94cc6..caa46e9bc 100644 --- a/rules/collect.smk +++ b/rules/collect.smk @@ -81,3 +81,22 @@ rule validate_elec_networks: **config["scenario"], kind=["production", "prices", "cross_border"], ), + + +rule plot_statistics: + input: + [ + expand( + RESULTS + + "statistics/figures/comparison/country_{country}/.statistics_{carrier}_plots", + country=config["plotting"].get("countries", "all"), + carrier=config["plotting"].get("carriers", ["all"]), + ), + expand( + RESULTS + + "statistics/figures/single/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}/country_{country}/.statistics_{carrier}_plots", + **config["scenario"], + country=config["plotting"].get("countries", "all"), + carrier=config["plotting"].get("carriers", ["all"]), + ), + ], diff --git a/rules/common.smk b/rules/common.smk index 5aa7ae53a..6f156455e 100644 --- a/rules/common.smk +++ b/rules/common.smk @@ -65,8 +65,11 @@ def has_internet_access(url="www.zenodo.org") -> bool: def input_eurostat(w): # 2016 includes BA, 2017 does not - report_year = config["energy"]["eurostat_report_year"] - return f"data/bundle-sector/eurostat-energy_balances-june_{report_year}_edition" + if config["energy"]["eurostat_report_year"] != 2023: + report_year = config["energy"]["eurostat_report_year"] + return f"data/bundle-sector/eurostat-energy_balances-june_{report_year}_edition" + else: + return "data/bundle-sector/eurostat-energy_balances-april_2023_edition" def solved_previous_horizon(wildcards): diff --git a/rules/postprocess.smk b/rules/postprocess.smk index 98b90fd14..cf49209a3 100644 --- a/rules/postprocess.smk +++ b/rules/postprocess.smk @@ -280,6 +280,92 @@ STATISTICS_BARPLOTS = [ "market_value", ] +STATISTICS = { + "capacity_factor": ("-", "p.u."), + "installed_capacity": (1e3, "GW"), + "optimal_capacity": (1e3, "GW"), + "capex": (1e9, "bn €"), + "opex": (1e9, "bn €"), + "total_cost": ("1e9", "bn €"), + "curtailment": (1e3, "GWh"), + "supply": (1e6, "TWh"), + "withdrawal": (1e6, "TWh"), + "energy_balance": (1e6, "TWh"), + "market_value": ("-", "€/MWh"), +} + + +rule save_statistics_csv: + params: + statistics=STATISTICS, + input: + network=RESULTS + + "postnetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc", + output: + **{ + f"{csv}": RESULTS + + "statistics/csv/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}/country_{country}/{carrier}_" + + f"{csv}.csv" + for carrier in config["plotting"].get("carriers", "all") + for csv in STATISTICS + }, + csv_touch=RESULTS + + "statistics/csv/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}/country_{country}/.statistics_{carrier}_csv", + script: + "../scripts/write_statistics.py" + + +rule plot_statistics_single: + params: + plotting=config["plotting"], + statistics=STATISTICS, + input: + **{ + f"{csv}": RESULTS + + "statistics/csv/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}/country_{country}/{carrier}_" + + f"{csv}.csv" + for carrier in config["plotting"].get("carriers", "all") + for csv in STATISTICS + }, + output: + **{ + f"{plot}": RESULTS + + "statistics/figures/single/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}/country_{country}/{carrier}_" + + f"{plot}.pdf" + for carrier in config["plotting"].get("carriers", "all") + for plot in STATISTICS + }, + barplots_touch=RESULTS + + "statistics/figures/single/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}/country_{country}/.statistics_{carrier}_plots", + script: + "../scripts/plot_statistics_single.py" + + +rule plot_statistics_comparison: + params: + plotting=config["plotting"], + statistics=STATISTICS, + input: + expand( + RESULTS + + "statistics/csv/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}/country_{country}/{carrier}_{csv}.csv", + **config["scenario"], + csv=STATISTICS, + allow_missing=True, + ), + output: + **{ + f"{plot}": RESULTS + + "statistics/figures/comparison/country_{country}/{carrier}_" + + f"{plot}.pdf" + for carrier in config["plotting"].get("carriers", "all") + for plot in STATISTICS + }, + barplots_touch=RESULTS + + "statistics/figures/comparison/country_{country}/.statistics_{carrier}_plots", + script: + "../scripts/plot_statistics_comparison.py" + rule plot_elec_statistics: params: diff --git a/rules/solve_myopic.smk b/rules/solve_myopic.smk index e8817de64..49df9d1ee 100644 --- a/rules/solve_myopic.smk +++ b/rules/solve_myopic.smk @@ -3,6 +3,40 @@ # SPDX-License-Identifier: MIT +rule build_existing_heating_distribution: + params: + baseyear=config["scenario"]["planning_horizons"][0], + sector=config["sector"], + existing_capacities=config["existing_capacities"], + input: + existing_heating="data/existing_infrastructure/existing_heating_raw.csv", + clustered_pop_layout=RESOURCES + "pop_layout_elec_s{simpl}_{clusters}.csv", + clustered_pop_energy_layout=RESOURCES + + "pop_weighted_energy_totals_s{simpl}_{clusters}.csv", + district_heat_share=RESOURCES + + "district_heat_share_elec_s{simpl}_{clusters}_{planning_horizons}.csv", + output: + existing_heating_distribution=RESOURCES + + "existing_heating_distribution_elec_s{simpl}_{clusters}_{planning_horizons}.csv", + wildcard_constraints: + planning_horizons=config["scenario"]["planning_horizons"][0], #only applies to baseyear + threads: 1 + resources: + mem_mb=2000, + log: + LOGS + + "build_existing_heating_distribution_elec_s{simpl}_{clusters}_{planning_horizons}.log", + benchmark: + ( + BENCHMARKS + + "build_existing_heating_distribution/elec_s{simpl}_{clusters}_{planning_horizons}" + ) + conda: + "../envs/environment.yaml" + script: + "../scripts/build_existing_heating_distribution.py" + + rule add_existing_baseyear: params: baseyear=config["scenario"]["planning_horizons"][0], diff --git a/scripts/build_ammonia_production.py b/scripts/build_ammonia_production.py index 1bcdf9ae5..77b330750 100644 --- a/scripts/build_ammonia_production.py +++ b/scripts/build_ammonia_production.py @@ -7,6 +7,7 @@ """ import country_converter as coco +import numpy as np import pandas as pd cc = coco.CountryConverter() @@ -30,8 +31,10 @@ ammonia.index = cc.convert(ammonia.index, to="iso2") years = [str(i) for i in range(2013, 2018)] - countries = ammonia.index.intersection(snakemake.params.countries) - ammonia = ammonia.loc[countries, years].astype(float) + + ammonia = ammonia[years] + ammonia.replace("--", np.nan, inplace=True) + ammonia = ammonia.astype(float) # convert from ktonN to ktonNH3 ammonia *= 17 / 14 diff --git a/scripts/build_energy_totals.py b/scripts/build_energy_totals.py index c67bb49d6..2d7ca6427 100644 --- a/scripts/build_energy_totals.py +++ b/scripts/build_energy_totals.py @@ -16,6 +16,7 @@ import pandas as pd from _helpers import mute_print from tqdm import tqdm +import os cc = coco.CountryConverter() logger = logging.getLogger(__name__) @@ -120,36 +121,92 @@ def build_eurostat(input_eurostat, countries, report_year, year): """ Return multi-index for all countries' energy data in TWh/a. """ - filenames = { + if report_year != 2023: + filenames = { 2016: f"/{year}-Energy-Balances-June2016edition.xlsx", 2017: f"/{year}-ENERGY-BALANCES-June2017edition.xlsx", } - with mute_print(): - dfs = pd.read_excel( - input_eurostat + filenames[report_year], - sheet_name=None, - skiprows=1, - index_col=list(range(4)), - ) - - # sorted_index necessary for slicing - lookup = eurostat_codes - labelled_dfs = { - lookup[df.columns[0]]: df - for df in dfs.values() - if lookup[df.columns[0]] in countries - } - df = pd.concat(labelled_dfs, sort=True).sort_index() - - # drop non-numeric and country columns - non_numeric_cols = df.columns[df.dtypes != float] - country_cols = df.columns.intersection(lookup.keys()) - to_drop = non_numeric_cols.union(country_cols) - df.drop(to_drop, axis=1, inplace=True) + with mute_print(): + dfs = pd.read_excel( + input_eurostat + filenames[report_year], + sheet_name=None, + skiprows=1, + index_col=list(range(4)), + ) - # convert ktoe/a to TWh/a - df *= 11.63 / 1e3 + # sorted_index necessary for slicing + lookup = eurostat_codes + labelled_dfs = { + lookup[df.columns[0]]: df + for df in dfs.values() + if lookup[df.columns[0]] in countries + } + df = pd.concat(labelled_dfs, sort=True).sort_index() + # drop non-numeric and country columns + non_numeric_cols = df.columns[df.dtypes != float] + country_cols = df.columns.intersection(lookup.keys()) + to_drop = non_numeric_cols.union(country_cols) + df.drop(to_drop, axis=1, inplace=True) + + # convert ktoe/a to TWh/a + df *= 11.63 / 1e3 + + else: + # read in every country file in countries + eurostat = pd.DataFrame() + countries = [country if country != 'GB' else 'UK' for country in countries] + for country in countries: + filename = f"/{country}-Energy-balance-sheets-April-2023-edition.xlsb" + if os.path.exists(input_eurostat + filename): + df = pd.read_excel( + input_eurostat + filename, + engine='pyxlsb', + sheet_name=str(year), + skiprows=4, + index_col=list(range(4))) + # replace entry 'Z' with 0 + df.replace('Z', 0, inplace=True) + # write 'International aviation' to the 2nd level of the multiindex + index_number = (df.index.get_level_values(1) == 'International aviation').argmax() + new_index = ('-', 'International aviation', 'International aviation', 'ktoe') + modified_index = list(df.index) + modified_index[index_number] = new_index + df.index = pd.MultiIndex.from_tuples(modified_index, names=df.index.names) + # drop the annoying subhead line + df.drop(df[df[year] == year].index, inplace=True) + # replace 'Z' with 0 + df = df.replace('Z', 0) + # add country to the multiindex + new_tuple = [(country, *idx) for idx in df.index] + new_mindex = pd.MultiIndex.from_tuples(new_tuple, names=['country', None, 'name', None, 'unit']) + df.index = new_mindex + # make numeric values where possible + df = df.apply(pd.to_numeric, errors='coerce') + # drop non-numeric columns + non_numeric_cols = df.columns[df.dtypes != float] + df.drop(non_numeric_cols, axis=1, inplace=True) + # concatenate the dataframes + eurostat = pd.concat([eurostat, df], axis=0) + + eurostat.drop(["Unnamed: 4", year, "Unnamed: 6"], axis=1, inplace=True) + # Renaming some indices + rename = { + 'Households': 'Residential', + 'Commercial & public services': 'Services', + 'Domestic navigation': 'Domestic Navigation' + } + for name, rename in rename.items(): + eurostat.index = eurostat.index.set_levels( + eurostat.index.levels[3].where(eurostat.index.levels[3] != name, rename), + level=3) + new_index = eurostat.index.set_levels(eurostat.index.levels[2].where(eurostat.index.levels[2] != 'International maritime bunkers', 'Bunkers'), level=2) + eurostat.index = new_index + + eurostat.rename(columns={'Total': 'Total all products'}, inplace=True) + eurostat.index = eurostat.index.set_levels(eurostat.index.levels[0].where(eurostat.index.levels[0] != 'UK', 'GB'), level=0) + + df = eurostat * 11.63 / 1e3 return df @@ -736,6 +793,129 @@ def build_transport_data(countries, population, idees): return transport_data +def rescale(idees_countries, energy, eurostat): + ''' + Takes JRC IDEES data from 2015 and rescales it by the ratio of the + eurostat data and the 2015 eurostat data. + missing data: ['passenger car efficiency', 'passenger cars'] + ''' + # read in the eurostat data for 2015 + eurostat_2015 = build_eurostat(input_eurostat, countries, 2023, 2015)[["Total all products", "Electricity"]] + eurostat_year = eurostat[["Total all products", "Electricity"]] + # calculate the ratio of the two data sets + ratio = eurostat_year / eurostat_2015 + ratio = ratio.droplevel([1,4]) + ratio.rename(columns={"Total all products": "total", "Electricity": "ele"}, inplace=True) + + residential_total = [ + "total residential space", + "total residential water", + "total residential cooking", + "total residential", + "derived heat residential", + "thermal uses residential", + ] + residential_ele = [ + "electricity residential space", + "electricity residential water", + "electricity residential cooking", + "electricity residential", + ] + + service_total = [ + "total services space", + "total services water", + "total services cooking", + "total services", + "derived heat services", + "thermal uses services", + ] + service_ele = [ + "electricity services space", + "electricity services water", + "electricity services cooking", + "electricity services", + ] + + agri_total = [ + "total agriculture heat", + "total agriculture machinery", + "total agriculture", + ] + agri_ele = [ + "total agriculture electricity", + ] + + road_total = [ + "total road", + "total passenger cars", + "total other road passenger", + "total light duty road freight", + ] + road_ele = [ + "electricity road", + "electricity passenger cars", + "electricity other road passenger", + "electricity light duty road freight", + ] + + rail_total = [ + "total rail", + "total rail passenger", + "total rail freight", + ] + rail_ele = [ + "electricity rail", + "electricity rail passenger", + "electricity rail freight", + ] + + avia_inter = [ + 'total aviation passenger', + 'total aviation freight', + 'total international aviation passenger', + 'total international aviation freight', + 'total international aviation' + ] + avia_domestic = [ + 'total domestic aviation passenger', + 'total domestic aviation freight', + 'total domestic aviation', + ] + navigation = [ + "total domestic navigation", + ] + + for country in idees_countries: + res = ratio.loc[(country, slice(None), 'Residential')] + energy.loc[country, residential_total] *= res[['total']].iloc[0,0] + energy.loc[country, residential_ele] *= res[['ele']].iloc[0,0] + + ser = ratio.loc[(country, slice(None), 'Services')] + energy.loc[country, service_total] *= ser[['total']].iloc[0,0] + energy.loc[country, service_ele] *= ser[['ele']].iloc[0,0] + + agri = ratio.loc[(country, slice(None), 'Agriculture & forestry')] + energy.loc[country, agri_total] *= agri[['total']].iloc[0,0] + energy.loc[country, agri_ele] *= agri[['ele']].iloc[0,0] + + road = ratio.loc[(country, slice(None), 'Road')] + energy.loc[country, road_total] *= road[['total']].iloc[0,0] + energy.loc[country, road_ele] *= road[['ele']].iloc[0,0] + + rail = ratio.loc[(country, slice(None), 'Rail')] + energy.loc[country, rail_total] *= rail[['total']].iloc[0,0] + energy.loc[country, rail_ele] *= rail[['ele']].iloc[0,0] + + avi_d = ratio.loc[(country, slice(None), 'Domestic aviation')] + avi_i = ratio.loc[(country, 'International aviation', slice(None))] + energy.loc[country, avia_inter] *= avi_i[['total']].iloc[0,0] + energy.loc[country, avia_domestic] *= avi_d[['total']].iloc[0,0] + + nav = ratio.loc[(country, slice(None), 'Domestic Navigation')] + energy.loc[country, navigation] *= nav[['total']].iloc[0,0] + + return energy if __name__ == "__main__": if "snakemake" not in globals(): @@ -758,12 +938,22 @@ def build_transport_data(countries, population, idees): input_eurostat = snakemake.input.eurostat eurostat = build_eurostat(input_eurostat, countries, report_year, data_year) swiss = build_swiss(data_year) - idees = build_idees(idees_countries, data_year) + # data from idees only exists for 2015 + if data_year > 2015: + # read in latest data and rescale later + idees = build_idees(idees_countries, 2015) + else: + idees = build_idees(idees_countries, data_year) energy = build_energy_totals(countries, eurostat, swiss, idees) + + if data_year > 2015: + energy = rescale(idees_countries, energy, eurostat) + energy.to_csv(snakemake.output.energy_name) - district_heat_share = build_district_heat_share(countries, idees) + # use rescaled idees data to calculate district heat share + district_heat_share = build_district_heat_share(countries, energy.loc[idees_countries]) district_heat_share.to_csv(snakemake.output.district_heat_share) base_year_emissions = params["base_emissions_year"] diff --git a/scripts/build_existing_heating_distribution.py b/scripts/build_existing_heating_distribution.py index 78518597f..a5e54a26e 100644 --- a/scripts/build_existing_heating_distribution.py +++ b/scripts/build_existing_heating_distribution.py @@ -14,7 +14,13 @@ def build_existing_heating(): - # retrieve existing heating capacities + # Add existing heating capacities, data comes from the study + # "Mapping and analyses of the current and future (2020 - 2030) + # heating/cooling fuel deployment (fossil/renewables) " + # https://energy.ec.europa.eu/publications/mapping-and-analyses-current-and-future-2020-2030-heatingcooling-fuel-deployment-fossilrenewables-1_en + # file: "WP2_DataAnnex_1_BuildingTechs_ForPublication_201603.xls" -> "existing_heating_raw.csv". + # data is for buildings only (i.e. NOT district heating) and represents the year 2012 + # TODO start from original file # Add existing heating capacities, data comes from the study # "Mapping and analyses of the current and future (2020 - 2030) diff --git a/scripts/build_industrial_energy_demand_per_country_today.py b/scripts/build_industrial_energy_demand_per_country_today.py index d1c672f19..65569b555 100644 --- a/scripts/build_industrial_energy_demand_per_country_today.py +++ b/scripts/build_industrial_energy_demand_per_country_today.py @@ -73,7 +73,7 @@ def industrial_energy_demand_per_country(country, year, jrc_dir): def get_subsector_data(sheet): df = df_dict[sheet][year].groupby(fuels).sum() - df["ammonia"] = 0.0 + df["hydrogen"] = 0.0 df["other"] = df["all"] - df.loc[df.index != "all"].sum() @@ -94,51 +94,34 @@ def get_subsector_data(sheet): return df -def add_ammonia_energy_demand(demand): - # MtNH3/a - fn = snakemake.input.ammonia_production - ammonia = pd.read_csv(fn, index_col=0)[str(year)] / 1e3 +def separate_basic_chemicals(demand, production): - def get_ammonia_by_fuel(x): - fuels = { - "gas": params["MWh_CH4_per_tNH3_SMR"], - "electricity": params["MWh_elec_per_tNH3_SMR"], - } - - return pd.Series({k: x * v for k, v in fuels.items()}) - - ammonia_by_fuel = ammonia.apply(get_ammonia_by_fuel).T - ammonia_by_fuel = ammonia_by_fuel.unstack().reindex( - index=demand.index, fill_value=0.0 - ) - - ammonia = pd.DataFrame({"ammonia": ammonia * params["MWh_NH3_per_tNH3"]}).T + ammonia = pd.DataFrame({"hydrogen": production["Ammonia"] * params["MWh_H2_per_tNH3_electrolysis"], + "electricity" : production["Ammonia"] * params["MWh_elec_per_tNH3_electrolysis"]}).T + chlorine = pd.DataFrame({"hydrogen": production["Chlorine"] * params["MWh_H2_per_tCl"], + "electricity" : production["Chlorine"] * params["MWh_elec_per_tCl"]}).T + methanol = pd.DataFrame({"gas": production["Methanol"] * params["MWh_CH4_per_tMeOH"], + "electricity" : production["Methanol"] * params["MWh_elec_per_tMeOH"]}).T demand["Ammonia"] = ammonia.unstack().reindex(index=demand.index, fill_value=0.0) + demand["Chlorine"] = chlorine.unstack().reindex(index=demand.index, fill_value=0.0) + demand["Methanol"] = methanol.unstack().reindex(index=demand.index, fill_value=0.0) - demand["Basic chemicals (without ammonia)"] = ( - demand["Basic chemicals"] - ammonia_by_fuel + demand["HVC"] = ( + demand["Basic chemicals"] - demand["Ammonia"] - demand["Methanol"] - demand["Chlorine"] ) - demand["Basic chemicals (without ammonia)"].clip(lower=0, inplace=True) - demand.drop(columns="Basic chemicals", inplace=True) + demand["HVC"].clip(lower=0, inplace=True) + return demand -def add_non_eu28_industrial_energy_demand(countries, demand): +def add_non_eu28_industrial_energy_demand(countries, demand, production): non_eu28 = countries.difference(eu28) if non_eu28.empty: return demand - # output in MtMaterial/a - fn = snakemake.input.industrial_production_per_country - production = pd.read_csv(fn, index_col=0) / 1e3 - - # recombine HVC, Chlorine and Methanol to Basic chemicals (without ammonia) - chemicals = ["HVC", "Chlorine", "Methanol"] - production["Basic chemicals (without ammonia)"] = production[chemicals].sum(axis=1) - production.drop(columns=chemicals, inplace=True) eu28_production = production.loc[countries.intersection(eu28)].sum() eu28_energy = demand.groupby(level=1).sum() @@ -182,9 +165,13 @@ def industrial_energy_demand(countries, year): demand = industrial_energy_demand(countries.intersection(eu28), year) - demand = add_ammonia_energy_demand(demand) + # output in MtMaterial/a + production = pd.read_csv(snakemake.input.industrial_production_per_country, + index_col=0) / 1e3 + + demand = separate_basic_chemicals(demand, production) - demand = add_non_eu28_industrial_energy_demand(countries, demand) + demand = add_non_eu28_industrial_energy_demand(countries, demand, production) # for format compatibility demand = demand.stack(dropna=False).unstack(level=[0, 2]) diff --git a/scripts/build_industrial_energy_demand_per_node.py b/scripts/build_industrial_energy_demand_per_node.py index 55c10c5d6..84f8679ad 100644 --- a/scripts/build_industrial_energy_demand_per_node.py +++ b/scripts/build_industrial_energy_demand_per_node.py @@ -19,23 +19,28 @@ planning_horizons=2030, ) - # import EU ratios df as csv + # import ratios fn = snakemake.input.industry_sector_ratios - industry_sector_ratios = pd.read_csv(fn, index_col=0) + sector_ratios = pd.read_csv(fn, + header=[0,1], + index_col=0) - # material demand per node and industry (kton/a) + # material demand per node and industry (Mton/a) fn = snakemake.input.industrial_production_per_node - nodal_production = pd.read_csv(fn, index_col=0) + nodal_production = pd.read_csv(fn, index_col=0) / 1e3 # energy demand today to get current electricity fn = snakemake.input.industrial_energy_demand_per_node_today nodal_today = pd.read_csv(fn, index_col=0) - # final energy consumption per node and industry (TWh/a) - nodal_df = nodal_production.dot(industry_sector_ratios.T) + nodal_sector_ratios = pd.concat({node: sector_ratios[node[:2]] for node in nodal_production.index}, + axis=1) + + nodal_production_stacked = nodal_production.stack() + nodal_production_stacked.index.names = [None,None] - # convert GWh to TWh and ktCO2 to MtCO2 - nodal_df *= 0.001 + # final energy consumption per node and industry (TWh/a) + nodal_df = (nodal_sector_ratios.multiply(nodal_production_stacked)).T.groupby(level=0).sum() rename_sectors = { "elec": "electricity", diff --git a/scripts/build_industrial_production_per_country.py b/scripts/build_industrial_production_per_country.py index 0aea4f150..afbe4de72 100644 --- a/scripts/build_industrial_production_per_country.py +++ b/scripts/build_industrial_production_per_country.py @@ -261,7 +261,7 @@ def separate_basic_chemicals(demand, year): demand["Basic chemicals"].clip(lower=0.0, inplace=True) # assume HVC, methanol, chlorine production proportional to non-ammonia basic chemicals - distribution_key = demand["Basic chemicals"] / demand["Basic chemicals"].sum() + distribution_key = demand["Basic chemicals"] / params["basic_chemicals_without_NH3_production_today"] / 1e3 demand["HVC"] = params["HVC_production_today"] * 1e3 * distribution_key demand["Chlorine"] = params["chlorine_production_today"] * 1e3 * distribution_key demand["Methanol"] = params["methanol_production_today"] * 1e3 * distribution_key diff --git a/scripts/build_industry_sector_ratios.py b/scripts/build_industry_sector_ratios.py index 457050026..805cf3e71 100644 --- a/scripts/build_industry_sector_ratios.py +++ b/scripts/build_industry_sector_ratios.py @@ -408,15 +408,15 @@ def chemicals_industry(): df.loc["methane", sector] -= ammonia_total * params["MWh_CH4_per_tNH3_SMR"] df.loc["elec", sector] -= ammonia_total * params["MWh_elec_per_tNH3_SMR"] - # subtract chlorine demand + # subtract chlorine demand (in MtCl/a) chlorine_total = params["chlorine_production_today"] - df.loc["hydrogen", sector] -= chlorine_total * params["MWh_H2_per_tCl"] - df.loc["elec", sector] -= chlorine_total * params["MWh_elec_per_tCl"] + df.loc["hydrogen", sector] -= chlorine_total * params["MWh_H2_per_tCl"] * 1e3 + df.loc["elec", sector] -= chlorine_total * params["MWh_elec_per_tCl"] * 1e3 - # subtract methanol demand + # subtract methanol demand (in MtMeOH/a) methanol_total = params["methanol_production_today"] - df.loc["methane", sector] -= methanol_total * params["MWh_CH4_per_tMeOH"] - df.loc["elec", sector] -= methanol_total * params["MWh_elec_per_tMeOH"] + df.loc["methane", sector] -= methanol_total * params["MWh_CH4_per_tMeOH"] * 1e3 + df.loc["elec", sector] -= methanol_total * params["MWh_elec_per_tMeOH"] * 1e3 # MWh/t material df.loc[sources, sector] = df.loc[sources, sector] / s_out diff --git a/scripts/build_industry_sector_ratios_intermediate.py b/scripts/build_industry_sector_ratios_intermediate.py new file mode 100644 index 000000000..86f88218a --- /dev/null +++ b/scripts/build_industry_sector_ratios_intermediate.py @@ -0,0 +1,77 @@ +# -*- coding: utf-8 -*- +# SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors +# +# SPDX-License-Identifier: MIT +""" +Build specific energy consumption by carrier and industries and by country, +that interpolates between the current average energy consumption (from 2015-2020) +and the ideal future best-in-class consumption. +""" + +import pandas as pd + +from prepare_sector_network import get + +def build_industry_sector_ratios_intermediate(): + + # in TWh/a + demand = pd.read_csv(snakemake.input.industrial_energy_demand_per_country_today, + header=[0,1], + index_col=0) + + # in Mt/a + production = pd.read_csv(snakemake.input.industrial_production_per_country, + index_col=0) / 1e3 + production = production.unstack().swaplevel() + + # in MWh/t + future_sector_ratios = pd.read_csv(snakemake.input.industry_sector_ratios, + index_col=0) + + production.index.names = [None,None] + + today_sector_ratios = demand.div(production, axis=1) + + today_sector_ratios.drop(columns=today_sector_ratios.columns[today_sector_ratios.isna().all()], + inplace=True) + + rename = pd.Series(today_sector_ratios.index, + today_sector_ratios.index) + rename["waste"] = "biomass" + rename["electricity"] = "elec" + rename["solid"] = "coke" + rename["gas"] = "methane" + rename["other"] = "biomass" + rename["liquid"] = "naphtha" + + today_sector_ratios.rename(rename, + inplace=True) + + + fraction_future = get(params["sector_ratios_fraction_future"], year) + + intermediate_sector_ratios = {} + + for ct in today_sector_ratios.columns.unique(level=0): + + intermediate_sector_ratio = future_sector_ratios.copy() + + intermediate_sector_ratio.loc[today_sector_ratios[ct].index,today_sector_ratios[ct].columns] = (fraction_future*intermediate_sector_ratio.loc[today_sector_ratios[ct].index,today_sector_ratios[ct].columns] + + (1 - fraction_future)*today_sector_ratios[ct]) + intermediate_sector_ratios[ct] = intermediate_sector_ratio + + intermediate_sector_ratios = pd.concat(intermediate_sector_ratios, axis=1) + + intermediate_sector_ratios.to_csv(snakemake.output.industry_sector_ratios) + +if __name__ == "__main__": + if "snakemake" not in globals(): + from _helpers import mock_snakemake + + snakemake = mock_snakemake("build_industry_sector_ratios_intermediate") + + year = int(snakemake.wildcards.planning_horizons[-4:]) + + params = snakemake.params.industry + + build_industry_sector_ratios_intermediate() diff --git a/scripts/plot_statistics_comparison.py b/scripts/plot_statistics_comparison.py new file mode 100644 index 000000000..ce20340b7 --- /dev/null +++ b/scripts/plot_statistics_comparison.py @@ -0,0 +1,108 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# SPDX-FileCopyrightText: : 2017-2023 The PyPSA-Eur Authors +# +# SPDX-License-Identifier: MIT + +import re + +import matplotlib.pyplot as plt +import pandas as pd +import seaborn as sns +from _helpers import configure_logging +from plot_summary import rename_techs + +sns.set_theme("paper", style="whitegrid") +STACKED = { + "capacity_factor": False, + "market_value": False, +} + + +def rename_index(df): + # rename index and drop duplicates + specific = df.index.map(lambda x: f"{x[1]}({x[0]})") + generic = df.index.get_level_values("carrier") + duplicated = generic.duplicated(keep=False) + index = specific.where(duplicated, generic) + df = df.set_axis(index) + # rename columns and drop duplicates + columns = df.columns.str.split("_", expand=True) + columns = [ + columns.get_level_values(level).unique() + for level in range(columns.nlevels) + if not columns.get_level_values(level).duplicated(keep=False).all() + ] + columns = pd.MultiIndex.from_arrays(columns) + df.columns = columns.map("\n".join) + return df + + +def plot_static_comparison(df, ax, stacked=False): + factor, unit = conversion[output] + df = df[df != 0] + df = df.dropna(axis=0, how="all").fillna(0) + if df.empty: + return + c = tech_colors[df.index.get_level_values("carrier").map(rename_techs)] + df = df.pipe(rename_index).T + df = df.div(float(factor)) if factor != "-" else df + df.plot.bar(color=c.values, ax=ax, stacked=stacked, legend=False, ylabel=unit) + ax.legend( + bbox_to_anchor=(1, 1), + loc="upper left", + ) + ax.grid(axis="x") + + +def read_csv(input, output): + try: + # filter required csv to plot the wanted output + files = list(filter(lambda x: output in x, input)) + pattern = r"elec_.*?(\d{4})" + network_labels = [re.search(pattern, f).group() for f in files] + df = pd.concat( + [ + pd.read_csv(f, skiprows=2).set_index(["component", "carrier"]) + for f in files + ], + axis=1, + keys=network_labels, + ) + # get plot label and drop from index + label = df.columns.get_level_values(1).unique()[0] + df.columns = df.columns.droplevel(1) + except Exception as e: + print(f"Error reading csv file for {output}: {e}") + df = pd.DataFrame() + return df + + +if __name__ == "__main__": + if "snakemake" not in globals(): + from _helpers import mock_snakemake + + snakemake = mock_snakemake( + "plot_statistics_comparison", + country="all", + carrier="electricity", + ) + configure_logging(snakemake) + + tech_colors = pd.Series(snakemake.params.plotting["tech_colors"]) + conversion = pd.Series(snakemake.params.statistics) + + for output in snakemake.output.keys(): + if "touch" in output: + with open(snakemake.output.barplots_touch, "a"): + pass + continue + fig, ax = plt.subplots() + df = read_csv(snakemake.input, output) + if df.empty: + fig.savefig(snakemake.output[output]) + continue + + plot_static_comparison(df, ax, stacked=STACKED.get(output, True)) + + fig.savefig(snakemake.output[output], bbox_inches="tight") diff --git a/scripts/plot_statistics_single.py b/scripts/plot_statistics_single.py new file mode 100644 index 000000000..97491f58e --- /dev/null +++ b/scripts/plot_statistics_single.py @@ -0,0 +1,75 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# SPDX-FileCopyrightText: : 2017-2023 The PyPSA-Eur Authors +# +# SPDX-License-Identifier: MIT + +import matplotlib.pyplot as plt +import pandas as pd +import seaborn as sns +from _helpers import configure_logging +from plot_summary import rename_techs + +sns.set_theme("paper", style="whitegrid") + + +def rename_index(ds): + specific = ds.index.map(lambda x: f"{x[1]}\n({x[0]})") + generic = ds.index.get_level_values("carrier") + duplicated = generic.duplicated(keep=False) + index = specific.where(duplicated, generic) + return ds.set_axis(index) + + +def plot_static_single(ds, ax): + factor, unit = conversion[output] + ds = ds.dropna() + c = tech_colors[ds.index.get_level_values("carrier").map(rename_techs)] + ds = ds.pipe(rename_index) + ds = ds.div(float(factor)) if factor != "-" else ds + ds.T.plot.barh(color=c.values, ax=ax, xlabel=unit) + ax.grid(axis="y") + + +def read_csv(input): + try: + df = pd.read_csv(input, skiprows=2) + df = df.set_index(["component", "carrier"]).iloc[:, 0] + except Exception as e: + print(f"An error occurred reading file {input}: {e}") + df = pd.DataFrame() + return df + + +if __name__ == "__main__": + if "snakemake" not in globals(): + from _helpers import mock_snakemake + + snakemake = mock_snakemake( + "plot_statistics_single", + simpl="", + ll="v1.5", + clusters="5", + opts="", + sector_opts="24h-T-H-B-I-A-dist1", + planning_horizons="2040", + country="all", + carrier="heat", + ) + configure_logging(snakemake) + + tech_colors = pd.Series(snakemake.params.plotting["tech_colors"]) + conversion = pd.Series(snakemake.params.statistics) + + for output in snakemake.output.keys(): + if "touch" in output: + with open(snakemake.output.barplots_touch, "a"): + pass + continue + fig, ax = plt.subplots() + ds = read_csv(snakemake.input[output]) + if ds.empty: + fig.savefig(snakemake.output[output]) + continue + plot_static_single(ds, ax) + fig.savefig(snakemake.output[output], bbox_inches="tight") diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index be8aea532..049f23181 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -172,6 +172,13 @@ def define_spatial(nodes, options): spatial.coal.nodes = ["EU coal"] spatial.coal.locations = ["EU"] + if options["regional_coal_demand"]: + spatial.coal.demand_locations = nodes + spatial.coal.industry = nodes + " coal for industry" + else: + spatial.coal.demand_locations = ["EU"] + spatial.coal.industry = ["EU coal for industry"] + # lignite spatial.lignite = SimpleNamespace() spatial.lignite.nodes = ["EU lignite"] @@ -3048,19 +3055,42 @@ def add_industry(n, costs): mwh_coal_per_mwh_coke = 1.366 # from eurostat energy balance p_set = ( - industrial_demand["coal"].sum() - + mwh_coal_per_mwh_coke * industrial_demand["coke"].sum() + industrial_demand["coal"] + + mwh_coal_per_mwh_coke * industrial_demand["coke"] ) / nhours + p_set.rename(lambda x: x + " coal for industry", inplace=True) + + if not options["regional_coal_demand"]: + p_set = p_set.sum() + + n.madd( + "Bus", + spatial.coal.industry, + location=spatial.coal.demand_locations, + carrier="coal for industry", + unit="MWh_LHV", + ) + n.madd( "Load", - spatial.coal.nodes, - suffix=" for industry", - bus=spatial.coal.nodes, + spatial.coal.industry, + bus=spatial.coal.industry, carrier="coal for industry", p_set=p_set, ) + n.madd( + "Link", + spatial.coal.industry, + bus0=spatial.coal.nodes, + bus1=spatial.coal.industry, + bus2="co2 atmosphere", + carrier="coal for industry", + p_nom_extendable=True, + efficiency2=costs.at["coal", "CO2 intensity"], + ) + def add_waste_heat(n): # TODO options? @@ -3378,7 +3408,11 @@ def define_clustering(attributes, aggregate_dict): for c in n.iterate_components(components): df = c.df - cols = df.columns[df.columns.str.contains("bus") | (df.columns == "carrier")] + cols = df.columns[ + df.columns.str.contains("bus") + | (df.columns == "carrier") + | (df.columns == "nice_name") + ] # rename columns and index df[cols] = df[cols].apply( @@ -3507,10 +3541,12 @@ def set_temporal_aggregation(n, opts, solver_name): return n -def lossy_bidirectional_links(n, carrier, efficiencies={}): +def lossy_bidirectional_links(n, carrier, efficiencies={}, subset=None): "Split bidirectional links into two unidirectional links to include transmission losses." - carrier_i = n.links.query("carrier == @carrier").index + if subset is None: + subset = n.links.index + carrier_i = n.links.query("carrier == @carrier").index.intersection(subset) if ( not any((v != 1.0) or (v >= 0) for v in efficiencies.values()) diff --git a/scripts/write_statistics.py b/scripts/write_statistics.py new file mode 100644 index 000000000..98d9a5c3a --- /dev/null +++ b/scripts/write_statistics.py @@ -0,0 +1,118 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# SPDX-FileCopyrightText: : 2017-2023 The PyPSA-Eur Authors +# +# SPDX-License-Identifier: MIT + +import matplotlib.pyplot as plt +import pandas as pd +import pypsa +import seaborn as sns +from _helpers import configure_logging +from pypsa.statistics import get_carrier + + +# grouperfunctions = hier schreiben und dann in statistics. +def groupby_country_and_carrier(n, c, nice_names=False): + df = n.df(c) + bus = "bus1" if "bus" not in n.df(c) else "bus" + country = df[bus].map(n.buses.location).map(n.buses.country).rename("country") + carrier = get_carrier(n, c, nice_names) + return [country, carrier] + + +def call_with_handle(func, **kwargs): + try: + ds = func(**kwargs) + except Exception as e: + print(f"An error occurred: {e}") + ds = pd.Series() + pass + return ds + + +if __name__ == "__main__": + if "snakemake" not in globals(): + from _helpers import mock_snakemake + + snakemake = mock_snakemake( + "save_statistics_csv", + simpl="", + ll="v1.5", + clusters="5", + opts="", + sector_opts="24h-T-H-B-I-A-dist1", + planning_horizons="2040", + country="all", + carrier="electricity", + ) + # configure_logging(snakemake) + config = snakemake.config + + n = pypsa.Network(snakemake.input.network) + + kwargs = {"nice_names": False} + + wildcards = dict(snakemake.wildcards) + carrier = wildcards.pop("carrier") + country = wildcards.pop("country") + + if carrier == "all": + pass + elif carrier in config["plotting"].get("carrier_groups", []): + bus_carrier = config["plotting"]["carrier_groups"][carrier] + kwargs["bus_carrier"] = bus_carrier + elif n.buses.carrier.str.contains(carrier).any(): + if carrier not in n.buses.carrier.unique(): + carrier = [ + bus_carrier + for bus_carrier in n.buses.carrier.unique() + if carrier in bus_carrier + ] + bus_carrier = carrier + kwargs["bus_carrier"] = bus_carrier + else: + raise ValueError( + f"Carrier {carrier} not found in network or carrier group in config." + ) + + if country == "all" or not country: + kwargs["groupby"] = get_carrier + else: + kwargs["groupby"] = groupby_country_and_carrier + + for output in snakemake.output.keys(): + if "touch" in output: + continue + if output == "energy_balance": + supply = call_with_handle(n.statistics.supply, **kwargs) + withdrawal = call_with_handle(n.statistics.withdrawal, **kwargs) + ds = ( + pd.concat([supply, withdrawal.mul(-1)]) + .groupby(supply.index.names) + .sum() + ) + ds.attrs = supply.attrs + ds.attrs["name"] = "Energy Balance" + elif output == "total_cost": + opex = call_with_handle(n.statistics.opex, **kwargs) + capex = call_with_handle(n.statistics.capex, **kwargs) + ds = opex.add(capex, fill_value=0) + ds.attrs["name"] = "Total Cost" + else: + func = eval(f"n.statistics.{output}") + ds = call_with_handle(func, **kwargs) + + if ds.empty: + print( + f"Empty series for {output} with bus carrier {bus_carrier} and country {country}." + ) + pd.Series().to_csv(snakemake.output[output]) + continue + if country and country != "all": + ds = ds.xs(country, level="country") + pd.Series(ds.attrs).to_csv(snakemake.output[output], header=False) + ds.dropna().to_csv(snakemake.output[output], mode="a") + # touch file + with open(snakemake.output.csv_touch, "a"): + pass