diff --git a/Snakefile b/Snakefile index 10fa6e04..37fd0f9f 100644 --- a/Snakefile +++ b/Snakefile @@ -3,9 +3,9 @@ configfile: "config.yaml" COSTS="data/costs.csv" wildcard_constraints: - lv="[0-9\.]+|inf", - simpl="[a-zA-Z0-9]*", - clusters="[0-9]+m?", + lv="[0-9\.]+|inf|all", + simpl="[a-zA-Z0-9]*|all", + clusters="[0-9]+m?|all", sectors="[+a-zA-Z0-9]+", opts="[-+a-zA-Z0-9]*" @@ -256,42 +256,23 @@ rule plot_network: ext="results/plots/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_{attr}_ext.{ext}" script: "scripts/plot_network.py" -# rule plot_costs: -# input: 'results/summaries/costs2-summary.csv' -# output: -# expand('results/plots/costs_{cost}_{resarea}_{sectors}_{opt}', -# **dict(chain(config['scenario'].items(), (('{param}'))) -# touch('results/plots/scenario_plots') -# params: -# tmpl="results/plots/costs_[cost]_[resarea]_[sectors]_[opt]" -# exts=["pdf", "png"] -# scripts: "scripts/plot_costs.py" +def summary_networks(w): + # It's mildly hacky to include the separate costs input as first entry + return ([COSTS] + + expand("results/networks/{network}_s{simpl}_{clusters}_lv{lv}_{opts}.nc", + network=w.network, + **{k: config["scenario"][k] if getattr(w, k) == "all" else getattr(w, k) + for k in ["simpl", "clusters", "lv", "opts"]})) -# rule scenario_comparison: -# input: -# expand('results/plots/network_{cost}_{sectors}_{opts}_{attr}.pdf', -# version=config['version'], -# attr=['p_nom'], -# **config['scenario']) -# output: -# html='results/plots/scenario_{param}.html' -# params: -# tmpl="network_[cost]_[resarea]_[sectors]_[opts]_[attr]", -# plot_dir='results/plots' -# script: "scripts/scenario_comparison.py" - -# rule extract_summaries: -# input: -# expand("results/networks/{cost}_{sectors}_{opts}.nc", -# **config['scenario']) -# output: -# **{n: "results/summaries/{}-summary.csv".format(n) -# for n in ['costs', 'costs2', 'e_curtailed', 'e_nom_opt', 'e', 'p_nom_opt']} -# params: -# scenario_tmpl="[cost]_[resarea]_[sectors]_[opts]", -# scenarios=config['scenario'] -# script: "scripts/extract_summaries.py" +rule make_summary: + input: summary_networks + output: directory("results/summaries/{network}_s{simpl}_{clusters}_lv{lv}_{opts}") + script: "scripts/make_summary.py" +rule plot_summary: + input: directory("results/summaries/{network}_s{simpl}_{clusters}_lv{lv}_{opts}") + output: "results/plots/summary_{summary}_{network}_s{simpl}_{clusters}_lv{lv}_{opts}.{ext}" + script: "scripts/plot_summary.py" # Local Variables: # mode: python diff --git a/config.yaml b/config.yaml index e0c9f7b6..82988f13 100644 --- a/config.yaml +++ b/config.yaml @@ -205,6 +205,7 @@ plotting: 'offwind' : "c" 'offshore wind' : "c" "hydro" : "#3B5323" + "hydro+PHS" : "#3B5323" "hydro reservoir" : "#3B5323" "ror" : "#78AB46" "run of river" : "#78AB46" @@ -267,6 +268,9 @@ plotting: "transport" : "grey" "electricity" : "k" "transport fuel cell" : "#AAAAAA" + # _helpers.load_network requirements + "AC-AC" : "k" + "AC line" : "k" nice_names: # OCGT: "Gas" # OCGT marginal: "Gas (marginal)" diff --git a/scripts/make_summary.py b/scripts/make_summary.py index 50e0b3d3..d9d55220 100644 --- a/scripts/make_summary.py +++ b/scripts/make_summary.py @@ -1,30 +1,16 @@ - +import os from six import iteritems - -import sys - -sys.path = ['/home/vres/lib/python3.5/site-packages'] + sys.path - +from itertools import product import pandas as pd import pypsa -from vresutils.costdata import annuity - -from prepare_network import generate_periodic_profiles - -from add_electricity import load_costs - -import yaml +from add_electricity import load_costs, update_transmission_costs idx = pd.IndexSlice opt_name = {"Store": "e", "Line" : "s", "Transformer" : "s"} - -#separator to find group name -find_by = " " - def assign_carriers(n): if "carrier" not in n.loads: @@ -35,7 +21,7 @@ def assign_carriers(n): if "carrier" not in n.lines: n.lines["carrier"] = "AC" - if n.stores.loc["EU gas Store","carrier"] == "": + if "EU gas store" in n.stores.index and n.stores.loc["EU gas Store","carrier"] == "": n.stores.loc["EU gas Store","carrier"] = "gas Store" @@ -68,17 +54,6 @@ def calculate_costs(n,label,costs): costs.loc[idx[c.list_name,"marginal",list(marginal_costs_grouped.index)],label] = marginal_costs_grouped.values - #add back in costs of links if there is a line volume limit - if label[1] != "opt": - costs.loc[("links-added","capital","transmission lines"),label] = ((costs_db.at['HVDC overhead', 'capital_cost']*n.links.length + costs_db.at['HVDC inverter pair', 'capital_cost'])*n.links.p_nom_opt)[n.links.carrier == "DC"].sum() - costs.loc[("lines-added","capital","transmission lines"),label] = costs_db.at["HVAC overhead", "capital_cost"]*(n.lines.length*n.lines.s_nom_opt).sum() - - - #add back in all hydro - #costs.loc[("storage_units","capital","hydro"),label] = (0.01)*2e6*n.storage_units.loc[n.storage_units.group=="hydro","p_nom"].sum() - #costs.loc[("storage_units","capital","PHS"),label] = (0.01)*2e6*n.storage_units.loc[n.storage_units.group=="PHS","p_nom"].sum() - #costs.loc[("generators","capital","ror"),label] = (0.02)*3e6*n.generators.loc[n.generators.group=="ror","p_nom"].sum() - return costs @@ -203,6 +178,8 @@ def calculate_metrics(n,label,metrics): if hasattr(n,"line_volume_limit"): metrics.at["line_volume_limit",label] = n.line_volume_limit + + if hasattr(n,"line_volume_limit_dual"): metrics.at["line_volume_shadow",label] = n.line_volume_limit_dual if "CO2Limit" in n.global_constraints.index: @@ -284,83 +261,59 @@ def calculate_weighted_prices(n,label,weighted_prices): +# BROKEN don't use +# +# def calculate_market_values(n, label, market_values): +# # Warning: doesn't include storage units -def calculate_market_values(n, label, market_values): - # Warning: doesn't include storage units +# n.buses["suffix"] = n.buses.index.str[2:] +# suffix = "" +# buses = n.buses.index[n.buses.suffix == suffix] - n.buses["suffix"] = n.buses.index.str[2:] +# ## First do market value of generators ## +# generators = n.generators.index[n.buses.loc[n.generators.bus,"suffix"] == suffix] +# techs = n.generators.loc[generators,"carrier"].value_counts().index +# market_values = market_values.reindex(market_values.index | techs) - suffix = "" +# for tech in techs: +# gens = generators[n.generators.loc[generators,"carrier"] == tech] +# dispatch = n.generators_t.p[gens].groupby(n.generators.loc[gens,"bus"],axis=1).sum().reindex(columns=buses,fill_value=0.) +# revenue = dispatch*n.buses_t.marginal_price[buses] +# market_values.at[tech,label] = revenue.sum().sum()/dispatch.sum().sum() - buses = n.buses.index[n.buses.suffix == suffix] +# ## Now do market value of links ## + +# for i in ["0","1"]: +# all_links = n.links.index[n.buses.loc[n.links["bus"+i],"suffix"] == suffix] +# techs = n.links.loc[all_links,"carrier"].value_counts().index +# market_values = market_values.reindex(market_values.index | techs) + +# for tech in techs: +# links = all_links[n.links.loc[all_links,"carrier"] == tech] +# dispatch = n.links_t["p"+i][links].groupby(n.links.loc[links,"bus"+i],axis=1).sum().reindex(columns=buses,fill_value=0.) +# revenue = dispatch*n.buses_t.marginal_price[buses] +# market_values.at[tech,label] = revenue.sum().sum()/dispatch.sum().sum() + +# return market_values - ## First do market value of generators ## +# OLD CODE must be adapted - generators = n.generators.index[n.buses.loc[n.generators.bus,"suffix"] == suffix] - - techs = n.generators.loc[generators,"carrier"].value_counts().index - - market_values = market_values.reindex(market_values.index | techs) +# def calculate_price_statistics(n, label, price_statistics): - for tech in techs: - gens = generators[n.generators.loc[generators,"carrier"] == tech] +# price_statistics = price_statistics.reindex(price_statistics.index|pd.Index(["zero_hours","mean","standard_deviation"])) +# n.buses["suffix"] = n.buses.index.str[2:] +# suffix = "" +# buses = n.buses.index[n.buses.suffix == suffix] - dispatch = n.generators_t.p[gens].groupby(n.generators.loc[gens,"bus"],axis=1).sum().reindex(columns=buses,fill_value=0.) - - revenue = dispatch*n.buses_t.marginal_price[buses] - - market_values.at[tech,label] = revenue.sum().sum()/dispatch.sum().sum() - - - - ## Now do market value of links ## - - for i in ["0","1"]: - all_links = n.links.index[n.buses.loc[n.links["bus"+i],"suffix"] == suffix] - - techs = n.links.loc[all_links,"carrier"].value_counts().index - - market_values = market_values.reindex(market_values.index | techs) - - for tech in techs: - links = all_links[n.links.loc[all_links,"carrier"] == tech] - - dispatch = n.links_t["p"+i][links].groupby(n.links.loc[links,"bus"+i],axis=1).sum().reindex(columns=buses,fill_value=0.) - - revenue = dispatch*n.buses_t.marginal_price[buses] - - market_values.at[tech,label] = revenue.sum().sum()/dispatch.sum().sum() - - return market_values - - -def calculate_price_statistics(n, label, price_statistics): - - - price_statistics = price_statistics.reindex(price_statistics.index|pd.Index(["zero_hours","mean","standard_deviation"])) - - n.buses["suffix"] = n.buses.index.str[2:] - - suffix = "" - - buses = n.buses.index[n.buses.suffix == suffix] - - - threshold = 0.1 #higher than phoney marginal_cost of wind/solar - - df = pd.DataFrame(data=0.,columns=buses,index=n.snapshots) - - df[n.buses_t.marginal_price[buses] < threshold] = 1. - - price_statistics.at["zero_hours", label] = df.sum().sum()/(df.shape[0]*df.shape[1]) - - price_statistics.at["mean", label] = n.buses_t.marginal_price[buses].unstack().mean() - - price_statistics.at["standard_deviation", label] = n.buses_t.marginal_price[buses].unstack().std() - - return price_statistics +# threshold = 0.1 #higher than phoney marginal_cost of wind/solar +# df = pd.DataFrame(data=0.,columns=buses,index=n.snapshots) +# df[n.buses_t.marginal_price[buses] < threshold] = 1. +# price_statistics.at["zero_hours", label] = df.sum().sum()/(df.shape[0]*df.shape[1]) +# price_statistics.at["mean", label] = n.buses_t.marginal_price[buses].unstack().mean() +# price_statistics.at["standard_deviation", label] = n.buses_t.marginal_price[buses].unstack().std() +# return price_statistics outputs = ["costs", @@ -370,69 +323,66 @@ outputs = ["costs", "supply_energy", "prices", "weighted_prices", - "price_statistics", - "market_values", + # "price_statistics", + # "market_values", "metrics", ] def make_summaries(networks_dict): - columns = pd.MultiIndex.from_tuples(networks_dict.keys(),names=["cluster","lv","opt"]) + columns = pd.MultiIndex.from_tuples(networks_dict.keys(),names=["simpl","clusters","lv","opts"]) - df = {} + dfs = {} for output in outputs: - df[output] = pd.DataFrame(columns=columns,dtype=float) + dfs[output] = pd.DataFrame(columns=columns,dtype=float) for label, filename in iteritems(networks_dict): print(label, filename) + if not os.path.exists(filename): + print("does not exist!!") + continue n = pypsa.Network(filename) assign_carriers(n) + Nyears = n.snapshot_weightings.sum()/8760. + costs = load_costs(Nyears, snakemake.input[0], + snakemake.config['costs'], snakemake.config['electricity']) + update_transmission_costs(n, costs) + for output in outputs: - df[output] = globals()["calculate_" + output](n, label, df[output]) + dfs[output] = globals()["calculate_" + output](n, label, dfs[output]) - return df + return dfs -def to_csv(df): - - for key in df: - df[key].to_csv(snakemake.output[key]) +def to_csv(dfs): + dir = snakemake.output[0] + os.makedirs(dir, exist_ok=True) + for key, df in iteritems(dfs): + df.to_csv(os.path.join(dir, f"{key}.csv")) if __name__ == "__main__": # Detect running outside of snakemake and mock snakemake for testing - if 'snakemake' not in globals(): - from vresutils import Dict - import yaml - snakemake = Dict() - with open('config.yaml') as f: - snakemake.config = yaml.load(f) - snakemake.input = Dict() - snakemake.input['heat_demand_name'] = 'data/heating/daily_heat_demand.h5' - snakemake.output = Dict() - name = "37-lv" - clusters = [37] - lvs = [1.0,1.25,2.0,3.0] - opts = ["Co2L-3H-T-H"] - for item in outputs: - snakemake.output[item] = snakemake.config['summary_dir'] + '/{name}/csvs/{item}.csv'.format(name=name,item=item) + def expand_from_wildcard(key): + w = getattr(snakemake.wildcards, key) + return snakemake.config["scenario"][key] if w == "all" else [w] + + networks_dict = {(simpl,clusters,lv,opts) : ('results/networks/elec_s{simpl}_{clusters}_lv{lv}_{opts}.nc' + .format(simpl=simpl, + clusters=clusters, + opts=opts, + lv=lv)) + for simpl in expand_from_wildcard("simpl") + for clusters in expand_from_wildcard("clusters") + for lv in expand_from_wildcard("lv") + for opts in expand_from_wildcard("opts")} - networks_dict = {(cluster,lv,opt) : - 'results/networks/elec_s_{cluster}_lv{lv}_{opt}.nc'\ - .format(cluster=cluster, - opt=opt, - lv=lv)\ - for cluster in clusters \ - for opt in opts \ - for lv in lvs} print(networks_dict) - costs_db = load_costs(Nyears=1.,tech_costs="data/costs.csv",config=snakemake.config["costs"],elec_config=snakemake.config['electricity']) + dfs = make_summaries(networks_dict) - df = make_summaries(networks_dict) - - to_csv(df) + to_csv(dfs) diff --git a/scripts/plot_summary.py b/scripts/plot_summary.py index 83461f15..de2f0c24 100644 --- a/scripts/plot_summary.py +++ b/scripts/plot_summary.py @@ -1,56 +1,47 @@ - - - +import os import pandas as pd - -#allow plotting without Xwindows -import matplotlib -matplotlib.use('Agg') - import matplotlib.pyplot as plt - - #consolidate and rename def rename_techs(label): - if label[:8] == "central ": - label = label[8:] - if label[:6] == "urban ": - label = label[6:] + if label.startswith("central "): + label = label[len("central "):] + elif label.startswith("urban "): + label = label[len("urban "):] if "retrofitting" in label: label = "building retrofitting" - if "H2" in label: + elif "H2" in label: label = "hydrogen storage" - if "CHP" in label: + elif "CHP" in label: label = "CHP" - if "water tank" in label: + elif "water tank" in label: label = "water tanks" - if label=="water tanks": + elif label == "water tanks": label = "hot water storage" - if "gas" in label and label != "gas boiler": + elif "gas" in label and label != "gas boiler": label = "natural gas" - if "solar thermal" in label: + elif "solar thermal" in label: label = "solar thermal" - if label == "solar": + elif label == "solar": label = "solar PV" - if label == "heat pump": + elif label == "heat pump": label = "air heat pump" - if label == "Sabatier": + elif label == "Sabatier": label = "methanation" - if label == "offwind": + elif label == "offwind": label = "offshore wind" - if label == "onwind": + elif label == "onwind": label = "onshore wind" - if label == "ror": + elif label == "ror": label = "hydroelectricity" - if label == "hydro": + elif label == "hydro": label = "hydroelectricity" - if label == "PHS": + elif label == "PHS": label = "hydroelectricity" - if label == "co2 Store": + elif label == "co2 Store": label = "DAC" - if "battery" in label: + elif "battery" in label: label = "battery storage" return label @@ -58,11 +49,10 @@ def rename_techs(label): preferred_order = pd.Index(["transmission lines","hydroelectricity","hydro reservoir","run of river","pumped hydro storage","onshore wind","offshore wind","solar PV","solar thermal","building retrofitting","ground heat pump","air heat pump","resistive heater","CHP","OCGT","gas boiler","gas","natural gas","methanation","hydrogen storage","battery storage","hot water storage"]) -def plot_costs(): - - - cost_df = pd.read_csv(snakemake.input.costs,index_col=list(range(3)),header=[0,1,2]) +def plot_costs(infn, fn=None): + ## For now ignore the simpl header + cost_df = pd.read_csv(infn,index_col=list(range(3)),header=[1,2,3]) df = cost_df.groupby(cost_df.index.get_level_values(2)).sum() @@ -109,12 +99,13 @@ def plot_costs(): fig.tight_layout() - fig.savefig(snakemake.output.costs,transparent=True) + if fn is not None: + fig.savefig(fn, transparent=True) -def plot_energy(): +def plot_energy(infn, fn=None): - energy_df = pd.read_csv(snakemake.input.energy,index_col=list(range(2)),header=[0,1,2]) + energy_df = pd.read_csv(infn, index_col=list(range(2)),header=[1,2,3]) df = energy_df.groupby(energy_df.index.get_level_values(1)).sum() @@ -161,25 +152,15 @@ def plot_energy(): fig.tight_layout() - fig.savefig(snakemake.output.energy,transparent=True) + if fn is not None: + fig.savefig(fn, transparent=True) if __name__ == "__main__": - # Detect running outside of snakemake and mock snakemake for testing - if 'snakemake' not in globals(): - from vresutils import Dict - import yaml - snakemake = Dict() - with open('config.yaml') as f: - snakemake.config = yaml.load(f) - snakemake.input = Dict() - snakemake.output = Dict() - name = "37-lv" + summary = snakemake.wildcards.summary + try: + func = globals()[f"plot_{summary}"] + except KeyError: + raise RuntimeError(f"plotting function for {summary} has not been defined") - for item in ["costs","energy"]: - snakemake.input[item] = snakemake.config['summary_dir'] + '/{name}/csvs/{item}.csv'.format(name=name,item=item) - snakemake.output[item] = snakemake.config['summary_dir'] + '/{name}/graphs/{item}.pdf'.format(name=name,item=item) - - plot_costs() - - plot_energy() + func(os.path.join(snakemake.input[0], f"{summary}.csv"), snakemake.output[0])