diff --git a/etrago/analyze/plot.py b/etrago/analyze/plot.py
index 07dd75d6..ce7488e8 100644
--- a/etrago/analyze/plot.py
+++ b/etrago/analyze/plot.py
@@ -51,6 +51,7 @@
     import tilemapbase
 
     from etrago.execute import import_gen_from_links
+    from etrago.tools.utilities import find_buses_area
 
 __copyright__ = (
     "Flensburg University of Applied Sciences, "
@@ -2474,6 +2475,8 @@ def plot_grid(
         * 'q_flow_max': maximal reactive flows
         * 'dlr': energy above nominal capacity
         * 'grey': plot all lines and DC links grey colored
+        * 'interest_area': plot all AC buses inside the interest area larger than buses
+          outside the interest area
 
     bus_sizes : float, optional
         Size of buses. The default is 0.001.
@@ -2907,6 +2910,20 @@ def plot_grid(
         bus_sizes[bus_sizes != "AC"] = 0
         bus_sizes[bus_sizes == "AC"] = 1 * bus_scaling
         bus_scaling = bus_sizes
+    elif bus_colors == "interest_area":
+        bus_scaling = bus_sizes
+        # only plot AC buses
+        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
+        # only plot buses inside interest area
+        buses_interest_area = find_buses_area(self, "AC")
+        buses_outside = [_ for _ in bus_sizes.index if _ not in buses_interest_area]
+        bus_sizes.loc[buses_outside] = bus_sizes.loc[buses_outside] * 0.3
+        bus_scaling = bus_sizes
+        bus_colors = coloring()["AC"]
     else:
         logger.warning("bus_color {} undefined".format(bus_colors))
 
diff --git a/etrago/appl.py b/etrago/appl.py
index 9f4eca97..224e25ed 100644
--- a/etrago/appl.py
+++ b/etrago/appl.py
@@ -112,6 +112,7 @@
     "extra_functionality": {},  # Choose function name or {}
     # Spatial Complexity:
     "delete_dispensable_ac_buses": True,  # bool. Find and delete expendable buses
+    "interest_area": False,  # False, path to shapefile or list of nuts names of the area that is excluded from the clustering. By default the buses inside remain the same, but the parameter "n_cluster_interest_area" inside "network clustering" defines if it should be clustered to a certain number of buses.
     "network_clustering_ehv": {
         "active": False,  # choose if clustering of HV buses to EHV buses is activated
         "busmap": False,  # False or path to stored busmap
@@ -119,8 +120,9 @@
     "network_clustering": {
         "active": True,  # choose if clustering is activated
         "method": "kmedoids-dijkstra",  # choose clustering method: kmeans or kmedoids-dijkstra
-        "n_clusters_AC": 30,  # total number of resulting AC nodes (DE+foreign)
+        "n_clusters_AC": 30,  # total number of resulting AC nodes (DE+foreign-interest_area)
         "cluster_foreign_AC": False,  # take foreign AC buses into account, True or False
+        "n_cluster_interest_area": False, # False or number of buses.
         "method_gas": "kmedoids-dijkstra",  # choose clustering method: kmeans or kmedoids-dijkstra
         "n_clusters_gas": 14,  # total number of resulting CH4 nodes (DE+foreign)
         "cluster_foreign_gas": False,  # take foreign CH4 buses into account, True or False
@@ -383,6 +385,16 @@ def run_etrago(args, json_path):
         connected lines are merged. This reduces the spatial complexity without
         losing any accuracy.
         Default: True.
+
+    interest_area: False, list, string
+        Area of especial interest that will be not clustered, except when
+        n_cluster_interest_area is provided. It is by default set to false.
+        When an interest_area is provided, the given value for n_clusters_AC
+        will mean the total of AC buses outside the area.The area can be
+        provided in two ways: list of nuts names e.G.
+        ["Cuxhaven", "Bremerhaven", "Bremen"] or a string with a path to a
+        shape file (.shp).
+        
     network_clustering_ehv : dict
         Choose if you want to apply an extra high voltage clustering to the
         electrical network.
@@ -414,12 +426,14 @@ def run_etrago(args, json_path):
             * "kmeans": considers geographical locations of buses
             * "kmedoids-dijkstra":  considers electrical distances between
             buses
-
             Default: "kmedoids-dijkstra".
-        * "n_clusters_AC" : int
+        * "n_clusters_AC" : int, False
             Defines total number of resulting AC nodes including DE and foreign
             nodes if `cluster_foreign_AC` is set to True, otherwise only DE
-            nodes.
+            nodes. When using the interest_area parameter, n_clusters_AC could
+            be set to False, which means that only the buses inside the 
+            provided area are clustered.
+
             Default: 30.
         * "cluster_foreign_AC" : bool
             If set to False, the AC buses outside Germany are not clustered
@@ -428,13 +442,18 @@ def run_etrago(args, json_path):
             as well and included in number of clusters specified through
             ``'n_clusters_AC'``.
             Default: False.
+
+        * "n_cluster_interest_area": False, int
+            Number of buses to cluster all the electrical buses in the area
+            of interest. Method provided in the arg "method" is used. If
+            it is set to False, the area of interest is not clustered.
+            Default: False.
         * "method_gas" : str
             Method used for gas clustering. You can choose between two
             clustering methods:
             * "kmeans": considers geographical locations of buses
             * "kmedoids-dijkstra":  considers 'electrical' distances between
             buses
-
             Default: "kmedoids-dijkstra".
         * "n_clusters_gas" : int
             Defines total number of resulting CH4 nodes including DE and
@@ -694,17 +713,6 @@ def run_etrago(args, json_path):
     # skip snapshots
     etrago.skip_snapshots()
 
-    # Temporary drop DLR as it is currently not working with sclopf
-    if etrago.args["method"]["type"] != "lopf":
-        etrago.network.lines_t.s_max_pu = pd.DataFrame(
-            index=etrago.network.snapshots,
-            columns=etrago.network.lines.index,
-            data=1.0,
-        )
-
-    etrago.network.lines.loc[etrago.network.lines.r == 0.0, "r"] = 10
-
-    # start linear optimal powerflow calculations
     etrago.optimize()
 
     # conduct lopf with full complex timeseries for dispatch disaggregation
diff --git a/etrago/args.json b/etrago/args.json
index 63c1fc1d..45977fef 100644
--- a/etrago/args.json
+++ b/etrago/args.json
@@ -61,6 +61,7 @@
     "generator_noise": 789456,
     "extra_functionality": {},
     "delete_dispensable_ac_buses": true,
+    "interest_area": false,
     "network_clustering_ehv": {
         "active": false,
         "busmap": false
@@ -70,6 +71,7 @@
         "method": "kmedoids-dijkstra",
         "n_clusters_AC": 30,
         "cluster_foreign_AC": false,
+        "n_cluster_interest_area": false,
         "method_gas": "kmedoids-dijkstra",
         "n_clusters_gas": 17,
         "cluster_foreign_gas": false,
diff --git a/etrago/cluster/electrical.py b/etrago/cluster/electrical.py
index 25102865..35f6dc60 100755
--- a/etrago/cluster/electrical.py
+++ b/etrago/cluster/electrical.py
@@ -48,7 +48,10 @@
         strategies_lines,
         strategies_one_ports,
     )
-    from etrago.tools.utilities import set_control_strategies
+    from etrago.tools.utilities import (
+        set_control_strategies,
+        select_elec_network,
+    )
 
     logger = logging.getLogger(__name__)
 
@@ -196,6 +199,7 @@ def find_de_closest(network, bus_ne):
         & (network.buses.carrier == "AC")
         & (network.buses.v_nom > 110)
     ].index
+
     busmap3 = pd.DataFrame(columns=["elec_bus", "carrier", "cluster"])
     for bus_ne in network2.buses.index:
         carry = network2.buses.loc[bus_ne, "carrier"]
@@ -278,7 +282,6 @@ def cluster_on_extra_high_voltage(etrago, busmap, with_time=True):
     busmap : dict
         Maps old bus_ids to new bus_ids including all sectors.
     """
-
     network_c = Network()
 
     network, busmap = adjust_no_electric_network(
@@ -298,7 +301,7 @@ def cluster_on_extra_high_voltage(etrago, busmap, with_time=True):
 
     # keep attached lines
     lines = network.lines.copy()
-    mask = lines.bus0.isin(buses.index)
+    mask = lines.bus0.isin(buses.index) & lines.bus1.isin(buses.index)
     lines = lines.loc[mask, :]
 
     # keep attached transformer
@@ -471,132 +474,6 @@ def ehv_clustering(self):
         logger.info("Network clustered to EHV-grid")
 
 
-def select_elec_network(etrago, apply_on="grid_model"):
-    """
-    Selects the electric network based on the clustering settings specified
-    in the Etrago object.
-
-    Parameters
-    ----------
-    etrago : Etrago
-        An instance of the Etrago class
-    apply_on: str
-        gives information about the objective of the output network. If
-        "grid_model" is provided, the value assigned in the args for
-        ["network_clustering"]["cluster_foreign_AC""] will define if the
-        foreign buses will be included in the network. if "market_model" is
-        provided, foreign buses will be always included.
-
-    Returns
-    -------
-    Tuple containing:
-        elec_network : pypsa.Network
-            Contains the electric network
-        n_clusters : int
-            number of clusters used in the clustering process.
-    """
-    if apply_on == "grid_model":
-        elec_network = etrago.network.copy()
-    elif apply_on == "market_model":
-        elec_network = etrago.network_tsa.copy()
-    else:
-        logger.warning(
-            """Parameter apply_on must be either 'grid_model' or 'market_model'
-            """
-        )
-    settings = etrago.args["network_clustering"]
-
-    if apply_on == "grid_model":
-        include_foreign = settings["cluster_foreign_AC"]
-    elif apply_on == "market_model":
-        include_foreign = True
-    else:
-        raise ValueError(
-            """Parameter apply_on must be either 'grid_model' or 'market_model'
-            """
-        )
-
-    if include_foreign:
-        elec_network.buses = elec_network.buses[
-            elec_network.buses.carrier == "AC"
-        ]
-        elec_network.links = elec_network.links[
-            (elec_network.links.carrier == "AC")
-            | (elec_network.links.carrier == "DC")
-        ]
-        n_clusters = settings["n_clusters_AC"]
-    else:
-        AC_filter = elec_network.buses.carrier.values == "AC"
-
-        foreign_buses = elec_network.buses[
-            (elec_network.buses.country != "DE")
-            & (elec_network.buses.carrier == "AC")
-        ]
-
-        num_neighboring_country = len(
-            foreign_buses[foreign_buses.index.isin(elec_network.loads.bus)]
-        )
-
-        elec_network.buses = elec_network.buses[
-            AC_filter & (elec_network.buses.country.values == "DE")
-        ]
-        n_clusters = settings["n_clusters_AC"] - num_neighboring_country
-
-    # Dealing with generators
-    elec_network.generators = elec_network.generators[
-        elec_network.generators.bus.isin(elec_network.buses.index)
-    ]
-
-    for attr in elec_network.generators_t:
-        elec_network.generators_t[attr] = elec_network.generators_t[attr].loc[
-            :,
-            elec_network.generators_t[attr].columns.isin(
-                elec_network.generators.index
-            ),
-        ]
-
-    # Dealing with loads
-    elec_network.loads = elec_network.loads[
-        elec_network.loads.bus.isin(elec_network.buses.index)
-    ]
-
-    for attr in elec_network.loads_t:
-        elec_network.loads_t[attr] = elec_network.loads_t[attr].loc[
-            :,
-            elec_network.loads_t[attr].columns.isin(elec_network.loads.index),
-        ]
-
-    # Dealing with storage_units
-    elec_network.storage_units = elec_network.storage_units[
-        elec_network.storage_units.bus.isin(elec_network.buses.index)
-    ]
-
-    for attr in elec_network.storage_units_t:
-        elec_network.storage_units_t[attr] = elec_network.storage_units_t[
-            attr
-        ].loc[
-            :,
-            elec_network.storage_units_t[attr].columns.isin(
-                elec_network.storage_units.index
-            ),
-        ]
-
-    # Dealing with stores
-    elec_network.stores = elec_network.stores[
-        elec_network.stores.bus.isin(elec_network.buses.index)
-    ]
-
-    for attr in elec_network.stores_t:
-        elec_network.stores_t[attr] = elec_network.stores_t[attr].loc[
-            :,
-            elec_network.stores_t[attr].columns.isin(
-                elec_network.stores.index
-            ),
-        ]
-
-    return elec_network, n_clusters
-
-
 def unify_foreign_buses(etrago):
     """
     Unifies foreign AC buses into clusters using the k-medoids algorithm with
@@ -787,7 +664,9 @@ def preprocessing(etrago, apply_on="grid_model"):
     else:
         busmap_foreign = pd.Series(name="foreign", dtype=str)
 
-    network_elec, n_clusters = select_elec_network(etrago, apply_on=apply_on)
+    network_elec, n_clusters, network_area = select_elec_network(
+        etrago, apply_on=apply_on
+    )
 
     if settings["method"] == "kmedoids-dijkstra":
         lines_col = network_elec.lines.columns
@@ -817,8 +696,16 @@ def preprocessing(etrago, apply_on="grid_model"):
         weight.index = weight.index.astype(str)
     else:
         weight = weighting_for_scenario(network=network_elec, save=False)
-
-    return network_elec, weight, n_clusters, busmap_foreign
+        weight_area = weighting_for_scenario(network=network_area, save=False)
+
+    return (
+        network_elec,
+        weight,
+        n_clusters,
+        busmap_foreign,
+        network_area,
+        weight_area,
+    )
 
 
 def postprocessing(
@@ -1009,6 +896,9 @@ def calc_availability_factor(gen):
     gen = network.generators[network.generators.carrier != "load shedding"][
         ["bus", "carrier", "p_nom"]
     ].copy()
+    if len(gen) == 0:
+        return pd.Series()
+
     gen["cf"] = gen.apply(calc_availability_factor, axis=1)
     gen["weight"] = gen["p_nom"] * gen["cf"]
 
@@ -1044,6 +934,49 @@ def calc_availability_factor(gen):
     return weight
 
 
+def include_busmap_area(etrago, busmap, medoid_idx, network_area, weight_area):
+    args = etrago.args["network_clustering"]
+
+    if not etrago.args["interest_area"]:
+        return busmap, medoid_idx
+    if not args["n_cluster_interest_area"]:
+        for bus in network_area.buses.index:
+            busmap[bus] = bus
+        return busmap, medoid_idx
+    else:
+        if args["method"] == "kmeans":
+            busmap_area = kmean_clustering(
+                etrago,
+                network_area,
+                weight_area,
+                args["n_cluster_interest_area"],
+            )
+            busmap_area = (
+                busmap_area.astype(int) + busmap.apply(int).max() + 1
+            ).apply(str)
+
+        if args["method"] == "kmedoids-dijkstra":
+            busmap_area, medoid_idx_area = kmedoids_dijkstra_clustering(
+                etrago,
+                network_area.buses,
+                network_area.lines,
+                weight_area,
+                args["n_cluster_interest_area"],
+            )
+
+            medoid_idx_area.index = (
+                medoid_idx_area.index.astype(int) + busmap.apply(int).max() + 1
+            )
+            busmap_area = (
+                busmap_area.astype(int) + busmap.apply(int).max() + 1
+            ).apply(str)
+            medoid_idx = pd.concat([medoid_idx, medoid_idx_area])
+
+        busmap = pd.concat([busmap, busmap_area])
+
+        return busmap, medoid_idx
+
+
 def run_spatial_clustering(self):
     """
     Main method for running spatial clustering on the electrical network.
@@ -1065,16 +998,27 @@ def run_spatial_clustering(self):
         else:
             self.disaggregated_network = self.network.copy(with_time=False)
 
-        elec_network, weight, n_clusters, busmap_foreign = preprocessing(self)
+        (
+            elec_network,
+            weight,
+            n_clusters,
+            busmap_foreign,
+            network_area,
+            weight_area,
+        ) = preprocessing(self)
 
         if self.args["network_clustering"]["method"] == "kmeans":
             if not self.args["network_clustering"]["k_elec_busmap"]:
                 logger.info("Start k-means Clustering")
 
-                busmap = kmean_clustering(
+                busmap_elec = kmean_clustering(
                     self, elec_network, weight, n_clusters
                 )
                 medoid_idx = pd.Series(dtype=str)
+                busmap, medoid_idx = include_busmap_area(
+                    self, busmap_elec, medoid_idx, network_area, weight_area
+                )
+
             else:
                 busmap = pd.Series(dtype=str)
                 medoid_idx = pd.Series(dtype=str)
@@ -1083,7 +1027,7 @@ def run_spatial_clustering(self):
             if not self.args["network_clustering"]["k_elec_busmap"]:
                 logger.info("Start k-medoids Dijkstra Clustering")
 
-                busmap, medoid_idx = kmedoids_dijkstra_clustering(
+                busmap_elec, medoid_idx = kmedoids_dijkstra_clustering(
                     self,
                     elec_network.buses,
                     elec_network.lines,
@@ -1091,6 +1035,10 @@ def run_spatial_clustering(self):
                     n_clusters,
                 )
 
+                busmap, medoid_idx = include_busmap_area(
+                    self, busmap_elec, medoid_idx, network_area, weight_area
+                )
+
             else:
                 busmap = pd.Series(dtype=str)
                 medoid_idx = pd.Series(dtype=str)
diff --git a/etrago/cluster/gas.py b/etrago/cluster/gas.py
index b31f4e1e..919ae075 100644
--- a/etrago/cluster/gas.py
+++ b/etrago/cluster/gas.py
@@ -43,7 +43,10 @@
         kmedoids_dijkstra_clustering,
         sum_with_inf,
     )
-    from etrago.tools.utilities import set_control_strategies
+    from etrago.tools.utilities import (
+        find_buses_area,
+        set_control_strategies,
+    )
 
 logger = logging.getLogger(__name__)
 
@@ -111,8 +114,8 @@ def preprocessing(etrago):
         )
     ]
 
-    # select buses dependent on whether they should be clustered in
-    # (only DE or DE+foreign)
+    # select buses dependent on whether they should be clustered in (only DE
+    # or DE+foreign DE except specified area)
     if not settings["cluster_foreign_gas"]:
         network_ch4.buses = network_ch4.buses.loc[
             ch4_filter & (network_ch4.buses["country"].values == "DE")
@@ -133,6 +136,12 @@ def preprocessing(etrago):
         network_ch4.buses = network_ch4.buses.loc[ch4_filter]
         n_clusters = settings["n_clusters_gas"]
 
+    # Exclude buses in the area that should not be clustered
+    busmap_area = find_buses_area(etrago, "CH4")
+    network_ch4.buses = network_ch4.buses[
+        ~network_ch4.buses.index.isin(busmap_area)
+    ]
+
     def weighting_for_scenario(ch4_buses, save=None):
         """
         Calculate CH4-bus weightings dependant on the connected
@@ -226,7 +235,8 @@ def weighting_for_scenario(ch4_buses, save=None):
         weight_ch4.loc[loaded_weights.index] = loaded_weights
     else:
         weight_ch4 = weighting_for_scenario(network_ch4.buses, save=False)
-    return network_ch4, weight_ch4.squeeze(axis=1), n_clusters
+
+    return network_ch4, weight_ch4.squeeze(axis=1), n_clusters, busmap_area
 
 
 def kmean_clustering_gas(etrago, network_ch4, weight, n_clusters):
@@ -288,9 +298,12 @@ def get_h2_clusters(etrago, busmap_ch4):
         to its corresponding cluster ID.
     """
     # Mapping of H2 buses to new CH4 cluster IDs
+    map_ch4_h2 = etrago.ch4_h2_mapping[
+        etrago.ch4_h2_mapping.index.isin(busmap_ch4.index)
+    ]
     busmap_h2 = pd.Series(
-        busmap_ch4.loc[etrago.ch4_h2_mapping.index].values,
-        index=etrago.ch4_h2_mapping.values,
+        busmap_ch4.loc[map_ch4_h2.index].values,
+        index=map_ch4_h2.values,
     )
 
     # Create unique H2 cluster IDs
@@ -304,7 +317,9 @@ def get_h2_clusters(etrago, busmap_ch4):
     return busmap
 
 
-def gas_postprocessing(etrago, busmap, medoid_idx=None):
+def gas_postprocessing(
+    etrago, busmap, medoid_idx=None, busmap_area=pd.Series()
+):
     """
     Performs the postprocessing for the gas grid clustering based on the
     provided busmap
@@ -359,6 +374,10 @@ def gas_postprocessing(etrago, busmap, medoid_idx=None):
     if "H2_grid" in etrago.network.buses.carrier.unique():
         busmap = get_h2_clusters(etrago, busmap)
 
+    if len(busmap_area) > 0:
+        for bus_area in busmap_area.values:
+            busmap[bus_area] = bus_area
+
     # Add all other buses to busmap
     missing_idx = list(
         etrago.network.buses[
@@ -967,7 +986,7 @@ def run_spatial_clustering_gas(self):
             method = settings["method_gas"]
             logger.info(f"Start {method} clustering GAS")
 
-            gas_network, weight, n_clusters = preprocessing(self)
+            gas_network, weight, n_clusters, busmap_area = preprocessing(self)
 
             if method == "kmeans":
                 if settings["k_gas_busmap"]:
@@ -1011,7 +1030,9 @@ def run_spatial_clustering_gas(self):
                     "spatial clustering method for the gas network"
                 )
                 raise ValueError(msg)
-            self.network, busmap = gas_postprocessing(self, busmap, medoid_idx)
+            self.network, busmap = gas_postprocessing(
+                self, busmap, medoid_idx, busmap_area
+            )
 
             self.update_busmap(busmap)
 
diff --git a/etrago/cluster/spatial.py b/etrago/cluster/spatial.py
index 9f782c1d..0dbbabe3 100755
--- a/etrago/cluster/spatial.py
+++ b/etrago/cluster/spatial.py
@@ -47,6 +47,7 @@
         buses_of_vlvl,
         connected_grid_lines,
         connected_transformer,
+        select_elec_network,
     )
 
     logger = logging.getLogger(__name__)
@@ -115,6 +116,10 @@ def strategies_buses():
 def strategies_lines():
     return {
         "geom": nan_links,
+        "cables": np.sum,
+        "topo": "first",
+        "country": "first",
+        "total_cables": np.sum,
     }
 
 
@@ -381,7 +386,7 @@ def shortest_path(paths, graph):
     return df
 
 
-def busmap_by_shortest_path(etrago, fromlvl, tolvl, cpu_cores=4):
+def busmap_by_shortest_path(network, fromlvl, tolvl, cpu_cores=4):
     """
     Creates a busmap for the EHV-Clustering between voltage levels based
     on dijkstra shortest path. The result is automatically written to the
@@ -411,15 +416,15 @@ def busmap_by_shortest_path(etrago, fromlvl, tolvl, cpu_cores=4):
     """
 
     # data preperation
-    s_buses = buses_grid_linked(etrago.network, fromlvl)
-    lines = connected_grid_lines(etrago.network, s_buses)
-    transformer = connected_transformer(etrago.network, s_buses)
-    mask = transformer.bus1.isin(buses_of_vlvl(etrago.network, tolvl))
+    s_buses = buses_grid_linked(network, fromlvl)
+    lines = connected_grid_lines(network, s_buses)
+    transformer = connected_transformer(network, s_buses)
+    mask = transformer.bus1.isin(buses_of_vlvl(network, tolvl))
 
-    dc = etrago.network.links[etrago.network.links.carrier == "DC"]
+    dc = network.links[network.links.carrier == "DC"]
     dc.index = "DC_" + dc.index
     lines_plus_dc = pd.concat([lines, dc])
-    lines_plus_dc = lines_plus_dc[etrago.network.lines.columns]
+    lines_plus_dc = lines_plus_dc[network.lines.columns]
     lines_plus_dc["carrier"] = "AC"
 
     # temporary end points, later replaced by bus1 pendant
@@ -454,18 +459,16 @@ def busmap_by_shortest_path(etrago, fromlvl, tolvl, cpu_cores=4):
     df.target = df.target.map(
         dict(
             zip(
-                etrago.network.transformers.bus0,
-                etrago.network.transformers.bus1,
+                network.transformers.bus0,
+                network.transformers.bus1,
             )
         )
     )
 
     # append to busmap buses only connected to transformer
-    transformer = etrago.network.transformers
+    transformer = network.transformers
     idx = list(
-        set(buses_of_vlvl(etrago.network, fromlvl)).symmetric_difference(
-            set(s_buses)
-        )
+        set(buses_of_vlvl(network, fromlvl)).symmetric_difference(set(s_buses))
     )
     mask = transformer.bus0.isin(idx)
 
@@ -478,7 +481,7 @@ def busmap_by_shortest_path(etrago, fromlvl, tolvl, cpu_cores=4):
     df = pd.concat([df, toappend], ignore_index=True, axis=0)
 
     # append all other buses
-    buses = etrago.network.buses[etrago.network.buses.carrier == "AC"]
+    buses = network.buses[network.buses.carrier == "AC"]
     mask = buses.index.isin(df.source)
 
     assert (buses[~mask].v_nom.astype(int).isin(tolvl)).all()
@@ -513,7 +516,6 @@ def busmap_ehv_clustering(etrago):
     busmap : dict
         Maps old bus_ids to new bus_ids.
     """
-
     if etrago.args["network_clustering_ehv"]["busmap"] is False:
         cpu_cores = etrago.args["network_clustering"]["CPU_cores"]
         if cpu_cores == "max":
@@ -521,16 +523,33 @@ def busmap_ehv_clustering(etrago):
         else:
             cpu_cores = int(cpu_cores)
 
-        busmap = busmap_by_shortest_path(
-            etrago,
-            fromlvl=[110],
-            tolvl=[220, 380, 400, 450],
-            cpu_cores=cpu_cores,
-        )
-        pd.DataFrame(busmap.items(), columns=["bus0", "bus1"]).to_csv(
-            "ehv_elecgrid_busmap_result.csv",
-            index=False,
-        )
+        if etrago.args["interest_area"] is False:
+            busmap = busmap_by_shortest_path(
+                etrago.network,
+                fromlvl=[110],
+                tolvl=[220, 380, 400, 450],
+                cpu_cores=cpu_cores,
+            )
+            pd.DataFrame(busmap.items(), columns=["bus0", "bus1"]).to_csv(
+                "ehv_elecgrid_busmap_result.csv",
+                index=False,
+            )
+        else:
+            network, _, area = select_elec_network(
+                etrago, apply_on="grid_model-ehv"
+            )
+            busmap = busmap_by_shortest_path(
+                network,
+                fromlvl=[110],
+                tolvl=[220, 380, 400, 450],
+                cpu_cores=cpu_cores,
+            )
+            for bus in area.buses.index:
+                busmap[bus] = bus
+            pd.DataFrame(busmap.items(), columns=["bus0", "bus1"]).to_csv(
+                "ehv_elecgrid_busmap_result.csv",
+                index=False,
+            )
     else:
         busmap = pd.read_csv(etrago.args["network_clustering_ehv"]["busmap"])
         busmap = pd.Series(
@@ -731,6 +750,17 @@ def kmedoids_dijkstra_clustering(
 
     settings = etrago.args["network_clustering"]
 
+    if n_clusters is False:
+        busmap = pd.Series(
+            range(len(buses)), index=buses.index, name="final_assignment"
+        )
+        busmap.index.name = "bus_id"
+        busmap = busmap.apply(str)
+
+        medoid_idx = pd.Series(busmap.index, index=busmap.values, name=0)
+
+        return busmap, medoid_idx
+
     # n_jobs was deprecated for the function fit(). scikit-learn recommends
     # to use threadpool_limits:
     # https://scikit-learn.org/stable/computing/parallelism.html
diff --git a/etrago/execute/__init__.py b/etrago/execute/__init__.py
index 7d00e861..5f70b80a 100644
--- a/etrago/execute/__init__.py
+++ b/etrago/execute/__init__.py
@@ -427,6 +427,8 @@ def optimize(self):
     None.
 
     """
+    # Verify feasibility
+    self.adjust_before_optimization()
 
     if self.args["method"]["market_optimization"]["active"]:
         self.market_optimization()
diff --git a/etrago/network.py b/etrago/network.py
index c69d0e0b..bb6d67e2 100644
--- a/etrago/network.py
+++ b/etrago/network.py
@@ -110,6 +110,7 @@
     set_random_noise,
     set_trafo_costs,
     update_busmap,
+    adjust_before_optimization,
 )
 
 logger = logging.getLogger(__name__)
@@ -236,6 +237,12 @@ def __init__(
 
             self.get_clustering_data(csv_folder_name)
 
+            conn = db.connection(section=self.args["db"])
+            self.engine = conn
+
+            session = sessionmaker(bind=conn)
+            self.session = session()
+
         else:
             logger.error("Set args or csv_folder_name")
 
@@ -372,6 +379,8 @@ def __init__(
 
     sclopf = iterate_sclopf
 
+    adjust_before_optimization = adjust_before_optimization
+
     def dc_lines(self):
         return self.filter_links_by_carrier("DC", like=False)
 
diff --git a/etrago/tools/utilities.py b/etrago/tools/utilities.py
index 3214981e..cdee1335 100755
--- a/etrago/tools/utilities.py
+++ b/etrago/tools/utilities.py
@@ -2973,3 +2973,258 @@ def manual_fixes_datamodel(etrago):
             ].index,
             inplace=True,
         )
+
+def select_elec_network(etrago, apply_on="grid_model"):
+    """
+    Creates networks to be used on the clustering based on settings specified
+    in the args.
+
+    Parameters
+    ----------
+    etrago : Etrago
+        An instance of the Etrago class
+    apply_on: str
+        gives information about the objective of the output network. If
+        "grid_model" is provided, the value assigned in the args for
+        ["network_clustering"]["cluster_foreign_AC""] will define if the
+        foreign buses will be included in the network. if "market_model" is
+        provided, foreign buses will be always included.
+
+    Returns
+    -------
+    Tuple containing:
+        elec_network : pypsa.Network
+            Contains the electric network
+        n_clusters : int
+            number of clusters used in the clustering process.
+        area_network : pypsa.Network
+            Contains the electric network in the area of interest defined in
+            network_clustering - interest_area.
+    """
+    settings = etrago.args["network_clustering"]
+
+    if apply_on in ["grid_model", "grid_model-ehv"]:
+        # Find buses in the area that should not be clustered
+        buses_area = find_buses_area(etrago, "AC")
+
+        elec_network_buses = etrago.network.buses[
+            (~etrago.network.buses.index.isin(buses_area))
+            & (etrago.network.buses.carrier == "AC")
+        ].index
+        if apply_on == "grid_model-ehv":
+            n_clusters = pd.NA
+        #Exclude foreign buses when set to don't include them in clustering
+        elif settings["cluster_foreign_AC"]:
+            n_clusters = settings["n_clusters_AC"]
+        else:
+            foreign_buses = etrago.network.buses[
+                (etrago.network.buses.country != "DE")
+                & (etrago.network.buses.carrier == "AC")
+            ]
+
+            num_neighboring_country = len(
+                foreign_buses[
+                    foreign_buses.index.isin(etrago.network.loads.bus)
+                ]
+            )
+
+            elec_network_buses = elec_network_buses[
+                ~elec_network_buses.isin(foreign_buses.index)
+            ]
+            if settings["n_clusters_AC"] is False:
+                n_clusters = False
+            else:
+                n_clusters = (
+                    settings["n_clusters_AC"] - num_neighboring_country
+                )
+                assert (
+                    n_clusters > 1
+                ), f"""'n_clusters_AC' must be greater than the number of
+                foreign countries({num_neighboring_country})"""
+
+        elec_network = network_based_on_buses(
+            etrago.network, elec_network_buses
+        )
+        area_network = network_based_on_buses(etrago.network, buses_area)
+
+    elif apply_on == "market_model":
+        elec_network_buses = etrago.network_tsa.buses[
+            etrago.network_tsa.buses.carrier == "AC"
+        ].index
+        elec_network = network_based_on_buses(
+            etrago.network_tsa, elec_network_buses
+        )
+        area_network = pypsa.Network()
+
+    else:
+        logger.warning(
+            """Parameter apply_on must be either 'grid_model' or 'market_model'
+            """
+        )
+
+    return elec_network, n_clusters, area_network
+
+
+def network_based_on_buses(network, buses):
+    """
+    Extract all the elements in a network related to the supplied list of
+    buses and return it like a new network.
+
+    Parameters
+    ----------
+    network : pypsa.Network
+        Original network that contains the buses of interest and other buses.
+    buses : Pandas.Series
+        Series that contains the name of all the buses that the new network
+        will contain.
+
+    Returns
+    -------
+    elec_network : pypsa.Network
+        network containing only electrical elements attached to the supplied
+        list of buses.
+
+    """
+    elec_network = network.copy()
+    elec_network.buses = elec_network.buses[
+        elec_network.buses.index.isin(buses)
+    ]
+    # Dealing with links
+    elec_network.links = elec_network.links[
+        (
+            (elec_network.links.carrier == "AC")
+            | (elec_network.links.carrier == "DC")
+        )
+        & (elec_network.links.bus0.isin(elec_network.buses.index))
+        & (elec_network.links.bus1.isin(elec_network.buses.index))
+    ]
+
+    # Dealing with generators
+    elec_network.generators = elec_network.generators[
+        elec_network.generators.bus.isin(elec_network.buses.index)
+    ]
+
+    for attr in elec_network.generators_t:
+        elec_network.generators_t[attr] = elec_network.generators_t[attr].loc[
+            :,
+            elec_network.generators_t[attr].columns.isin(
+                elec_network.generators.index
+            ),
+        ]
+
+    # Dealing with loads
+    elec_network.loads = elec_network.loads[
+        elec_network.loads.bus.isin(elec_network.buses.index)
+    ]
+
+    for attr in elec_network.loads_t:
+        elec_network.loads_t[attr] = elec_network.loads_t[attr].loc[
+            :,
+            elec_network.loads_t[attr].columns.isin(elec_network.loads.index),
+        ]
+
+    # Dealing with storage_units
+    elec_network.storage_units = elec_network.storage_units[
+        elec_network.storage_units.bus.isin(elec_network.buses.index)
+    ]
+
+    for attr in elec_network.storage_units_t:
+        elec_network.storage_units_t[attr] = elec_network.storage_units_t[
+            attr
+        ].loc[
+            :,
+            elec_network.storage_units_t[attr].columns.isin(
+                elec_network.storage_units.index
+            ),
+        ]
+
+    # Dealing with stores
+    elec_network.stores = elec_network.stores[
+        elec_network.stores.bus.isin(elec_network.buses.index)
+    ]
+
+    for attr in elec_network.stores_t:
+        elec_network.stores_t[attr] = elec_network.stores_t[attr].loc[
+            :,
+            elec_network.stores_t[attr].columns.isin(
+                elec_network.stores.index
+            ),
+        ]
+    return elec_network
+
+
+def find_buses_area(etrago, carrier):
+    """
+    Find buses of a specified carrier in a defined area. Usually used to
+    findout the buses that sould not be clustered.
+    """
+    settings = etrago.args
+
+    if settings["interest_area"]:
+        if isinstance(settings["interest_area"], list):
+            con = etrago.engine
+            query = "SELECT gen, geometry FROM boundaries.vg250_krs"
+
+            de_areas = gpd.read_postgis(query, con, geom_col="geometry")
+            de_areas = de_areas[
+                de_areas["gen"].isin(settings["interest_area"])
+            ]
+        elif isinstance(settings["interest_area"], str):
+            de_areas = gpd.read_file(settings["interest_area"])
+        else:
+            raise Exception(
+                "not supported format supplied to 'interest_area' argument"
+            )
+
+        try:
+            buses_area = gpd.GeoDataFrame(
+                etrago.network.buses, geometry="geom", crs=4326
+            )
+        except:
+            buses_area = etrago.network.buses[["x", "y", "carrier"]]
+            buses_area["geom"] = buses_area.apply(
+                lambda x: Point(x["x"], x["y"]), axis=1
+            )
+            buses_area = gpd.GeoDataFrame(
+                buses_area, geometry="geom", crs=4326
+            )
+
+        buses_area = gpd.clip(buses_area, de_areas)
+        buses_area = buses_area[buses_area.carrier == carrier]
+
+    else:
+        buses_area = pd.DataFrame()
+
+    return buses_area.index
+
+def adjust_before_optimization(self):
+
+    def check_e_initial(etrago):
+        stores = etrago.network.stores
+        stores_t = etrago.network.stores_t
+        for st in stores_t["e_max_pu"].columns:
+            e_initial_pu = stores.at[st, "e_initial"] / stores.at[st, "e_nom"]
+            min_e = stores_t["e_min_pu"].iloc[0, :][st]
+            max_e = stores_t["e_max_pu"].iloc[0, :][st]
+            if (e_initial_pu >= min_e) & (e_initial_pu <= max_e):
+                continue
+            else:
+                stores.at[st, "e_initial"] = (
+                    stores.at[st, "e_nom"] * (min_e + max_e) / 2
+                )
+
+        return stores
+
+    # Temporary drop DLR as it is currently not working with sclopf
+    if self.args["method"]["type"] != "lopf":
+        self.network.lines_t.s_max_pu = pd.DataFrame(
+            index=self.network.snapshots,
+            columns=self.network.lines.index,
+            data=1.0,
+        )
+
+    self.network.storage_units.cyclic_state_of_charge = True
+
+    self.network.lines.loc[self.network.lines.r == 0.0, "r"] = 10
+
+    self.network.stores = check_e_initial(self)
\ No newline at end of file