pypsa-eur/scripts/paper_graphics.py

403 lines
13 KiB
Python

import pandas as pd
#allow plotting without Xwindows
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import pypsa
import numpy as np
from plot_summary import rename_techs, preferred_order
from make_summary import assign_carriers
#from sector/scripts/paper_graphics-co2_sweep.py
from matplotlib.patches import Circle, Ellipse
from matplotlib.legend_handler import HandlerPatch
import cartopy.crs as ccrs
override_component_attrs = pypsa.descriptors.Dict({k : v.copy() for k,v in pypsa.components.component_attrs.items()})
override_component_attrs["Link"].loc["bus2"] = ["string",np.nan,np.nan,"2nd bus","Input (optional)"]
override_component_attrs["Link"].loc["bus3"] = ["string",np.nan,np.nan,"3rd bus","Input (optional)"]
override_component_attrs["Link"].loc["efficiency2"] = ["static or series","per unit",1.,"2nd bus efficiency","Input (optional)"]
override_component_attrs["Link"].loc["efficiency3"] = ["static or series","per unit",1.,"3rd bus efficiency","Input (optional)"]
override_component_attrs["Link"].loc["p2"] = ["series","MW",0.,"2nd bus output","Output"]
override_component_attrs["Link"].loc["p3"] = ["series","MW",0.,"3rd bus output","Output"]
def rename_techs_tyndp(tech):
tech = rename_techs(tech)
if "heat pump" in tech or "resistive heater" in tech:
return "power-to-heat"
elif tech in ["methanation","hydrogen storage","helmeth"]:
return "power-to-gas"
elif tech in ["OCGT","CHP","gas boiler"]:
return "gas-to-power/heat"
elif "solar" in tech:
return "solar"
elif tech == "Fischer-Tropsch":
return "power-to-liquid"
elif "offshore wind" in tech:
return "offshore wind"
else:
return tech
def make_handler_map_to_scale_circles_as_in(ax, dont_resize_actively=False):
fig = ax.get_figure()
def axes2pt():
return np.diff(ax.transData.transform([(0,0), (1,1)]), axis=0)[0] * (72./fig.dpi)
ellipses = []
if not dont_resize_actively:
def update_width_height(event):
dist = axes2pt()
for e, radius in ellipses: e.width, e.height = 2. * radius * dist
fig.canvas.mpl_connect('resize_event', update_width_height)
ax.callbacks.connect('xlim_changed', update_width_height)
ax.callbacks.connect('ylim_changed', update_width_height)
def legend_circle_handler(legend, orig_handle, xdescent, ydescent,
width, height, fontsize):
w, h = 2. * orig_handle.get_radius() * axes2pt()
e = Ellipse(xy=(0.5*width-0.5*xdescent, 0.5*height-0.5*ydescent), width=w, height=w)
ellipses.append((e, orig_handle.get_radius()))
return e
return {Circle: HandlerPatch(patch_func=legend_circle_handler)}
def make_legend_circles_for(sizes, scale=1.0, **kw):
return [Circle((0,0), radius=(s/scale)**0.5, **kw) for s in sizes]
def assign_location(n):
for c in n.iterate_components(n.one_port_components|n.branch_components):
ifind = pd.Series(c.df.index.str.find(" ",start=4),c.df.index)
for i in ifind.value_counts().index:
#these have already been assigned defaults
if i == -1:
continue
names = ifind.index[ifind == i]
c.df.loc[names,'location'] = names.str[:i]
def plot_map(components=["links","stores","storage_units","generators"],bus_size_factor=1.7e10,suffix="all"):
n = pypsa.Network(snakemake.input.network,
override_component_attrs=override_component_attrs)
assign_location(n)
#Drop non-electric buses so they don't clutter the plot
n.buses.drop(n.buses.index[n.buses.carrier != "AC"],inplace=True)
costs = pd.DataFrame(index=n.buses.index)
for comp in components:
getattr(n,comp)["nice_group"] = getattr(n,comp).carrier.map(rename_techs_tyndp)
attr = "e_nom_opt" if comp == "stores" else "p_nom_opt"
costs = pd.concat((costs,(getattr(n,comp).capital_cost*getattr(n,comp)[attr]).groupby((getattr(n,comp).location,getattr(n,comp).nice_group)).sum().unstack().fillna(0.)),axis=1)
print(comp,costs)
costs = costs.groupby(costs.columns,axis=1).sum()
costs.drop(list(costs.columns[(costs == 0.).all()]),axis=1,inplace=True)
new_columns = (preferred_order&costs.columns).append(costs.columns.difference(preferred_order))
costs = costs[new_columns]
#print(costs)
#print(costs.sum())
costs = costs.stack()#.sort_index()
#print(costs)
fig, ax = plt.subplots(subplot_kw={"projection":ccrs.PlateCarree()})
fig.set_size_inches(7,6)
linewidth_factor=2e3
ac_color="gray"
dc_color="m"
#hack because impossible to drop buses...
n.buses.loc["EU gas",["x","y"]] = n.buses.loc["DE0 0",["x","y"]]
n.links.drop(n.links.index[(n.links.carrier != "DC") & (n.links.carrier != "B2B")],inplace=True)
line_widths_exp = pd.concat(dict(Line=(n.lines.s_nom_opt-n.lines.s_nom), Link=(n.links.p_nom_opt-n.links.p_nom)))
#PDF has minimum width, so set these to zero
line_threshold = 500.
line_widths_exp[line_widths_exp < line_threshold] = 0.
#drop non-bus
to_drop = costs.index.levels[0]^n.buses.index
print("dropping non-buses",to_drop)
costs.drop(to_drop,level=0,inplace=True,axis=0)
#make sure they are removed from index
costs.index = pd.MultiIndex.from_tuples(costs.index.values)
n.plot(bus_sizes=costs/bus_size_factor,
bus_colors=snakemake.config['plotting']['tech_colors'],
line_colors=dict(Line=ac_color, Link=dc_color),
line_widths=line_widths_exp/linewidth_factor,
ax=ax)
handles = make_legend_circles_for([5e9, 1e9], scale=bus_size_factor, facecolor="gray")
labels = ["{} bEUR/a".format(s) for s in (5, 1)]
l2 = ax.legend(handles, labels,
loc="upper left", bbox_to_anchor=(0.01, 1.01),
labelspacing=1.0,
framealpha=1.,
title='System cost',
handler_map=make_handler_map_to_scale_circles_as_in(ax))
ax.add_artist(l2)
handles = []
labels = []
for s in (10, 5):
handles.append(plt.Line2D([0],[0],color=ac_color,
linewidth=s*1e3/linewidth_factor))
labels.append("{} GW".format(s))
l1 = l1_1 = ax.legend(handles, labels,
loc="upper left", bbox_to_anchor=(0.30, 1.01),
framealpha=1,
labelspacing=0.8, handletextpad=1.5,
title='Transmission reinforcement')
ax.add_artist(l1_1)
#ax.set_title("Scenario {} with {} transmission".format(snakemake.config['plotting']['scenario_names'][flex],"optimal" if line_limit == "opt" else "no"))
fig.tight_layout()
fig.savefig("{}/{}/graphs/{}-CO0-spatial-costs-{}.pdf".format(snakemake.config['summary_dir'],
snakemake.config['run'],
snakemake.input.scenario.replace(".","p"),
suffix),transparent=True)
def plot_map_without():
n = pypsa.Network(snakemake.input.network,
override_component_attrs=override_component_attrs)
assign_location(n)
#Drop non-electric buses so they don't clutter the plot
n.buses.drop(n.buses.index[n.buses.carrier != "AC"],inplace=True)
fig, ax = plt.subplots(subplot_kw={"projection":ccrs.PlateCarree()})
fig.set_size_inches(7,6)
linewidth_factor=2e3
ac_color="gray"
dc_color="m"
#hack because impossible to drop buses...
n.buses.loc["EU gas",["x","y"]] = n.buses.loc["DE0 0",["x","y"]]
n.links.drop(n.links.index[(n.links.carrier != "DC") & (n.links.carrier != "B2B")],inplace=True)
line_widths_exp = pd.concat(dict(Line=(n.lines.s_nom), Link=(n.links.p_nom)))
#PDF has minimum width, so set these to zero
line_threshold = 0.
line_widths_exp[line_widths_exp < line_threshold] = 0.
line_widths_exp[line_widths_exp > 1e4] = 1e4
n.plot(bus_sizes=10,
bus_colors="k",
line_colors=dict(Line=ac_color, Link=dc_color),
line_widths=line_widths_exp/linewidth_factor,
ax=ax)
handles = []
labels = []
for s in (10, 5):
handles.append(plt.Line2D([0],[0],color=ac_color,
linewidth=s*1e3/linewidth_factor))
labels.append("{} GW".format(s))
l1 = l1_1 = ax.legend(handles, labels,
loc="upper left", bbox_to_anchor=(0.05, 1.01),
framealpha=1,
labelspacing=0.8, handletextpad=1.5,
title='Today\'s transmission')
ax.add_artist(l1_1)
#ax.set_title("Scenario {} with {} transmission".format(snakemake.config['plotting']['scenario_names'][flex],"optimal" if line_limit == "opt" else "no"))
fig.tight_layout()
fig.savefig("{}/{}/graphs/today.pdf".format(snakemake.config['summary_dir'],
snakemake.config['run']),
transparent=True)
def plot_series(carrier="AC"):
n = pypsa.Network(snakemake.input.network,
override_component_attrs=override_component_attrs)
assign_location(n)
assign_carriers(n)
buses = n.buses.index[n.buses.carrier == carrier]
supply = pd.DataFrame(index=n.snapshots)
for c in n.iterate_components(n.branch_components):
for i in range(2):
supply = pd.concat((supply,(-1)*c.pnl["p"+str(i)].loc[:,c.df.index[c.df["bus" + str(i)].isin(buses)]].groupby(c.df.carrier,axis=1).sum()),axis=1)
for c in n.iterate_components(n.one_port_components):
comps = c.df.index[c.df.bus.isin(buses)]
supply = pd.concat((supply,((c.pnl["p"].loc[:,comps]).multiply(c.df.loc[comps,"sign"])).groupby(c.df.carrier,axis=1).sum()),axis=1)
supply = supply.groupby(rename_techs_tyndp,axis=1).sum()
both = supply.columns[(supply < 0.).any() & (supply > 0.).any()]
positive_supply = supply[both]
negative_supply = supply[both]
positive_supply[positive_supply < 0.] = 0.
negative_supply[negative_supply > 0.] = 0.
supply[both] = positive_supply
suffix = " charging"
negative_supply.columns = negative_supply.columns + suffix
supply = pd.concat((supply,negative_supply),axis=1)
fig, ax = plt.subplots()
fig.set_size_inches((8,5))
#14-21.2 for flaute
#19-26.1 for flaute
start = "2013-01-19"
stop = "2013-01-26"
threshold = 10e3
to_drop = supply.columns[(abs(supply) < threshold).all()]
print("dropping",to_drop)
supply.drop(columns=to_drop,inplace=True)
supply.index.name=None
supply = supply/1e3
supply.rename(columns={"electricity" : "electric demand",
"heat" : "heat demand"},
inplace=True)
preferred_order = pd.Index(["electric demand","transmission lines","hydroelectricity","hydro reservoir","run of river","pumped hydro storage","CHP","onshore wind","offshore wind","solar PV","solar thermal","building retrofitting","ground heat pump","air heat pump","resistive heater","OCGT","gas boiler","gas","natural gas","methanation","hydrogen storage","battery storage","hot water storage"])
new_columns = (preferred_order&supply.columns).append(supply.columns.difference(preferred_order))
supply.loc[start:stop,new_columns].plot(ax=ax,kind="area",stacked=True,linewidth=0.,color=[snakemake.config['plotting']['tech_colors'][i.replace(suffix,"")] for i in new_columns])
handles,labels = ax.get_legend_handles_labels()
handles.reverse()
labels.reverse()
new_handles = []
new_labels = []
for i,item in enumerate(labels):
if "charging" not in item:
new_handles.append(handles[i])
new_labels.append(labels[i])
ax.legend(new_handles,new_labels,ncol=3,loc="upper left")
ax.set_xlim([start,stop])
ax.set_ylim([-1300,1900])
ax.grid(True)
ax.set_ylabel("Power [GW]")
fig.tight_layout()
fig.savefig("{}/series-{}-{}-{}.pdf".format(snakemake.config['summary_dir'],carrier,start,stop),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.input.scenario = "lv1.0" #lv1.0, lv1.25, lvopt
snakemake.config["run"] = "190503-es2050-lv"
snakemake.input.network = "{}{}/postnetworks/elec_s_181_{}__Co2L0-3H-T-H-B-I-solar3.nc".format(snakemake.config['results_dir'],
snakemake.config['run'],
snakemake.input.scenario)
plot_map(components=["generators","links","stores","storage_units"],bus_size_factor=1.5e10,suffix="all")
plot_map_without()
#plot_map(components=["generators"],bus_size_factor=1.7e10,suffix="generators")
#plot_map(components=["links","stores","storage_units"],bus_size_factor=4e9,suffix="rest")
#plot_series(carrier="AC")
#plot_series(carrier="heat")