diff --git a/etrago/tools/calc_results.py b/etrago/tools/calc_results.py
index fca593a9..f3e3338d 100755
--- a/etrago/tools/calc_results.py
+++ b/etrago/tools/calc_results.py
@@ -22,7 +22,8 @@
 calc_results.py defines methods to calculate results of eTraGo
 """
 import os
-if 'READTHEDOCS' not in os.environ:
+
+if "READTHEDOCS" not in os.environ:
     import time
     import logging
 
@@ -31,16 +32,18 @@
 
     logger = logging.getLogger(__name__)
 
-__copyright__ = ("Flensburg University of Applied Sciences, "
-                 "Europa-Universität Flensburg, "
-                 "Centre for Sustainable Energy Systems, "
-                 "DLR-Institute for Networked Energy Systems")
+__copyright__ = (
+    "Flensburg University of Applied Sciences, "
+    "Europa-Universität Flensburg, "
+    "Centre for Sustainable Energy Systems, "
+    "DLR-Institute for Networked Energy Systems"
+)
 __license__ = "GNU Affero General Public License Version 3 (AGPL-3.0)"
 __author__ = "ulfmueller, s3pp, wolfbunke, mariusves, lukasol"
 
 
 def _calc_storage_expansion(self):
-    """ Function that calulates storage expansion in MW
+    """Function that calulates storage expansion in MW
 
     Returns
     -------
@@ -48,14 +51,18 @@ def _calc_storage_expansion(self):
         storage expansion in MW
 
     """
-    return (self.network.storage_units.p_nom_opt -
-            self.network.storage_units.p_nom_min
-            )[self.network.storage_units.p_nom_extendable]\
-                .groupby(self.network.storage_units.carrier).sum()
+    return (
+        (
+            self.network.storage_units.p_nom_opt
+            - self.network.storage_units.p_nom_min
+        )[self.network.storage_units.p_nom_extendable]
+        .groupby(self.network.storage_units.carrier)
+        .sum()
+    )
 
 
 def _calc_store_expansion(self):
-    """ Function that calulates store expansion in MW
+    """Function that calulates store expansion in MW
 
     Returns
     -------
@@ -63,13 +70,13 @@ def _calc_store_expansion(self):
         store expansion in MW
 
     """
-    return (self.network.stores.e_nom_opt -
-            self.network.stores.e_nom_min
-            )[self.network.stores.e_nom_extendable]
+    return (self.network.stores.e_nom_opt - self.network.stores.e_nom_min)[
+        self.network.stores.e_nom_extendable
+    ]
 
 
 def _calc_sectorcoupling_link_expansion(self):
-    """ Function that calulates expansion of sectorcoupling links in MW
+    """Function that calulates expansion of sectorcoupling links in MW
 
     Returns
     -------
@@ -81,21 +88,21 @@ def _calc_sectorcoupling_link_expansion(self):
 
     links = [0, 0, 0, 0]
 
-    l1 = ext_links[ext_links.carrier=='H2_to_power']
-    l2 = ext_links[ext_links.carrier=='power_to_H2']
-    l3 = ext_links[ext_links.carrier=='H2_to_CH4']
-    l4 = ext_links[ext_links.carrier=='CH4_to_H2']
+    l1 = ext_links[ext_links.carrier == "H2_to_power"]
+    l2 = ext_links[ext_links.carrier == "power_to_H2"]
+    l3 = ext_links[ext_links.carrier == "H2_to_CH4"]
+    l4 = ext_links[ext_links.carrier == "CH4_to_H2"]
 
-    links[0] = (l1.p_nom_opt-l1.p_nom_min).sum()
-    links[1] = (l2.p_nom_opt-l2.p_nom_min).sum()
-    links[2] = (l3.p_nom_opt-l3.p_nom_min).sum()
-    links[3] = (l4.p_nom_opt-l4.p_nom_min).sum()
+    links[0] = (l1.p_nom_opt - l1.p_nom_min).sum()
+    links[1] = (l2.p_nom_opt - l2.p_nom_min).sum()
+    links[2] = (l3.p_nom_opt - l3.p_nom_min).sum()
+    links[3] = (l4.p_nom_opt - l4.p_nom_min).sum()
 
     return links
 
 
 def _calc_network_expansion(self):
-    """ Function that calulates electrical network expansion in MW
+    """Function that calulates electrical network expansion in MW
 
     Returns
     -------
@@ -106,21 +113,20 @@ def _calc_network_expansion(self):
 
     network = self.network
 
-    lines = (network.lines.s_nom_opt -
-            network.lines.s_nom_min
-            )[network.lines.s_nom_extendable]
+    lines = (network.lines.s_nom_opt - network.lines.s_nom_min)[
+        network.lines.s_nom_extendable
+    ]
 
     ext_links = network.links[network.links.p_nom_extendable]
-    ext_dc_lines = ext_links[ext_links.carrier=='DC']
+    ext_dc_lines = ext_links[ext_links.carrier == "DC"]
 
-    dc_links = (ext_dc_lines.p_nom_opt -
-            ext_dc_lines.p_nom_min)
+    dc_links = ext_dc_lines.p_nom_opt - ext_dc_lines.p_nom_min
 
     return lines, dc_links
 
 
 def calc_investment_cost(self):
-    """ Function that calulates overall annualized investment costs.
+    """Function that calulates overall annualized investment costs.
 
     Returns
     -------
@@ -141,30 +147,40 @@ def calc_investment_cost(self):
     ext_lines = network.lines[network.lines.s_nom_extendable]
     ext_trafos = network.transformers[network.transformers.s_nom_extendable]
     ext_links = network.links[network.links.p_nom_extendable]
-    ext_dc_lines = ext_links[ext_links.carrier=='DC']
+    ext_dc_lines = ext_links[ext_links.carrier == "DC"]
 
     if not ext_lines.empty:
-        network_costs[0] = ((ext_lines.s_nom_opt-ext_lines.s_nom_min
-                             )*ext_lines.capital_cost).sum()
+        network_costs[0] = (
+            (ext_lines.s_nom_opt - ext_lines.s_nom_min)
+            * ext_lines.capital_cost
+        ).sum()
 
     if not ext_trafos.empty:
-        network_costs[0] = network_costs[0]+((
-            ext_trafos.s_nom_opt-ext_trafos.s_nom
-            )*ext_trafos.capital_cost).sum()
+        network_costs[0] = (
+            network_costs[0]
+            + (
+                (ext_trafos.s_nom_opt - ext_trafos.s_nom)
+                * ext_trafos.capital_cost
+            ).sum()
+        )
 
     if not ext_dc_lines.empty:
-        network_costs[1] = ((ext_dc_lines.p_nom_opt-ext_dc_lines.p_nom_min
-                             )*ext_dc_lines.capital_cost).sum()
+        network_costs[1] = (
+            (ext_dc_lines.p_nom_opt - ext_dc_lines.p_nom_min)
+            * ext_dc_lines.capital_cost
+        ).sum()
 
     # links in other sectors / coupling different sectors
 
     link_costs = 0
 
-    ext_links = ext_links[ext_links.carrier!='DC']
+    ext_links = ext_links[ext_links.carrier != "DC"]
 
     if not ext_links.empty:
-        link_costs = ((ext_links.p_nom_opt-ext_links.p_nom_min
-                             )*ext_links.capital_cost).sum()
+        link_costs = (
+            (ext_links.p_nom_opt - ext_links.p_nom_min)
+            * ext_links.capital_cost
+        ).sum()
 
     # storage and store costs
 
@@ -174,12 +190,10 @@ def calc_investment_cost(self):
     ext_store = network.stores[network.stores.e_nom_extendable]
 
     if not ext_storage.empty:
-        sto_costs[0] = (ext_storage.p_nom_opt*
-                         ext_storage.capital_cost).sum()
+        sto_costs[0] = (ext_storage.p_nom_opt * ext_storage.capital_cost).sum()
 
     if not ext_store.empty:
-        sto_costs[1] = (ext_store.e_nom_opt*
-                         ext_store.capital_cost).sum()
+        sto_costs[1] = (ext_store.e_nom_opt * ext_store.capital_cost).sum()
 
     return network_costs, link_costs, sto_costs
 
@@ -196,21 +210,35 @@ def calc_marginal_cost(self):
 
     """
     network = self.network
-    gen = network.generators_t.p.mul(
-        network.snapshot_weightings.objective, axis=0).sum(axis=0).mul(
-            network.generators.marginal_cost).sum()
-    link = abs(network.links_t.p0).mul(
-        network.snapshot_weightings.objective, axis=0).sum(axis=0).mul(
-            network.links.marginal_cost).sum()
-    stor = network.storage_units_t.p.mul(
-        network.snapshot_weightings.objective, axis=0).sum(axis=0).mul(
-            network.storage_units.marginal_cost).sum()
+    gen = (
+        network.generators_t.p.mul(
+            network.snapshot_weightings.objective, axis=0
+        )
+        .sum(axis=0)
+        .mul(network.generators.marginal_cost)
+        .sum()
+    )
+    link = (
+        abs(network.links_t.p0)
+        .mul(network.snapshot_weightings.objective, axis=0)
+        .sum(axis=0)
+        .mul(network.links.marginal_cost)
+        .sum()
+    )
+    stor = (
+        network.storage_units_t.p.mul(
+            network.snapshot_weightings.objective, axis=0
+        )
+        .sum(axis=0)
+        .mul(network.storage_units.marginal_cost)
+        .sum()
+    )
     marginal_cost = gen + link + stor
     return marginal_cost
 
 
 def calc_etrago_results(self):
