add option to define country as subregion and restore them back

Snakefile

@@ -227,6 +227,7 @@ rule build_shapes:
         build_shape_options=config["build_shape_options"],
         crs=config["crs"],
         countries=config["countries"],
+        subregion=config.get("subregion",{}),
     input:
         # naturalearth='data/bundle/naturalearth/ne_10m_admin_0_countries.shp',
         # eez='data/bundle/eez/World_EEZ_v8_2014.shp',
@@ -239,6 +240,8 @@ rule build_shapes:
         offshore_shapes="resources/" + RDIR + "shapes/offshore_shapes.geojson",
         africa_shape="resources/" + RDIR + "shapes/africa_shape.geojson",
         gadm_shapes="resources/" + RDIR + "shapes/gadm_shapes.geojson",
+        subregion_shapes=("resources/" + RDIR + "shapes/subregion_shapes.geojson"
+                          if config.get("subregion", {}).get("define_by_gadm", False) else []),
     log:
         "logs/" + RDIR + "build_shapes.log",
     benchmark:
@@ -553,7 +556,7 @@ rule add_electricity:
 rule simplify_network:
     params:
         renewable=config["renewable"],
-        geo_crs=config["crs"]["geo_crs"],
+        crs=config["crs"],
         cluster_options=config["cluster_options"],
         countries=config["countries"],
         build_shape_options=config["build_shape_options"],
@@ -567,6 +570,11 @@ rule simplify_network:
         tech_costs=COSTS,
         regions_onshore="resources/" + RDIR + "bus_regions/regions_onshore.geojson",
         regions_offshore="resources/" + RDIR + "bus_regions/regions_offshore.geojson",
+        country_shapes=("resources/" + RDIR + "shapes/country_shapes.geojson"
+                          if config.get("subregion") else []),
+        subregion_shapes=("resources/" + RDIR + "shapes/subregion_shapes.geojson"
+                          or config.get("subregion",{}).get("path_custom_shapes") 
+                          if config.get("subregion") else []),
     output:
         network="networks/" + RDIR + "elec_s{simpl}.nc",
         regions_onshore="resources/"

scripts/build_shapes.py

@@ -1318,6 +1318,56 @@ def gadm(
 
     return df_gadm
 
+def crop_country(gadm_shapes, subregion_config):
+    country_shapes_new = gpd.GeoDataFrame(columns=["name", "geometry"]).set_index(
+        "name"
+    )
+
+    remain_gadm_shapes = gadm_shapes.copy(deep=True)
+    country_gadm = remain_gadm_shapes["country"].unique()
+
+    for sub_region in subregion_config:
+
+        region_GADM = [
+            country + "." + str(region) + "_1"
+            for country in subregion_config[sub_region]
+            if country in country_gadm
+            for region in subregion_config[sub_region][country]
+        ]
+
+        if not region_GADM:
+            continue
+
+        sub_country_geometry = gadm_shapes.loc[region_GADM].unary_union
+        sub_country_shapes = gpd.GeoDataFrame(
+            {
+                "name": sub_region,
+                "geometry": [sub_country_geometry],
+            }
+        ).set_index("name")
+
+        country_shapes_new = pd.concat([country_shapes_new, sub_country_shapes])
+
+        remain_gadm_shapes = remain_gadm_shapes.query("index != @region_GADM")
+
+    for country in remain_gadm_shapes.country.unique():
+        country_geometry_new = remain_gadm_shapes.query(
+            "country == @country"
+        ).unary_union
+        country_shapes_country = gpd.GeoDataFrame(
+            {
+                "name": country,
+                "geometry": [country_geometry_new],
+            }
+        ).set_index("name")
+
+        country_shapes_new = pd.concat([country_shapes_new, country_shapes_country])
+
+    return gpd.GeoDataFrame(
+        country_shapes_new,
+        crs=offshore_shapes.crs,
+        geometry=country_shapes_new.geometry,
+    )
 
 if __name__ == "__main__":
     if "snakemake" not in globals():
@@ -1379,4 +1429,10 @@ def gadm(
         nprocesses=nprocesses,
         simplify_gadm=simplify_gadm,
     )
+
+    if snakemake.params.get("subregion").get("define_by_gadm", False):
+        subregion_config = snakemake.params.subregion["define_by_gadm"]
+        subregion_shapes = crop_country(gadm_shapes, subregion_config)
+        subregion_shapes.to_file(snakemake.output.subregion_shapes)
+
     save_to_geojson(gadm_shapes, out.gadm_shapes)

scripts/simplify_network.py

@@ -799,7 +799,7 @@ def merge_into_network(n, threshold, aggregation_strategies=dict()):
     generators_mean_origin = n.generators.p_nom.mean()
     load_mean_origin = n.loads_t.p_set.mean().mean()
 
-    network_crs = snakemake.params.geo_crs
+    network_crs = snakemake.params.crs["geo_crs"]
 
     n.determine_network_topology()
 
@@ -959,6 +959,16 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()):
 
     return clustering.network, busmap
 
+def nearest_shape(n, path_shapes, distance_crs):
+    from shapely.geometry import Point
+
+    shapes = gpd.read_file(path_shapes, crs=distance_crs).set_index("name")["geometry"]
+
+    for i in n.buses.index:
+        point = Point(n.buses.loc[i,"x"], n.buses.loc[i,"y"])
+        n.buses.loc[i,"country"] = shapes.distance(point).sort_values().index[0]
+
+    return n
 
 if __name__ == "__main__":
     if "snakemake" not in globals():
@@ -1070,7 +1080,7 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()):
         distribution_cluster = snakemake.params.cluster_options["distribute_cluster"]
         focus_weights = snakemake.params.focus_weights
         gadm_layer_id = snakemake.params.build_shape_options["gadm_layer_id"]
-        geo_crs = snakemake.params.geo_crs
+        geo_crs = snakemake.params.crs["geo_crs"]
         renewable_config = snakemake.params.renewable
         solver_name = snakemake.config["solving"]["solver"]["name"]
 
@@ -1106,6 +1116,11 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()):
 
     update_p_nom_max(n)
 
+    subregion_shapes = snakemake.input.subregion_shapes
+    if subregion_shapes:
+        distance_crs = snakemake.params.crs["distance_crs"]
+        n = nearest_shape(n, subregion_shapes, distance_crs)
+
     p_threshold_drop_isolated = max(
         0.0, cluster_config.get("p_threshold_drop_isolated", 0.0)
     )
@@ -1130,6 +1145,10 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()):
         )
         busmaps.append(fetched_nodes_map)
 
+    if subregion_shapes:
+        country_shapes = snakemake.input.country_shapes
+        n = nearest_shape(n, country_shapes, distance_crs)
+
     n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards)))
     n.export_to_netcdf(snakemake.output.network)