pypsa-eur/scripts/plot_network.py

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"]
override_component_attrs["StorageUnit"].loc["p_dispatch"] = ["series","MW",0.,"Storage discharging.","Output"]
override_component_attrs["StorageUnit"].loc["p_store"] = ["series","MW",0.,"Storage charging.","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):

    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)


    if snakemake.wildcards.lv == "1.0":
        #should be zero
        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)))
    else:
        line_widths_exp = pd.concat(dict(Line=(n.lines.s_nom_opt-n.lines.s_nom_min), Link=(n.links.p_nom_opt-n.links.p_nom_min)))

    #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(snakemake.output.map,transparent=True)



def plot_h2_map():

    n = pypsa.Network(snakemake.input.network,
                      override_component_attrs=override_component_attrs)

    if "H2 pipeline" not in n.links.carrier.unique():
        return

    assign_location(n)

    bus_size_factor=1e5
    linewidth_factor=1e4
    #MW below which not drawn
    line_threshold=1e3
    bus_color="m"
    link_color="c"

    #Drop non-electric buses so they don't clutter the plot
    n.buses.drop(n.buses.index[n.buses.carrier != "AC"],inplace=True)

    elec = n.links.index[n.links.carrier == "H2 Electrolysis"]

    bus_sizes = pd.Series(0., index=n.buses.index)
    bus_sizes.loc[elec.str.replace(" H2 Electrolysis","")] = n.links.loc[elec,"p_nom_opt"].values/bus_size_factor

    #make a fake MultiIndex so that area is correct for legend
    bus_sizes.index = pd.MultiIndex.from_product([bus_sizes.index,["electrolysis"]])

    n.links.drop(n.links.index[n.links.carrier != "H2 pipeline"],inplace=True)

    link_widths = n.links.p_nom_opt/linewidth_factor
    link_widths[n.links.p_nom_opt < line_threshold] = 0.

    n.links.bus0 = n.links.bus0.str.replace(" H2","")
    n.links.bus1 = n.links.bus1.str.replace(" H2","")

    print(link_widths.sort_values())

    print(n.links[["bus0","bus1"]])

    fig, ax = plt.subplots(subplot_kw={"projection":ccrs.PlateCarree()})

    fig.set_size_inches(7,6)

    n.plot(bus_sizes=bus_sizes,
           bus_colors={"electrolysis" : bus_color},
           line_colors=dict(Link=link_color),
           line_widths={"Link" : link_widths},
           branch_components=["Link"],
           ax=ax)

    handles = make_legend_circles_for([50000, 10000], scale=bus_size_factor, facecolor=bus_color)
    labels = ["{} GW".format(s) for s in (50, 10)]
    l2 = ax.legend(handles, labels,
                       loc="upper left", bbox_to_anchor=(0.01, 1.01),
                       labelspacing=1.0,
                       framealpha=1.,
                       title='Electrolyzer capacity',
                       handler_map=make_handler_map_to_scale_circles_as_in(ax))
    ax.add_artist(l2)

    handles = []
    labels = []

    for s in (50, 10):
        handles.append(plt.Line2D([0],[0],color=link_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='H2 pipeline capacity')
    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(snakemake.output.map.replace("-costs-all","-h2_network"),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)


    if snakemake.wildcards.lv == "1.0":
        line_widths_exp = pd.concat(dict(Line=(n.lines.s_nom), Link=(n.links.p_nom)))
    else:
        line_widths_exp = pd.concat(dict(Line=(n.lines.s_nom_min), Link=(n.links.p_nom_min)))

    #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(snakemake.output.today,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.output = 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)
        snakemake.output.network = "{}{}/maps/elec_s_181_{}__Co2L0-3H-T-H-B-I-solar3.pdf".format(snakemake.config['results_dir'],
                                                                                                 snakemake.config['run'],
                                                                                                 snakemake.input.scenario)


    plot_map(components=["generators","links","stores","storage_units"],bus_size_factor=1.5e10)

    plot_h2_map()

    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")