-    """ Function that calculates main results of grid optimization
+    """Function that calculates main results of grid optimization
     and adds them to Etrago object.
 
     Returns
@@ -218,99 +246,137 @@ def calc_etrago_results(self):
     None.
 
     """
-    self.results = pd.DataFrame(columns=['unit', 'value'],
-                                index=['annual system costs',
-                                       'annual investment costs',
-                                       'annual marginal costs',
-                                       'annual electrical grid investment costs',
-                                       'annual ac grid investment costs',
-                                       'annual dc grid investment costs',
-                                       'annual links investment costs',
-                                       'annual storage+store investment costs',
-                                       'annual electrical storage investment costs',
-                                       'annual store investment costs',
-                                       'battery storage expansion',
-                                       'store expansion',
-                                       'H2 store expansion',
-                                       'CH4 store expansion',
-                                       'heat store expansion',
-                                       'storage+store expansion',
-                                       'fuel cell links expansion',
-                                       'electrolyzer links expansion',
-                                       'methanisation links expansion',
-                                       'Steam Methane Reformation links expansion',
-                                       'abs. electrical grid expansion',
-                                       'abs. electrical ac grid expansion',
-                                       'abs. electrical dc grid expansion',
-                                       'rel. electrical ac grid expansion',
-                                       'rel. electrical dc grid expansion'])
-
-    self.results.unit[self.results.index.str.contains('cost')] = 'EUR/a'
-    self.results.unit[self.results.index.str.contains('expansion')] = 'MW'
-    self.results.unit[self.results.index.str.contains('rel.')] = 'p.u.'
+    self.results = pd.DataFrame(
+        columns=["unit", "value"],
+        index=[
+            "annual system costs",
+            "annual investment costs",
+            "annual marginal costs",
+            "annual electrical grid investment costs",
+            "annual ac grid investment costs",
+            "annual dc grid investment costs",
+            "annual links investment costs",
+            "annual storage+store investment costs",
+            "annual electrical storage investment costs",
+            "annual store investment costs",
+            "battery storage expansion",
+            "store expansion",
+            "H2 store expansion",
+            "CH4 store expansion",
+            "heat store expansion",
+            "storage+store expansion",
+            "fuel cell links expansion",
+            "electrolyzer links expansion",
+            "methanisation links expansion",
+            "Steam Methane Reformation links expansion",
+            "abs. electrical grid expansion",
+            "abs. electrical ac grid expansion",
+            "abs. electrical dc grid expansion",
+            "rel. electrical ac grid expansion",
+            "rel. electrical dc grid expansion",
+        ],
+    )
+
+    self.results.unit[self.results.index.str.contains("cost")] = "EUR/a"
+    self.results.unit[self.results.index.str.contains("expansion")] = "MW"
+    self.results.unit[self.results.index.str.contains("rel.")] = "p.u."
 
     # system costs
 
-    self.results.value['annual ac grid investment costs'] = calc_investment_cost(self)[0][0]
-    self.results.value['annual dc grid investment costs'] = calc_investment_cost(self)[0][1]
-    self.results.value['annual electrical grid investment costs'] = sum(calc_investment_cost(self)[0])
-
-    self.results.value['annual links investment costs'] = calc_investment_cost(self)[1]
-
-    self.results.value['annual electrical storage investment costs'] = calc_investment_cost(self)[2][0]
-    self.results.value['annual store investment costs'] = calc_investment_cost(self)[2][1]
-    self.results.value['annual storage+store investment costs'] = sum(calc_investment_cost(self)[2])
-
-
-    self.results.value['annual investment costs'] = \
-        sum(calc_investment_cost(self)[0]) + calc_investment_cost(self)[1] + sum(calc_investment_cost(self)[2])
-    self.results.value['annual marginal costs'] = calc_marginal_cost(self)
-
-    self.results.value['annual system costs'] = \
-        self.results.value['annual investment costs'] + self.results.value['annual marginal costs']
+    self.results.value[
+        "annual ac grid investment costs"
+    ] = calc_investment_cost(self)[0][0]
+    self.results.value[
+        "annual dc grid investment costs"
+    ] = calc_investment_cost(self)[0][1]
+    self.results.value["annual electrical grid investment costs"] = sum(
+        calc_investment_cost(self)[0]
+    )
+
+    self.results.value["annual links investment costs"] = calc_investment_cost(
+        self
+    )[1]
+
+    self.results.value[
+        "annual electrical storage investment costs"
+    ] = calc_investment_cost(self)[2][0]
+    self.results.value["annual store investment costs"] = calc_investment_cost(
+        self
+    )[2][1]
+    self.results.value["annual storage+store investment costs"] = sum(
+        calc_investment_cost(self)[2]
+    )
+
+    self.results.value["annual investment costs"] = (
+        sum(calc_investment_cost(self)[0])
+        + calc_investment_cost(self)[1]
+        + sum(calc_investment_cost(self)[2])
+    )
+    self.results.value["annual marginal costs"] = calc_marginal_cost(self)
+
+    self.results.value["annual system costs"] = (
+        self.results.value["annual investment costs"]
+        + self.results.value["annual marginal costs"]
+    )
 
     # storage and store expansion
 
     network = self.network
 
     if not network.storage_units[network.storage_units.p_nom_extendable].empty:
-
-        self.results.value['battery storage expansion'] = \
-            _calc_storage_expansion(self).sum()
+        self.results.value[
+            "battery storage expansion"
+        ] = _calc_storage_expansion(self).sum()
 
         store = _calc_store_expansion(self)
-        self.results.value['store expansion'] = store.sum()
-        self.results.value['H2 store expansion'] = \
-            store[store.index.str.contains('H2')].sum()
-        self.results.value['CH4 store expansion'] = \
-            store[store.index.str.contains('CH4')].sum()
-        self.results.value['heat store expansion'] = \
-            store[store.index.str.contains('heat')].sum()
-
-        self.results.value['storage+store expansion'] = \
-            self.results.value['battery storage expansion'] + self.results.value['store expansion']
+        self.results.value["store expansion"] = store.sum()
+        self.results.value["H2 store expansion"] = store[
+            store.index.str.contains("H2")
+        ].sum()
+        self.results.value["CH4 store expansion"] = store[
+            store.index.str.contains("CH4")
+        ].sum()
+        self.results.value["heat store expansion"] = store[
+            store.index.str.contains("heat")
+        ].sum()
+
+        self.results.value["storage+store expansion"] = (
+            self.results.value["battery storage expansion"]
+            + self.results.value["store expansion"]
+        )
 
     # links expansion
 
     if not network.links[network.links.p_nom_extendable].empty:
-
         links = _calc_sectorcoupling_link_expansion(self)
-        self.results.value['fuel cell links expansion'] = links[0]
-        self.results.value['electrolyzer links expansion'] = links[1]
-        self.results.value['methanisation links expansion'] = links[2]
-        self.results.value['Steam Methane Reformation links expansion'] = links[3]
+        self.results.value["fuel cell links expansion"] = links[0]
+        self.results.value["electrolyzer links expansion"] = links[1]
+        self.results.value["methanisation links expansion"] = links[2]
+        self.results.value[
+            "Steam Methane Reformation links expansion"
+        ] = links[3]
 
     # grid expansion
 
     if not network.lines[network.lines.s_nom_extendable].empty:
-
-        self.results.value['abs. electrical ac grid expansion'] = _calc_network_expansion(self)[0].sum()
-        self.results.value['abs. electrical dc grid expansion'] = _calc_network_expansion(self)[1].sum()
-        self.results.value['abs. electrical grid expansion'] = self.results.value['abs. electrical ac grid expansion'] + self.results.value['abs. electrical dc grid expansion']
+        self.results.value[
+            "abs. electrical ac grid expansion"
+        ] = _calc_network_expansion(self)[0].sum()
+        self.results.value[
+            "abs. electrical dc grid expansion"
+        ] = _calc_network_expansion(self)[1].sum()
+        self.results.value["abs. electrical grid expansion"] = (
+            self.results.value["abs. electrical ac grid expansion"]
+            + self.results.value["abs. electrical dc grid expansion"]
+        )
 
         ext_lines = network.lines[network.lines.s_nom_extendable]
         ext_links = network.links[network.links.p_nom_extendable]
-        ext_dc_lines = ext_links[ext_links.carrier=='DC']
-
-        self.results.value['rel. electrical ac grid expansion'] = (_calc_network_expansion(self)[0].sum() / ext_lines.s_nom.sum())
-        self.results.value['rel. electrical dc grid expansion'] = (_calc_network_expansion(self)[1].sum() / ext_dc_lines.p_nom.sum())
\ No newline at end of file
+        ext_dc_lines = ext_links[ext_links.carrier == "DC"]
+
+        self.results.value["rel. electrical ac grid expansion"] = (
+            _calc_network_expansion(self)[0].sum() / ext_lines.s_nom.sum()
+        )
+        self.results.value["rel. electrical dc grid expansion"] = (
+            _calc_network_expansion(self)[1].sum() / ext_dc_lines.p_nom.sum()
+        )
diff --git a/etrago/tools/execute.py b/etrago/tools/execute.py
index f7a13412..8c2441bc 100755
--- a/etrago/tools/execute.py
+++ b/etrago/tools/execute.py
@@ -511,6 +511,107 @@ def dispatch_disaggregation(self):
         logger.info("Time for LOPF [min]: {}".format(round(z, 2)))
 
 
+def import_gen_from_links(network):
+    """
+    create gas generators from links in order to not lose them when
+    dropping non-electric carriers
+    """
+    # Discard all generators < 1kW
+    discard_gen = network.links[network.links["p_nom"] <= 0.001].index
+    network.links.drop(discard_gen, inplace=True)
+    for df in network.links_t:
+        if not network.links_t[df].empty:
+            network.links_t[df].drop(
+                columns=discard_gen.values, inplace=True, errors="ignore"
+            )
+
+    gas_to_add = network.links[
+        network.links.carrier.isin(
+            [
+                "central_gas_CHP",
+                "OCGT",
+                "H2_to_power",
+                "industrial_gas_CHP",
+            ]
+        )
+    ].copy()
+
+    # Drop generators from the links table
+    network.links.drop(gas_to_add.index, inplace=True)
+
+    gas_to_add.rename(columns={"bus1": "bus"}, inplace=True)
+
+    # Create generators' names like in network.generators
+    gas_to_add["Generator"] = (
+        gas_to_add["bus"] + " " + gas_to_add.index + gas_to_add["carrier"]
+    )
+    gas_to_add_orig = gas_to_add.copy()
+    gas_to_add.set_index("Generator", drop=True, inplace=True)
+    gas_to_add = gas_to_add[
+        gas_to_add.columns[
+            gas_to_add.columns.isin(network.generators.columns)
+        ]
+    ]
+
+    network.import_components_from_dataframe(gas_to_add, "Generator")
+
+    # Dealing with generators_t
+    columns_new = network.links_t.p1.columns[
+        network.links_t.p1.columns.isin(gas_to_add_orig.index)
+    ]
+
+    new_gen_t = network.links_t.p1[columns_new] * -1
+    new_gen_t.rename(columns=gas_to_add_orig["Generator"], inplace=True)
+    network.generators_t.p = network.generators_t.p.join(new_gen_t)
+
+    # Drop generators from the links_t table
+    for df in network.links_t:
+        if not network.links_t[df].empty:
+            network.links_t[df].drop(
+                columns=gas_to_add_orig.index,
+                inplace=True,
+                errors="ignore",
+            )
+
+    # Group generators per bus if needed
+    if not (
+        network.generators.groupby(["bus", "carrier"]).p_nom.count() == 1
+    ).all():
+        network.generators["weight"] = network.generators.p_nom
+        df, df_t = aggregategenerators(
+            network,
+            busmap=pd.Series(
+                index=network.buses.index, data=network.buses.index
+            ),
+            custom_strategies=strategies_generators(),
+        )
+
+        # Keep control arguments from generators
+        control = network.generators.groupby(
+            ["bus", "carrier"]
+        ).control.first()
+        control.index = (
+            control.index.get_level_values(0)
+            + " "
+            + control.index.get_level_values(1)
+        )
+        df.control = control
+
+        # Drop non-aggregated generators
+        network.mremove("Generator", network.generators.index)
+
+        # Insert aggregated generators and time series
+        network.import_components_from_dataframe(df, "Generator")
+
+        for attr, data in df_t.items():
+            if not data.empty:
+                network.import_series_from_dataframe(
+                    data, "Generator", attr
+                )
+
+    return
+
+
 def run_pf_post_lopf(self):
     """
     Function that runs pf_post_lopf according to arguments.
@@ -659,117 +760,6 @@ def drop_foreign_components(network):
 
         return foreign_bus, foreign_comp, foreign_series
 
-    def import_gen_from_links(network):
-        """
-        Creates gas generators from links in order to not lose them when
-        dropping non-electric carriers.
-
-        Parameters
-        ----------
-        network : pypsa.Network object
-            Container for all network components.
-
-        Returns
-        -------
-        None.
-
-        """
-
-        # Discard all generators < 1kW
-        discard_gen = network.links[network.links["p_nom"] <= 0.001].index
-        network.links.drop(discard_gen, inplace=True)
-        for df in network.links_t:
-            if not network.links_t[df].empty:
-                network.links_t[df].drop(
-                    columns=discard_gen.values, inplace=True, errors="ignore"
-                )
-
-        gas_to_add = network.links[
-            network.links.carrier.isin(
-                [
-                    "central_gas_CHP",
-                    "OCGT",
-                    "H2_to_power",
-                    "industrial_gas_CHP",
-                ]
-            )
-        ].copy()
-
-        # Drop generators from the links table
-        network.links.drop(gas_to_add.index, inplace=True)
-
-        gas_to_add.rename(columns={"bus1": "bus"}, inplace=True)
-
-        # Create generators' names like in network.generators
-        gas_to_add["Generator"] = (
-            gas_to_add["bus"] + " " + gas_to_add.index + gas_to_add["carrier"]
-        )
-        gas_to_add_orig = gas_to_add.copy()
-        gas_to_add.set_index("Generator", drop=True, inplace=True)
-        gas_to_add = gas_to_add[
-            gas_to_add.columns[
-                gas_to_add.columns.isin(network.generators.columns)
-            ]
-        ]
-
-        network.import_components_from_dataframe(gas_to_add, "Generator")
-
-        # Dealing with generators_t
-        columns_new = network.links_t.p1.columns[
-            network.links_t.p1.columns.isin(gas_to_add_orig.index)
-        ]
-
-        new_gen_t = network.links_t.p1[columns_new] * -1
-        new_gen_t.rename(columns=gas_to_add_orig["Generator"], inplace=True)
-        network.generators_t.p = network.generators_t.p.join(new_gen_t)
-
-        # Drop generators from the links_t table
-        for df in network.links_t:
-            if not network.links_t[df].empty:
-                network.links_t[df].drop(
-                    columns=gas_to_add_orig.index,
-                    inplace=True,
-                    errors="ignore",
-                )
-
-        # Group generators per bus if needed
-        if not (
-            network.generators.groupby(["bus", "carrier"]).p_nom.count() == 1
-        ).all():
-            network.generators["weight"] = network.generators.p_nom
-            df, df_t = aggregategenerators(
-                network,
-                busmap=pd.Series(
-                    index=network.buses.index, data=network.buses.index
-                ),
-                custom_strategies=strategies_generators(),
-            )
-
-            # Keep control arguments from generators
-            control = network.generators.groupby(
-                ["bus", "carrier"]
-            ).control.first()
-            control.index = (
-                control.index.get_level_values(0)
-                + " "
-                + control.index.get_level_values(1)
-            )
-            df.control = control
-
-            # Drop non-aggregated generators
-            network.mremove("Generator", network.generators.index)
-
-            # Insert aggregated generators and time series
-            network.import_components_from_dataframe(df, "Generator")
-
-            for attr, data in df_t.items():
-                if not data.empty:
-                    network.import_series_from_dataframe(
-                        data, "Generator", attr
-                    )
-
-        return
-
     x = time.time()
     network = etrago.network
     args = etrago.args
diff --git a/etrago/tools/plot.py b/etrago/tools/plot.py
index d3ae2068..0ab5b827 100644
--- a/etrago/tools/plot.py
+++ b/etrago/tools/plot.py
@@ -25,7 +25,11 @@
 import logging
 import os
 
+from etrago.tools.execute import import_gen_from_links
 from matplotlib import pyplot as plt
+from matplotlib.legend_handler import HandlerPatch
+from matplotlib.patches import Circle, Ellipse
+from pyproj import Proj, transform
 import matplotlib
 import matplotlib.patches as mpatches
 import numpy as np
@@ -33,10 +37,7 @@
 
 cartopy_present = True
 try:
-    import cartopy
     import cartopy.crs as ccrs
-    import cartopy.mpl.geoaxes
-    import requests
 except ImportError:
     cartopy_present = False
 from pypsa.plot import draw_map_cartopy
@@ -135,6 +136,7 @@ def coloring():
         "power_to_H2": "cyan",
         "H2_overground": "cyan",
         "H2_underground": "cyan",
+        "H2": "cyan",
         "dsm-cts": "dodgerblue",
         "dsm-ind-osm": "dodgerblue",
         "dsm-ind-sites": "dodgerblue",
@@ -143,6 +145,8 @@ def coloring():
         "central_resistive_heater": "blueviolet",
         "rural_heat_pump": "violet",
         "CH4": "yellow",
+        "CH4_biogas": "yellow",
+        "CH4_NG": "yellow",
         "CH4_to_H2": "yellowgreen",
         "industrial_gas_CHP": "olive",
         "rural_gas_boiler": "sandybrown",
@@ -164,6 +168,7 @@ def coloring():
         "H2_grid": "green",
         "H2_saltcavern": "darkgreen",
         "central_heat_store": "firebrick",
+        "heat": "firebrick",
         "rural_heat_store": "salmon",
         "AC": "blue",
         "nuclear": "palegreen",
@@ -191,6 +196,7 @@ def coloring():
         "pumped_hydro": "indigo",
         "BEV charger": "indigo",
         "BEV_charger": "indigo",
+        "others": "dimgrey",
     }
 
     return colors
@@ -226,7 +232,7 @@ def plot_line_loading_diff(networkA, networkB, timestep=0, osm=False):
         * 'zoom' : resolution of osm
 
     """
-    if osm != False:
+    if osm is not False:
         if set_epsg_network.counter == 0:
             set_epsg_network(networkA)
             set_epsg_network(networkB)
@@ -363,7 +369,7 @@ def network_expansion_diff(
         * 'zoom' : resolution of osm
 
     """
-    if osm != False:
+    if osm is not False:
         if set_epsg_network.counter == 0:
             set_epsg_network(networkA)
             set_epsg_network(networkB)
@@ -667,7 +673,7 @@ def plot_voltage(network, boundaries=[], osm=False):
     -------
     Plot
     """
-    if osm != False:
+    if osm is not False:
         if set_epsg_network.counter == 0:
             set_epsg_network(network)
         plot_osm(osm["x"], osm["y"], osm["zoom"])
@@ -774,6 +780,16 @@ def calc_dispatch_per_carrier(network, timesteps):
 
     """
 
+    import_gen_from_links(network)
+
+    ac_buses = network.buses[network.buses.carrier == "AC"].index
+    network.generators = network.generators[
+        network.generators.bus.isin(ac_buses)
+    ]
+    network.generators_t.p = network.generators_t.p.loc[
+        :, network.generators_t.p.columns.isin(network.generators.index)
+    ]
+
     index = [
         (network.generators.bus[idx], network.generators.carrier[idx])
         for idx in network.generators.index
@@ -802,7 +818,16 @@ def calc_dispatch_per_carrier(network, timesteps):
     return dist
 
 
-def calc_storage_expansion_per_bus(network):
+def calc_storage_expansion_per_bus(
+    network,
+    carriers=[
+        "battery",
+        "H2_overground",
+        "H2_underground",
+        "rural_heat_store",
+        "central_heat_store",
+    ],
+):
     """Function that calculates storage expansion per bus and technology
 
     Parameters
@@ -816,60 +841,10 @@ def calc_storage_expansion_per_bus(network):
         storage expansion per bus and technology
 
     """
-
-    batteries = network.storage_units[
-        network.storage_units.carrier == "battery"
-    ]
-    h2_overground = network.stores[network.stores.carrier == "H2_overground"]
-    h2_underground = network.stores[network.stores.carrier == "H2_underground"]
-    rural_heat = network.stores[network.stores.carrier == "rural_heat_store"]
-    central_heat = network.stores[
-        network.stores.carrier == "central_heat_store"
-    ]
-    # hydrogen = network.storage_units[network.storage_units.carrier ==
-    #                                  'extendable_hydrogen_storage']
-    battery_distribution = (
-        network.storage_units.p_nom_opt[batteries.index]
-        .groupby(network.storage_units.bus)
-        .sum()
-        .reindex(network.buses.index, fill_value=0.0)
-    )
-    h2_over_distribution = (
-        network.stores.e_nom_opt[h2_overground.index]
-        .groupby(network.stores.bus)
-        .sum()
-        .reindex(network.buses.index, fill_value=0.0)
-    )
-    h2_under_distribution = (
-        network.stores.e_nom_opt[h2_underground.index]
-        .groupby(network.stores.bus)
-        .sum()
-        .reindex(network.buses.index, fill_value=0.0)
-    )
-    rural_heat_distribution = (
-        network.stores.e_nom_opt[rural_heat.index]
-        .groupby(network.stores.bus)
-        .sum()
-        .reindex(network.buses.index, fill_value=0.0)
-    )
-    central_heat_distribution = (
-        network.stores.e_nom_opt[central_heat.index]
-        .groupby(network.stores.bus)
-        .sum()
-        .reindex(network.buses.index, fill_value=0.0)
-    )
-    # hydrogen_distribution =\
-    #     network.storage_units.p_nom_opt[hydrogen.index].groupby(
-    #         network.storage_units.bus).sum().reindex(
-    #             network.buses.index, fill_value=0.)
     index = [(idx, "battery") for idx in network.buses.index]
-    for c in [
-        "H2_overground",
-        "H2_underground",
-        "rural_heat_store",
-        "central_heat_store",
-    ]:
-        index.extend([(idx, c) for idx in network.buses.index])
+    for c in carriers:
+        if c != "battery":
+            index.extend([(idx, c) for idx in network.buses.index])
     # index.extend([(idx, 'hydrogen_storage') for idx in network.buses.index])
 
     dist = pd.Series(
@@ -877,25 +852,107 @@ def calc_storage_expansion_per_bus(network):
         dtype=float,
     )
 
-    dist.iloc[
-        dist.index.get_level_values("carrier") == "battery"
-    ] = battery_distribution.sort_index().values
-    dist.iloc[
-        dist.index.get_level_values("carrier") == "H2_overground"
-    ] = h2_over_distribution.sort_index().values
-    dist.iloc[
-        dist.index.get_level_values("carrier") == "H2_underground"
-    ] = h2_under_distribution.sort_index().values
-    dist.iloc[
-        dist.index.get_level_values("carrier") == "rural_heat_store"
-    ] = rural_heat_distribution.sort_index().values
-    dist.iloc[
-        dist.index.get_level_values("carrier") == "central_heat_store"
-    ] = central_heat_distribution.sort_index().values
-    # dist.iloc[dist.index.get_level_values('carrier') == 'hydrogen_storage'] = \
-    #         hydrogen_distribution.sort_index().values
-    # network.carriers.color['hydrogen_storage'] = 'orange'
-    # network.carriers.color['battery_storage'] = 'blue'
+    if "battery" in carriers:
+        batteries = network.storage_units[
+            network.storage_units.carrier == "battery"
+        ]
+        battery_distribution = (
+            (
+                network.storage_units.p_nom_opt[batteries.index]
+                - network.storage_units.p_nom_min[batteries.index]
+            )
+            .groupby(network.storage_units.bus)
+            .sum()
+            .reindex(network.buses.index, fill_value=0.0)
+        ).mul(6)
+
+        battery_distribution.index = pd.MultiIndex.from_tuples(
+            [(idx, "battery") for idx in battery_distribution.index]
+        )
+
+        dist.loc[
+            dist.index.get_level_values("carrier") == "battery"
+        ] = battery_distribution
+    if "H2_overground" in carriers:
+        h2_overground = network.stores[
+            network.stores.carrier == "H2_overground"
+        ]
+        h2_over_distribution = (
+            network.stores.e_nom_opt[h2_overground.index]
+            .groupby(network.stores.bus)
+            .sum()
+            .reindex(network.buses.index, fill_value=0.0)
+        )
+
+        h2_over_distribution.index = pd.MultiIndex.from_tuples(
+            [(idx, "H2_overground") for idx in h2_over_distribution.index]
+        )
+
+        dist.loc[
+            dist.index.get_level_values("carrier") == "H2_overground"
+        ] = h2_over_distribution
+
+    if "H2_overground" in carriers:
+        h2_underground = network.stores[
+            network.stores.carrier == "H2_underground"
+        ]
+        h2_under_distribution = (
+            network.stores.e_nom_opt[h2_underground.index]
+            .groupby(network.stores.bus)
+            .sum()
+            .reindex(network.buses.index, fill_value=0.0)
+        )
+
+        h2_under_distribution.index = pd.MultiIndex.from_tuples(
+            [(idx, "H2_underground") for idx in h2_under_distribution.index]
+        )
+
+        dist.loc[
+            dist.index.get_level_values("carrier") == "H2_underground"
+        ] = h2_under_distribution
+
+    if "rural_heat_store" in carriers:
+        rural_heat = network.stores[
+            network.stores.carrier == "rural_heat_store"
+        ]
+        rural_heat_distribution = (
+            network.stores.e_nom_opt[rural_heat.index]
+            .groupby(network.stores.bus)
+            .sum()
+            .reindex(network.buses.index, fill_value=0.0)
+        )
+
+        rural_heat_distribution.index = pd.MultiIndex.from_tuples(
+            [
+                (idx, "rural_heat_store")
+                for idx in rural_heat_distribution.index
+            ]
+        )
+
+        dist.loc[
+            dist.index.get_level_values("carrier") == "rural_heat_store"
+        ] = rural_heat_distribution
+    if "central_heat_store" in carriers:
+        central_heat = network.stores[
+            network.stores.carrier == "central_heat_store"
+        ]
+        central_heat_distribution = (
+            network.stores.e_nom_opt[central_heat.index]
+            .groupby(network.stores.bus)
+            .sum()
+            .reindex(network.buses.index, fill_value=0.0)
+        )
+
+        central_heat_distribution.index = pd.MultiIndex.from_tuples(
+            [
+                (idx, "central_heat_store")
+                for idx in central_heat_distribution.index
+            ]
+        )
+
+        dist.loc[
+            dist.index.get_level_values("carrier") == "central_heat_store"
+        ] = central_heat_distribution
 
     return dist
 
@@ -1060,11 +1117,11 @@ def nodal_gen_dispatch(
     None.
     """
 
-    if osm != False:
+    if osm is not False:
         if set_epsg_network.counter == 0:
             set_epsg_network(network)
         fig, ax = plot_osm(osm["x"], osm["y"], osm["zoom"])
-    elif (osm == False) and cartopy_present:
+    elif (osm is False) and cartopy_present:
         fig, ax = plt.subplots(
             subplot_kw={"projection": ccrs.PlateCarree()}, figsize=(5, 5)
         )
@@ -1184,6 +1241,16 @@ def nodal_production_balance(network, timesteps, scaling=0.00001):
 
     """
 
+    import_gen_from_links(network)
+
+    ac_buses = network.buses[network.buses.carrier == "AC"].index
+    network.generators = network.generators[
+        network.generators.bus.isin(ac_buses)
+    ]
+    network.generators_t.p = network.generators_t.p.loc[
+        :, network.generators_t.p.columns.isin(network.generators.index)
+    ]
+
     gen = (
         mul_weighting(network, network.generators_t.p)
         .groupby(network.generators.bus, axis=1)
@@ -1207,7 +1274,11 @@ def nodal_production_balance(network, timesteps, scaling=0.00001):
     )
 
     bus_sizes = residual_load.abs() * scaling
+    bus_sizes = pd.Series(data=bus_sizes, index=network.buses.index).fillna(0)
 
+    bus_colors = pd.Series(data=bus_colors, index=network.buses.index).fillna(
+        "grey"
+    )
     return bus_sizes, bus_colors
 
 
@@ -1357,6 +1428,7 @@ def storage_soc_sorted(network, filename=None):
     ----------
     network : PyPSA network container
         Holds topology of grid including results from powerflow analysis
+
     filename : path to folder
 
     Returns
@@ -1375,15 +1447,6 @@ def storage_soc_sorted(network, filename=None):
         & (network.storage_units.max_hours == 168)
     ]
 
-    cap_batt = (
-        network.storage_units.max_hours[sbatt]
-        * network.storage_units.p_nom_opt[sbatt]
-    ).sum()
-    cap_hydr = (
-        network.storage_units.max_hours[shydr]
-        * network.storage_units.p_nom_opt[shydr]
-    ).sum()
-
     fig, ax = plt.subplots(1, 1)
 
     if (
@@ -1460,7 +1523,7 @@ def mul_weighting(network, timeseries):
          timeseries considering snapshot_weightings
 
     """
-    return timeseries.mul(network.snapshot_weightings, axis=0)
+    return timeseries.mul(network.snapshot_weightings.generators, axis=0)
 
 
 def calc_ac_loading(network, timesteps):
@@ -1477,7 +1540,7 @@ def calc_ac_loading(network, timesteps):
     Returns
     -------
     pandas.Series
-        ACC line loading in MVA
+        AC line loading in MVA
 
     """
 
@@ -1503,6 +1566,7 @@ def calc_ac_loading(network, timesteps):
 def calc_dc_loading(network, timesteps):
     """Calculates loading of DC-lines
 
+
     Parameters
     ----------
     network : :class:`pypsa.Network
@@ -1527,13 +1591,13 @@ def calc_dc_loading(network, timesteps):
                 & (network.links.length == row["length"])
             ]
         ).empty:
-            l = network.links.index[
+            links_df = network.links.index[
                 (network.links.bus0 == row["bus1"])
                 & (network.links.bus1 == row["bus0"])
                 & (network.links.length == row["length"])
             ]
 
-            network.links.at[i, "linked_to"] = l.values[0]
+            network.links.at[i, "linked_to"] = links_df.values[0]
 
     network.links.linked_to = network.links.linked_to.astype(str)
     # Set p_nom_max and line_loading for one directional links
@@ -1557,14 +1621,16 @@ def calc_dc_loading(network, timesteps):
                 network.links.index == row["linked_to"]
             ].values[0],
         )
-
-    return (
+    dc_load = (
         mul_weighting(network, link_load)
         .loc[network.snapshots[timesteps]]
         .abs()
         .sum()[network.links.index]
         / p_nom_opt_max
     ).dropna()
+    dc_load.loc[network.links.carrier != "DC"] = 0
+
+    return dc_load
 
 
 def plotting_colors(network):
@@ -1607,7 +1673,7 @@ def calc_network_expansion(network, method="abs", ext_min=0.1):
 
     Returns
     -------
-    all_network : :class:`pypsa.Network
+    network : :class:`pypsa.Network
         Whole network including not extended lines
     extension_lines : pandas.Series
         AC-line expansion
@@ -1615,52 +1681,70 @@ def calc_network_expansion(network, method="abs", ext_min=0.1):
         DC-line expansion
 
     """
-    all_network = network.copy()
 
-    network.lines = network.lines[
-        network.lines.s_nom_extendable
+    network_c = network.copy()
+
+    network_c.lines = network_c.lines[
+        network_c.lines.s_nom_extendable
         & (
-            (network.lines.s_nom_opt - network.lines.s_nom_min)
-            / network.lines.s_nom
+            (network_c.lines.s_nom_opt - network_c.lines.s_nom_min)
+            / network_c.lines.s_nom
             >= ext_min
         )
     ]
-    network.links = network.links[
-        network.links.p_nom_extendable
+    network_c.links = network_c.links[
+        network_c.links.p_nom_extendable
+        & (network_c.links.carrier == "DC")
         & (
-            (network.links.p_nom_opt - network.links.p_nom_min)
-            / network.links.p_nom
+            (network_c.links.p_nom_opt - network_c.links.p_nom_min)
+            / network_c.links.p_nom
             >= ext_min
         )
     ]
 
-    for i, row in network.links.iterrows():
-        linked = network.links[
-            (row["bus1"] == network.links.bus0)
-            & (row["bus0"] == network.links.bus1)
+    for i, row in network_c.links.iterrows():
+        linked = network_c.links[
+            (row["bus1"] == network_c.links.bus0)
+            & (row["bus0"] == network_c.links.bus1)
         ]
         if not linked.empty:
             if row["p_nom_opt"] < linked.p_nom_opt.values[0]:
-                network.links.p_nom_opt[i] = linked.p_nom_opt.values[0]
+                network_c.links.p_nom_opt[i] = linked.p_nom_opt.values[0]
 
     if method == "rel":
         extension_lines = (
             100
-            * (network.lines.s_nom_opt - network.lines.s_nom_min)
-            / network.lines.s_nom
+            * (network_c.lines.s_nom_opt - network_c.lines.s_nom_min)
+            / network_c.lines.s_nom
+        )
+
+        extension_links = pd.DataFrame(
+            data=network_c.links, index=network_c.links.index
         )
 
         extension_links = (
             100
-            * (network.links.p_nom_opt - network.links.p_nom_min)
-            / (network.links.p_nom)
+            * (network_c.links.p_nom_opt - network_c.links.p_nom_min)
+            / (network_c.links.p_nom)
         )
+        extension_links = extension_links.fillna(0)
+
     if method == "abs":
-        extension_lines = network.lines.s_nom_opt - network.lines.s_nom_min
+        extension_lines = network_c.lines.s_nom_opt - network_c.lines.s_nom_min
+
+        extension_links = pd.DataFrame(
+            data=network_c.links, index=network_c.links.index
+        )
 
-        extension_links = network.links.p_nom_opt - network.links.p_nom_min
+        extension_links = network_c.links.p_nom_opt - network_c.links.p_nom_min
 
-    return all_network, extension_lines, extension_links
+    extension_lines = pd.Series(
+        data=extension_lines, index=network.lines.index
+    ).fillna(0)
+    extension_links = pd.Series(
+        data=extension_links, index=network.links.index
+    ).fillna(0)
+    return network, extension_lines, extension_links
 
 
 def plot_background_grid(network, ax):
@@ -1678,6 +1762,9 @@ def plot_background_grid(network, ax):
     None.
 
     """
+    link_widths = pd.Series(index=network.links.index, data=0)
+    link_widths.loc[network.links.carrier == "DC"] = 0.3
+
     if cartopy_present:
         network.plot(
             ax=ax,
@@ -1685,10 +1772,11 @@ def plot_background_grid(network, ax):
             link_colors="grey",
             bus_sizes=0,
             line_widths=0.5,
-            link_widths=0.3,  # 0.55,
+            link_widths=link_widths,
             geomap=True,
             projection=ccrs.PlateCarree(),
             color_geomap=True,
+            boundaries=[-2.5, 16, 46.8, 58],
         )
     else:
         network.plot(
@@ -1697,7 +1785,7 @@ def plot_background_grid(network, ax):
             link_colors="grey",
             bus_sizes=0,
             line_widths=0.5,
-            link_widths=0.3,  # 0.55,
+            link_widths=link_widths,
             geomap=False,
         )
 
@@ -1724,26 +1812,26 @@ def demand_side_management(self, buses, snapshots, agg="5h", used=False):
     """
     df = pd.DataFrame(index=self.network.snapshots[snapshots])
 
-    l = self.network.links[
+    link = self.network.links[
         (self.network.links.carrier == "dsm")
         & (self.network.links.bus0.isin(buses))
     ]
     s = self.network.stores[
         (self.network.stores.carrier == "dsm")
-        & (self.network.stores.bus.isin(l.bus1.values))
+        & (self.network.stores.bus.isin(link.bus1.values))
     ]
 
     df["p_min"] = (
-        self.network.links_t.p_min_pu[l.index]
-        .mul(l.p_nom, axis=1)
+        self.network.links_t.p_min_pu[link.index]
+        .mul(link.p_nom, axis=1)
         .sum(axis=1)
         .resample(agg)
         .mean()
         .iloc[snapshots]
     )
     df["p_max"] = (
-        self.network.links_t.p_max_pu[l.index]
-        .mul(l.p_nom, axis=1)
+        self.network.links_t.p_max_pu[link.index]
+        .mul(link.p_nom, axis=1)
         .sum(axis=1)
         .resample(agg)
         .mean()
@@ -1765,7 +1853,7 @@ def demand_side_management(self, buses, snapshots, agg="5h", used=False):
 
     if used:
         df["p"] = (
-            self.network.links_t.p0[l.index]
+            self.network.links_t.p0[link.index]
             .clip(lower=0)
             .sum(axis=1)
             .resample(agg)
@@ -2034,7 +2122,8 @@ def flexibility_usage(
     snapshots : list, optional
         Considered snapshots, if empty all are considered. The default is [].
     buses : list, optional
-        Considered components at AC buses, if empty all are considered. The default is [].
+        Considered components at AC buses, if empty all are considered.
+        The default is [].
     pre_path : str, optional
         State of and where you want to store the figure. The default is None.
 
@@ -2048,9 +2137,6 @@ def flexibility_usage(
     colors["h2_store"] = colors["H2_underground"]
     colors["heat"] = colors["central_heat_store"]
 
-    potential = pd.DataFrame(index=self.network.snapshots[snapshots])
-    used = pd.DataFrame(index=self.network.snapshots[snapshots])
-
     if not buses:
         buses = self.network.buses.index
 
@@ -2228,6 +2314,12 @@ def plot_grid(
     disaggregated=False,
     ext_min=0.1,
     ext_width=False,
+    legend_entries="all",
+    scaling_store_expansion={
+        "H2": 50,
+        "heat": 0.1,
+        "battery": 10,
+    },
 ):
     """Function that plots etrago.network and results for lines and buses
 
@@ -2242,6 +2334,8 @@ def plot_grid(
         * 'expansion_abs': absolute network expansion in MVA
         * 'expansion_rel': network expansion in p.u. of existing capacity
         * 'q_flow_max': maximal reactive flows
+        * 'dlr': energy above nominal capacity
+        * 'grey': plot all lines and DC links grey colored
 
     bus_sizes : float, optional
         Size of buses. The default is 0.001.
@@ -2252,10 +2346,15 @@ def plot_grid(
         * 'nodal_production_balance': net producer/consumer in selected timeteps
         * 'storage_expansion': storage expansion per bus and technology
         * 'storage_distribution': installed storage units per bus
+        * 'h2_battery_storage_expansion': storage expansion per bus and
+           technology for underground and overground H2 and batteries.
         * 'gen_dist': dispatch per carrier in selected timesteps
+        * 'PowerToH2': location and sizes of electrolizers
+        * 'flexibility_usage': use of DSM and BEV charger
 
     timesteps : array, optional
-        Timesteps consideredd in time depended plots. The default is range(2).
+        Timesteps consideredd in time depended plots. The default
+        is range(2).
     osm : bool or dict, e.g. {'x': [1,20], 'y': [47, 56], 'zoom' : 6}
         If not False, osm is set as background
         with the following settings as dict:
@@ -2273,9 +2372,15 @@ def plot_grid(
     ext_min: float
         Choose minimum relative line extension shown in plot in p.u..
     ext_width: float or bool
-        Choose if line_width respects line extension. Turn off with 'False' or
-        set linear factor to decremise extension line_width.
+        Choose if line_width respects line extension. Turn off with
+        'False' or set linear factor to decremise extension line_width.
         The default is False.
+    legend_entries : list, optional
+        Set the legends for buses to be plotted. The default is 'all'.
+    scaling_store_expansion : dict, optional
+        Set scaling values to be used per technology for the plots
+        storage_expansion and h2_battery_storage_expansion. The default is
+        {"H2": 50, "heat": 0.1, "battery": 10}
 
     Returns
     -------
@@ -2293,16 +2398,17 @@ def plot_grid(
 
     # Set default values
     flow = None
+    title = ""
     line_widths = 2
     link_widths = 0
 
     # Plot osm map in background
-    if osm != False:
+    if osm is not False:
         if network.srid == 4326:
             set_epsg_network(network)
         fig, ax = plot_osm(osm["x"], osm["y"], osm["zoom"])
 
-    elif (osm == False) and cartopy_present:
+    elif (osm is False) and cartopy_present:
         fig, ax = plt.subplots(
             subplot_kw={"projection": ccrs.PlateCarree()}, figsize=(5, 5)
         )
@@ -2310,22 +2416,28 @@ def plot_grid(
     else:
         fig, ax = plt.subplots(figsize=(5, 5))
 
+    fig.set_tight_layout(True)
+
     # Set line colors
     if line_colors == "line_loading":
-        title = (
-            "Mean loading from "
-            + str(network.snapshots[timesteps[0]])
-            + " to "
-            + str(network.snapshots[timesteps[-1]])
-        )
-        rep_snapshots = network.snapshot_weightings[
+        title = "Mean line loading"
+        rep_snapshots = network.snapshot_weightings["objective"][
             network.snapshots[timesteps]
         ].sum()
         line_colors = calc_ac_loading(network, timesteps).abs() / rep_snapshots
         link_colors = calc_dc_loading(network, timesteps).abs() / rep_snapshots
+        if ext_width is not False:
+            link_widths = link_colors.apply(
+                lambda x: 10 + (x / ext_width) if x != 0 else 0
+            )
+            line_widths = 10 + (line_colors / ext_width)
+        else:
+            link_widths = link_colors.apply(lambda x: 10 if x != 0 else 0)
+            line_widths = 10
         label = "line loading in p.u."
+        plot_background_grid(network, ax)
         # Only active flow direction is displayed!
-        flow = pd.Series(index=network.branches().index, dtype="float64")
+        flow = pd.Series(1, index=network.branches().index, dtype="float64")
         flow.iloc[flow.index.get_level_values("component") == "Line"] = (
             mul_weighting(network, network.lines_t.p0)
             .loc[network.snapshots[timesteps]]
@@ -2333,24 +2445,58 @@ def plot_grid(
             / network.lines.s_nom
             / rep_snapshots
         ).values
-        flow.iloc[flow.index.get_level_values("component") == "Link"] = (
-            calc_dc_loading(network, timesteps) / rep_snapshots
-        ).values
+
+        dc_loading = calc_dc_loading(network, timesteps) / rep_snapshots
+        dc_loading.index = pd.MultiIndex.from_tuples(
+            [("Link", name) for name in dc_loading.index],
+            names=["component", "name"],
+        )
+        flow.loc["Link", :] = dc_loading
+
+        flow = flow[
+            (flow.index.get_level_values("component") == "Line")
+            | (
+                flow.index.isin(
+                    link_widths[
+                        link_widths.index.isin(
+                            network.links[network.links.carrier == "DC"].index
+                        )
+                    ].index,
+                    level=1,
+                )
+            )
+        ]
+        flow[flow < 0] = -1
+        flow[flow > 0] = 1
+
     elif line_colors == "v_nom":
         title = "Voltage levels"
         label = "v_nom in kV"
         line_colors = network.lines.v_nom
-        link_colors = network.links.v_nom
+        link_colors = pd.Series(data=0, index=network.links.index)
+        plot_background_grid(network, ax)
     elif line_colors == "expansion_abs":
         title = "Network expansion"
-        label = "network expansion in MW"
+        label = "network expansion in GVA"
         all_network, line_colors, link_colors = calc_network_expansion(
             network, method="abs", ext_min=ext_min
         )
         plot_background_grid(all_network, ax)
-        if ext_width != False:
-            line_widths = 0.5 + (line_colors / ext_width)
-            link_widths = 0.5 + (link_colors / ext_width)
+
+        if ext_width is not False:
+            line_widths = line_colors / ext_width
+            link_widths = link_colors.apply(
+                lambda x: x / ext_width if x != 0 else 0
+            )
+        else:
+            dc_link = network.links.index[network.links.carrier == "DC"]
+            link_widths = pd.Series(0, index=network.links.index)
+            link_widths.loc[dc_link] = 1.5
+            line_widths = line_colors.apply(lambda x: 1.5 if x != 0 else 0)
+
+        link_colors = link_colors.mul(1e-3)
+        line_colors = line_colors.mul(1e-3)
+
     elif line_colors == "expansion_rel":
         title = "Network expansion"
         label = "network expansion in %"
@@ -2358,20 +2504,71 @@ def plot_grid(
             network, method="rel", ext_min=ext_min
         )
         plot_background_grid(all_network, ax)
-        if ext_width != False:
+        if ext_width is not False:
             line_widths = 0.5 + (line_colors / ext_width)
-            link_widths = 0.5 + (link_colors / ext_width)
+            link_widths = link_colors.apply(
+                lambda x: 0.5 + x / ext_width if x != 0 else 0
+            )
+        else:
+            dc_link = network.links.index[network.links.carrier == "DC"]
+            link_widths = pd.Series(0, index=network.links.index)
+            link_widths.loc[dc_link] = 2
+            line_widths = line_colors.apply(lambda x: 1.5 if x != 0 else 0)
     elif line_colors == "q_flow_max":
-        title = "Maximmal reactive power flows"
-        label = "flow in Mvar"
+        title = "Maximum reactive power flows"
+        label = "flow in pu"
         line_colors = abs(
             network.lines_t.q0.abs().max() / (network.lines.s_nom)
         )
+        if ext_width is not False:
+            line_widths = 0.5 + (line_colors / ext_width)
+        link_colors = pd.Series(data=0, index=network.links.index)
+        plot_background_grid(network, ax)
+    elif line_colors == "dlr":
+        title = "Dynamic line rating"
+        label = "TWh above nominal capacity"
+        plot_background_grid(network, ax)
+
+        # calc min capacity per line in the given period: Since lines with
+        # different original voltage level could be aggregated during the
+        # clustering, the security factors can be values in between the values
+        # provided in the args for branch_capacity_factor.
+        network.lines.s_max_pu = network.lines_t.s_max_pu.min()
+        line_loading = network.lines_t.p0.mul(
+            1 / (network.lines.s_nom_opt * network.lines.s_max_pu)
+        ).abs()
+        line_loading = line_loading.iloc[timesteps, :]
+        # keep only the capacity allowed by dlr
+        line_loading = line_loading - 1
+        dlr_usage = (
+            line_loading[line_loading > 0]
+            .fillna(0)
+            .mul(network.snapshot_weightings.generators, axis=0)
+            .sum()
+        )
+        dlr_usage = (
+            dlr_usage * network.lines.s_nom * network.lines.s_max_pu / 1000000
+        )
+        dlr_usage = dlr_usage.round(decimals=0)
+        line_colors = dlr_usage
+        if ext_width is not False:
+            line_widths = 0.2 + (line_colors / ext_width)
         link_colors = pd.Series(data=0, index=network.links.index)
+
+    elif line_colors == "grey":
+        title = ""
+        label = ""
+        line_colors = "grey"
+        link_colors = "grey"
+        plot_background_grid(network, ax)
+        link_widths = 0
+        line_widths = 0
+
     else:
         logger.warning("line_color {} undefined".format(line_colors))
 
     # Set bus colors
+    bus_legend = False
 
     if bus_colors == "nodal_production_balance":
         bus_scaling = bus_sizes
@@ -2381,10 +2578,50 @@ def plot_grid(
         bus_legend = "Nodal production balance"
         bus_unit = "TWh"
     elif bus_colors == "storage_expansion":
+        if not isinstance(scaling_store_expansion, dict):
+            raise Exception(
+                """To plot storage_expansion, the argument\
+            scaling_store_expansion must be a dictionary like:
+                            {"H2": 50,
+                            "heat": 0.1,
+                            "battery": 10}"""
+            )
         bus_scaling = bus_sizes
         bus_sizes = bus_scaling * calc_storage_expansion_per_bus(network)
+        for store_carrier in ["H2", "heat", "battery"]:
+            bus_sizes[
+                bus_sizes.index.get_level_values("carrier").str.contains(
+                    store_carrier
+                )
+            ] *= scaling_store_expansion[store_carrier]
         bus_legend = "Storage expansion"
         bus_unit = "GW"
+    elif bus_colors == "h2_battery_storage_expansion":
+        bus_scaling = bus_sizes
+        bus_sizes = bus_scaling * calc_storage_expansion_per_bus(
+            network, carriers=["battery", "H2_overground", "H2_underground"]
+        )
+        if (
+            ("battery" not in scaling_store_expansion.keys())
+            | ("H2_overground" not in scaling_store_expansion.keys())
+            | ("H2_underground" not in scaling_store_expansion.keys())
+        ):
+            raise Exception(
+                """To plot h2_battery_storage_expansion, the argument\
+            scaling_store_expansion must be a dictionary like:
+                            {"H2_overground": 1,
+                             "H2_underground": 1,
+                             "battery": 1,}"""
+            )
+
+        for store_carrier in ["battery", "H2_overground", "H2_underground"]:
+            bus_sizes[
+                bus_sizes.index.get_level_values("carrier").str.contains(
+                    store_carrier
+                )
+            ] *= scaling_store_expansion[store_carrier]
+        bus_legend = "Battery and H2 storage expansion"
+        bus_unit = "GW"
     elif bus_colors == "storage_distribution":
         bus_scaling = bus_sizes
         bus_sizes = (
@@ -2398,6 +2635,40 @@ def plot_grid(
         bus_sizes = bus_scaling * calc_dispatch_per_carrier(network, timesteps)
         bus_legend = "Dispatch"
         bus_unit = "TW"
+    elif bus_colors == "flexibility_usage":
+        bus_scaling = bus_sizes
+        flex_links = network.links[
+            network.links.carrier.isin(
+                [
+                    "dsm",
+                    "BEV charger",
+                ]
+            )
+        ]
+        flex_links["p0_sum"] = (
+            network.links_t.p0[flex_links.index]
+            .mul(network.snapshot_weightings.generators, axis=0)
+            .abs()
+            .sum()
+        )
+        flex_links["p1_sum"] = (
+            network.links_t.p1[flex_links.index]
+            .mul(network.snapshot_weightings.generators, axis=0)
+            .sum()
+        )
+        bus_sizes = (
+            bus_scaling * flex_links.groupby(["bus0", "carrier"]).p0_sum.sum()
+        )
+        bus_unit = "TWh"
+        bus_legend = "flexibility_usage"
+    elif bus_colors == "h2_storage_expansion":
+        bus_scaling = bus_sizes
+        bus_sizes = bus_scaling * calc_storage_expansion_per_bus(network)
+        bus_sizes = bus_sizes.reset_index()
+        bus_sizes = bus_sizes[bus_sizes.carrier.str.contains("H2")]
+        bus_sizes.set_index(["bus", "carrier"], inplace=True)
+        bus_legend = "Storage expansion"
+        bus_unit = "GW"
     elif (
         bus_colors == "PowerToH2"
     ):  # PowerToH2 plots p_nom_opt of links with carrier=power to H2
@@ -2409,10 +2680,19 @@ def plot_grid(
             .sum()
             .p_nom_opt
         )
+        if len(bus_sizes) == 0:
+            print("There is no PowerToH2 to plot")
         bus_colors = coloring()["power_to_H2"]
         bus_legend = "PowerToH2"
         bus_unit = "TW"
-
+    elif bus_colors == "grey":
+        bus_scaling = bus_sizes
+        bus_sizes = pd.Series(
+            data=network.buses.carrier, index=network.buses.index
+        )
+        bus_sizes[bus_sizes != "AC"] = 0
+        bus_sizes[bus_sizes == "AC"] = 1 * bus_scaling
+        bus_scaling = bus_sizes
     else:
         logger.warning("bus_color {} undefined".format(bus_colors))
 
@@ -2425,12 +2705,13 @@ def plot_grid(
             bus_sizes=bus_sizes,
             bus_colors=bus_colors,
             line_widths=line_widths,
-            link_widths=0,  # link_widths,
+            link_widths=link_widths,
             flow=flow,
             title=title,
             geomap=False,
             projection=ccrs.PlateCarree(),
             color_geomap=True,
+            boundaries=[-2.5, 16, 46.8, 58],
         )
     else:
         ll = network.plot(
@@ -2441,24 +2722,63 @@ def plot_grid(
             bus_sizes=bus_sizes,
             bus_colors=bus_colors,
             line_widths=line_widths,
-            link_widths=0,  # link_widths,
+            link_widths=link_widths,
             flow=flow,
             title=title,
             geomap=False,
+            boundaries=[-2.5, 16, 46.8, 58],
         )
+    l3 = None
 
     # legends for bus sizes and colors
-    if type(bus_sizes) != float:
-        handles = make_legend_circles_for(
-            [bus_sizes.min(), bus_sizes.max()], scale=1, facecolor="gray"
-        )
-        labels = [
-            ("{} " + bus_unit).format(s)
-            for s in (
-                round(bus_sizes.min() / bus_scaling / 1000, 0),
-                round(bus_sizes.max() / bus_scaling / 1000, 0),
+    if bus_legend:
+        handles = []
+        labels = []
+        if scaling_store_expansion:
+            if not isinstance(legend_entries, list):
+                if bus_legend == "Storage expansion":
+                    legend_entries = list(scaling_store_expansion.keys())
+                if bus_legend == "Battery and H2 storage expansion":
+                    legend_entries = [
+                        "battery",
+                        "H2_overground",
+                        "H2_underground",
+                    ]
+            for i in legend_entries:
+                try:
+                    max_value = bus_sizes[
+                        bus_sizes.index.get_level_values(
+                            "carrier"
+                        ).str.contains(i)
+                    ].max()
+                except KeyError:
+                    max_value = bus_sizes.max()
+                handles.append(
+                    make_legend_circles_for(
+                        [max_value],
+                        scale=1,
+                        facecolor=network.carriers.color[i],
+                    )[0]
+                )
+                labels.append(
+                    f"""
+                    {round(max_value/bus_scaling/scaling_store_expansion[i]/
+                           1000, 0).astype(int)} {bus_unit} """
+                    + i
+                )
+        else:
+            if len(bus_sizes) > 0:
+                max_value = bus_sizes.max()
+            else:
+                max_value = 0
+            labels.append(f"{round(max_value / bus_scaling /1000, 0)} GWh ")
+            handles.append(
+                make_legend_circles_for(
+                    [max_value],
+                    scale=1,
+                    facecolor="grey",
+                )[0]
             )
-        ]
 
         l2 = ax.legend(
             handles,
@@ -2469,61 +2789,91 @@ def plot_grid(
             framealpha=1.0,
             title=bus_legend,
             handler_map=make_handler_map_to_scale_circles_as_in(ax),
+            prop={"size": 8},
         )
         ax.add_artist(l2)
 
-        handles = []
-        if bus_legend == "Nodal production balance":
-            positive = mpatches.Patch(color="green", label="generation")
-            negative = mpatches.Patch(color="red", label="consumption")
-            handles = [positive, negative]
-        else:
-            for i in network.carriers.color.index:
-                patch = mpatches.Patch(
-                    color=network.carriers.color[i], label=i
-                )
-                handles.append(patch)
-
-        l3 = plt.legend(
-            handles=handles, loc="upper left", ncol=3, bbox_to_anchor=(-0.1, 0)
-        )
-        ax.add_artist(l3)
+        plt.setp(l2.get_title(), fontsize="9")
+
+        if not scaling_store_expansion:
+            handles = []
+            if bus_legend == "Nodal production balance":
+                positive = mpatches.Patch(color="green", label="generation")
+                negative = mpatches.Patch(color="red", label="consumption")
+                handles = [positive, negative]
+
+            elif bus_legend == "PowerToH2":
+                pth = mpatches.Patch(color="cyan", label="PowerToH2")
+                handles = [pth]
+            elif legend_entries != "all":
+                for i in legend_entries:
+                    patch = mpatches.Patch(
+                        color=network.carriers.color[i], label=i
+                    )
+                    handles.append(patch)
+            else:
+                for i in bus_sizes.index.get_level_values("carrier").unique():
+                    patch = mpatches.Patch(
+                        color=network.carriers.color[i], label=i
+                    )
+                    handles.append(patch)
+
+            l3 = plt.legend(
+                handles=handles,
+                loc="upper left",
+                ncol=2,
+                bbox_to_anchor=(0, 0),
+            )
+            ax.add_artist(l3)
+
+    if type(line_colors) != str:
+        # Set fixed boundaries if selected in parameters
+        if not boundaries:
+            boundaries = [
+                min(round(line_colors.min(), 1), round(link_colors.min(), 1)),
+                max(round(line_colors.max()), round(link_colors.max())),
+            ]
 
-    # Set fixed boundaries if selected in parameters
-    if not boundaries:
-        boundaries = [
-            min(line_colors.min(), link_colors.min()),
-            max(line_colors.max(), link_colors.max()),
+        # Create ticks for legend
+        v = [
+            round(x, 1) for x in np.linspace(boundaries[0], boundaries[1], 101)
         ]
+        for l_collection in ll:
+            l_collection.set_clim(boundaries[0], boundaries[1])
 
-    # Create ticks for legend
-    v = np.linspace(boundaries[0], boundaries[1], 101)
-
-    # colorbar for line heatmap
-    cb = plt.colorbar(
-        ll[1], boundaries=v, ticks=v[0:101:10], fraction=0.046, pad=0.04
-    )
-    # Set legend label
-    cb.set_label(label)
+        # colorbar for line heatmap
+        cb = plt.colorbar(
+            ll[1],
+            values=v,
+            ticks=v[0:101:10],
+            fraction=0.028,
+            pad=0.04,
+        )
+        # Set legend label
+        cb.set_label(label)
 
     # Show plot or save to file
     if filename is None:
-        if type(bus_sizes) != float:
+        if not isinstance(bus_sizes, (pd.Series, float)):
             logger.warning("Legend of bus sizes will change when zooming")
+        plt.tight_layout()
         plt.show()
     else:
         from matplotlib import pylab
 
-        pylab.savefig(filename, dpi=300, bbox_inches="tight")
+        if l3 is None:
+            pylab.savefig(filename, dpi=300, bbox_inches="tight")
+        else:
+            pylab.savefig(
+                filename, dpi=300, bbox_inches="tight", bbox_extra_artists=[l3]
+            )
         plt.close()
 
 
 set_epsg_network.counter = 0
 
-### the following functions are copied from pypsa-eur-sec ###
-### see here: https://github.com/PyPSA/pypsa-eur-sec/blob/master/scripts/plot_network.py
-from matplotlib.legend_handler import HandlerPatch
-from matplotlib.patches import Circle, Ellipse
+# the following functions are copied from pypsa-eur-sec. see:
+# https://github.com/PyPSA/pypsa-eur-sec/blob/master/scripts/plot_network.py
 
 
 def make_handler_map_to_scale_circles_as_in(ax, dont_resize_actively=False):
@@ -2731,7 +3081,8 @@ def plot_gas_generation(
     self : :class:`Etrago
         Overall container of Etrago
     t_resolution : str, optional
-        sets the resampling rate of timeseries data to allow for smoother line plots
+        sets the resampling rate of timeseries data to allow for smoother
+        line plots
     save_path : bool, optional
         Path to save the generated plot. The default is False.
 
@@ -2796,10 +3147,12 @@ def plot_gas_summary(self, t_resolution="20H", stacked=True, save_path=False):
     self : :class:`Etrago
         Overall container of Etrago
     t_resolution : str, optional
-        sets the resampling rate of timeseries data to allow for smoother line plots
+        sets the resampling rate of timeseries data to allow for smoother
+        line plots
     stacked : bool, optional
-        If True all TS data will be shown as stacked area plot. Total gas generation
-        will then also be plotted to check for matching demand and generation.
+        If True all TS data will be shown as stacked area plot. Total gas
+        generation will then also be plotted to check for matching demand and
+        generation.
     save_path : bool, optional
         Path to save the generated plot. The default is False.
 
@@ -2950,7 +3303,8 @@ def plot_h2_generation(self, t_resolution="20H", save_path=False):
     self : :class:`Etrago
         Overall container of Etrago
     t_resolution : str, optional
-        sets the resampling rate of timeseries data to allow for smoother line plots
+        sets the resampling rate of timeseries data to allow for smoother
+        line plots
     save_path : bool, optional
         Path to save the generated plot. The default is False.
 
@@ -3008,10 +3362,12 @@ def plot_h2_summary(self, t_resolution="20H", stacked=True, save_path=False):
     self : :class:`Etrago
         Overall container of Etrago
     t_resolution : str, optional
-        sets the resampling rate of timeseries data to allow for smoother line plots
+        sets the resampling rate of timeseries data to allow for smoother
+        line plots
     stacked : bool, optional
-        If True all TS data will be shown as stacked area plot. Total H2 generation
-        will then also be plotted to check for matching demand and generation.
+        If True all TS data will be shown as stacked area plot. Total H2
+        generation will then also be plotted to check for matching demand and
+        generation.
     save_path : bool, optional
         Path to save the generated plot. The default is False.
 
@@ -3020,7 +3376,6 @@ def plot_h2_summary(self, t_resolution="20H", stacked=True, save_path=False):
     None.
 
     """
-    colors = coloring()
 
     rel_h2_links = ["H2_feedin", "H2_to_CH4", "H2_to_power"]
     rel_h2_loads = ["H2_for_industry", "H2_hgv_load"]
@@ -3133,7 +3488,8 @@ def plot_heat_loads(self, t_resolution="20H", save_path=False):
     self : :class:`Etrago
         Overall container of Etrago
     t_resolution : str, optional
-        sets the resampling rate of timeseries data to allow for smoother line plots
+        sets the resampling rate of timeseries data to allow for smoother
+        line plots
     save_path : bool, optional
         Path to save the generated plot. The default is False.
 
@@ -3182,10 +3538,12 @@ def plot_heat_summary(self, t_resolution="20H", stacked=True, save_path=False):
     self : :class:`Etrago
         Overall container of Etrago
     t_resolution : str, optional
-        sets the resampling rate of timeseries data to allow for smoother line plots
+        sets the resampling rate of timeseries data to allow for smoother
+        line plots
     stacked : bool, optional
-        If True all TS data will be shown as stacked area plot. Total heat demand
-        will then also be plotted to check for matching generation and demand.
+        If True all TS data will be shown as stacked area plot. Total heat
+        demand will then also be plotted to check for matching generation and
+        demand.
     save_path : bool, optional
         Path to save the generated plot. The default is False.
 
@@ -3257,7 +3615,7 @@ def plot_heat_summary(self, t_resolution="20H", stacked=True, save_path=False):
         / 1e3
     )
 
-    if stacked == True:
+    if stacked:
         data = pd.DataFrame(-(data.sum(axis=1)))
         data = data.rename(columns={0: heat_gen_techs[0]})