Add the option to define, create and remove subregions

Snakefile

@@ -33,6 +33,14 @@ config.update({"git_commit": get_last_commit_message(".")})
 # convert country list according to the desired region
 config["countries"] = create_country_list(config["countries"])
 
+# extend countries based on subregions and remove the dropped subregions
+config["countries_plus"] = [
+    x
+    for x in config["countries"]
+    + list(config.get("subregion", {}).get("define_subregion_gadm", {}).keys())
+    if (x not in config.get("subregion", {}).get("drop_subregion", {}))
+]
+
 # create a list of iteration steps, required to solve the experimental design
 # each value is used as wildcard input e.g. solution_{unc}
 config["scenario"]["unc"] = [
@@ -128,7 +136,7 @@ rule make_all_summaries:
             + RDIR
             + "summaries/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}_{country}",
             **config["scenario"],
-            country=["all"] + config["countries"],
+            country=["all"] + config["countries_plus"],
         ),
 
 
@@ -140,7 +148,7 @@ rule plot_all_summaries:
             + "plots/summary_{summary}_elec_s{simpl}_{clusters}_ec_l{ll}_{opts}_{country}.{ext}",
             summary=["energy", "costs"],
             **config["scenario"],
-            country=["all"] + config["countries"],
+            country=["all"] + config["countries_plus"],
             ext=["png", "pdf"],
         ),
 
@@ -210,7 +218,7 @@ rule clean_osm_data:
 rule build_osm_network:
     params:
         build_osm_network=config.get("build_osm_network", {}),
-        countries=config["countries"],
+        countries=config["countries_plus"],
         crs=config["crs"],
     input:
         generators="resources/" + RDIR + "osm/clean/all_clean_generators.geojson",
@@ -237,6 +245,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',
@@ -268,7 +277,7 @@ rule base_network:
         links=config["links"],
         lines=config["lines"],
         hvdc_as_lines=config["electricity"]["hvdc_as_lines"],
-        countries=config["countries"],
+        countries=config["countries_plus"],
         base_network=config["base_network"],
     input:
         osm_buses="resources/" + RDIR + "base_network/all_buses_build_network.csv",
@@ -298,7 +307,7 @@ rule build_bus_regions:
     params:
         alternative_clustering=config["cluster_options"]["alternative_clustering"],
         crs=config["crs"],
-        countries=config["countries"],
+        countries=config["countries_plus"],
     input:
         country_shapes="resources/" + RDIR + "shapes/country_shapes.geojson",
         offshore_shapes="resources/" + RDIR + "shapes/offshore_shapes.geojson",
@@ -444,7 +453,7 @@ rule build_renewable_profiles:
     params:
         crs=config["crs"],
         renewable=config["renewable"],
-        countries=config["countries"],
+        countries=config["countries_plus"],
         alternative_clustering=config["cluster_options"]["alternative_clustering"],
     input:
         natura="resources/" + RDIR + "natura.tiff",
@@ -481,6 +490,8 @@ rule build_powerplants:
     params:
         geo_crs=config["crs"]["geo_crs"],
         countries=config["countries"],
+        countries_plus=config["countries_plus"],
+        subregion=config.get("subregion", {}),
         gadm_layer_id=config["build_shape_options"]["gadm_layer_id"],
         alternative_clustering=config["cluster_options"]["alternative_clustering"],
         powerplants_filter=config["electricity"]["powerplants_filter"],
@@ -558,8 +569,9 @@ rule simplify_network:
     params:
         renewable=config["renewable"],
         geo_crs=config["crs"]["geo_crs"],
+        distance_crs=config["crs"]["distance_crs"],
         cluster_options=config["cluster_options"],
-        countries=config["countries"],
+        countries=config["countries_plus"],
         build_shape_options=config["build_shape_options"],
         electricity=config["electricity"],
         costs=config["costs"],
@@ -604,7 +616,7 @@ if config["augmented_line_connection"].get("add_to_snakefile", False) == True:
             length_factor=config["lines"]["length_factor"],
             renewable=config["renewable"],
             geo_crs=config["crs"]["geo_crs"],
-            countries=config["countries"],
+            countries=config["countries_plus"],
             cluster_options=config["cluster_options"],
             focus_weights=config.get("focus_weights", None),
             #custom_busmap=config["enable"].get("custom_busmap", False)
@@ -689,7 +701,7 @@ if config["augmented_line_connection"].get("add_to_snakefile", False) == False:
             length_factor=config["lines"]["length_factor"],
             renewable=config["renewable"],
             geo_crs=config["crs"]["geo_crs"],
-            countries=config["countries"],
+            countries=config["countries_plus"],
             gadm_layer_id=config["build_shape_options"]["gadm_layer_id"],
             cluster_options=config["cluster_options"],
         input:
@@ -1023,7 +1035,7 @@ rule build_test_configs:
 
 rule make_statistics:
     params:
-        countries=config["countries"],
+        countries=config["countries_plus"],
         renewable_carriers=config["electricity"]["renewable_carriers"],
         renewable=config["renewable"],
         crs=config["crs"],

config.default.yaml

@@ -12,6 +12,17 @@ logging:
 countries: ["NG", "BJ"]
 # Can be replaced by country ["NG", "BJ"], continent ["Africa"] or user-specific region, see more at https://pypsa-earth.readthedocs.io/en/latest/configuration.html#top-level-configuration
 
+subregion: []
+# define_subregion_gadm:
+#   Lagos-Porto-Novo:
+#     NG: [25,28]
+#     BJ: [10,11]
+# Chose the subregion based on its GADM subdivision Map ID code.
+# Chose the country in https://geodata.ucdavis.edu/gadm/gadm4.1/gpkg and use https://ngageoint.github.io/geopackage-viewer-js/ and chose ADM_ADM_1 to identify the GADM ID
+# This will take away the land of the GADM ID code from the original country to the new country.
+# One subregion can be composed of multiple countries.
+# drop_subregion: [] # remove the country/subregion from the model.
+
 enable:
   retrieve_databundle: true  #  Recommended 'true', for the first run. Otherwise data might be missing.
   retrieve_cost_data: true  # true: retrieves cost data from technology data and saves in resources/costs.csv, false: uses cost data in data/costs.csv

scripts/add_electricity.py

@@ -740,6 +740,24 @@ def estimate_renewable_capacities_irena(
     capacities = capacities.query(
         "Year == @year and Technology in @tech_keys and Country in @countries"
     )
+    custom_capacity = estimate_renewable_capacities_config.get("custom_capacity", {})
+    if custom_capacity:
+        custom_capacities_df = (
+            pd.DataFrame.from_dict(custom_capacity)
+            .reset_index()
+            .melt(id_vars=["index"], value_vars=["Onshore", "PV"])
+            .dropna()
+        )
+        custom_capacities_df = custom_capacities_df.rename(
+            columns={"index": "Country", "variable": "Technology", "value": "Capacity"}
+        )
+
+        capacities = pd.concat([capacities, custom_capacities_df])
+
+    if all(country in n.buses.country.unique().tolist() for country in countries):
+        countries = n.buses.country.unique().tolist()
+        capacities = capacities.query("Country == @countries")
+
     capacities = capacities.groupby(["Technology", "Country"]).Capacity.sum()
 
     logger.info(

scripts/build_demand_profiles.py

@@ -178,6 +178,7 @@ def build_demand_profiles(
     admin_shapes,
     countries,
     scale,
+    substitute_country_load,
     start_date,
     end_date,
     out_path,
@@ -228,6 +229,21 @@ def build_demand_profiles(
     # filter load for analysed countries
     gegis_load = gegis_load.loc[gegis_load.region_code.isin(countries)]
 
+    # If Gegis data is unavailable for a specific country, it is advisable to use Gegis data from another region as an alternative. Later, it can be scaled
+    if substitute_country_load:
+        for country_a, country_b in substitute_country_load.items():
+            logger.info(f"substitute load data of {country_a} using {country_b}.")
+
+            if country_a in gegis_load.region_code.unique():
+                logger.info(f"dropping original load data of {country_a}.")
+                gegis_load = gegis_load.query("region_code != @country_a")
+
+            gegis_load_new = gegis_load.loc[gegis_load.region_code == country_b]
+            gegis_load_new.loc[:, "region_code"] = country_a
+            gegis_load_new.loc[:, "region_name"] = country_a
+
+            gegis_load = pd.concat([gegis_load, gegis_load_new])
+
     if isinstance(scale, dict):
         logger.info(f"Using custom scaling factor for load data.")
         DEFAULT_VAL = scale.get("DEFAULT", 1.0)
@@ -307,6 +323,9 @@ def upsample(cntry, group):
     countries = snakemake.params.countries
     admin_shapes = snakemake.input.gadm_shapes
     scale = snakemake.params.load_options.get("scale", 1.0)
+    substitute_country_load = snakemake.params.load_options.get(
+        "substitute_country_load", False
+    )
     start_date = snakemake.params.snapshots["start"]
     end_date = snakemake.params.snapshots["end"]
     out_path = snakemake.output[0]
@@ -318,6 +337,7 @@ def upsample(cntry, group):
         admin_shapes,
         countries,
         scale,
+        substitute_country_load,
         start_date,
         end_date,
         out_path,

scripts/build_powerplants.py

@@ -310,7 +310,11 @@ def replace_natural_gas_technology(df: pd.DataFrame):
     filepath_osm2pm_ppl = snakemake.output.powerplants_osm2pm
 
     n = pypsa.Network(snakemake.input.base_network)
-    countries_codes = n.buses.country.unique()
+
+    if snakemake.params.subregion:
+        countries_codes = snakemake.params.countries
+    else:
+        countries_codes = n.buses.country.unique()
     countries_names = list(map(two_digits_2_name_country, countries_codes))
 
     config["target_countries"] = countries_names
@@ -362,6 +366,19 @@ def replace_natural_gas_technology(df: pd.DataFrame):
         tree_i = kdtree.query(ppl.loc[ppl_i, ["lon", "lat"]].values)[1]
         ppl.loc[ppl_i, "bus"] = substation_i.append(pd.Index([np.nan]))[tree_i]
 
+    if snakemake.params.subregion:
+        subregion_include = list(snakemake.params.subregion["define_subregion_gadm"])
+        subregion_exclude = snakemake.params.subregion["drop_subregion"]
+
+        substation_i = n.buses.query(
+            "substation_lv and country.isin(@subregion_include)"
+        ).index
+        kdtree = KDTree(n.buses.loc[substation_i, ["x", "y"]].values)
+        ppl_i = ppl.query("Country.isin(@subregion_exclude)").index
+
+        tree_i = kdtree.query(ppl.loc[ppl_i, ["lon", "lat"]].values)[1]
+        ppl.loc[ppl_i, "bus"] = substation_i.append(pd.Index([np.nan]))[tree_i]
+
     if cntries_without_ppl:
         logger.warning(f"No powerplants known in: {', '.join(cntries_without_ppl)}")
 

scripts/build_shapes.py

@@ -1306,6 +1306,84 @@ def gadm(
     return df_gadm
 
 
+def crop_country(country_shapes, gadm_shapes, subregion_config):
+    country_shapes_new = gpd.GeoDataFrame(columns=["name", "geometry"]).set_index(
+        "name"
+    )
+
+    remain_gadm_shapes = gadm_shapes.copy(deep=True)
+
+    for sub_region in subregion_config:
+        region_GADM = [
+            country + "." + str(region) + "_1"
+            for country in subregion_config[sub_region]
+            for region in subregion_config[sub_region][country]
+        ]
+
+        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,
+    )
+
+
+def crop_offshore(offshore_shapes, gadm_shapes, subregion_config):
+    import cartopy.crs as ccrs
+    from build_bus_regions import custom_voronoi_partition_pts
+
+    offshore_regions_voronoi = gpd.GeoDataFrame(columns=["name", "geometry"]).set_index(
+        "name"
+    )
+
+    for country in offshore_shapes.index:
+        gadm_shapes_country = gadm_shapes.loc[gadm_shapes.country == country]
+        gadm_points = gadm_shapes_country.to_crs(ccrs.PlateCarree().proj4_init).centroid
+        gadm_points_gdf = pd.DataFrame(data={"x": gadm_points.x, "y": gadm_points.y})
+        offshore_geometry = custom_voronoi_partition_pts(
+            gadm_points_gdf, offshore_shapes.geometry.loc[country]
+        )
+
+        offshore_regions = gpd.GeoDataFrame(
+            {
+                "name": gadm_shapes_country.index,
+                "geometry": offshore_geometry,
+            },
+        ).set_index("name")
+
+        offshore_regions_voronoi = pd.concat(
+            [offshore_regions_voronoi, offshore_regions]
+        )
+
+    offshore_regions_voronoi["country"] = offshore_regions_voronoi.index.str[:2]
+
+    return crop_country(offshore_shapes, offshore_regions_voronoi, subregion_config)
+
+
 if __name__ == "__main__":
     if "snakemake" not in globals():
         from _helpers import mock_snakemake
@@ -1340,19 +1418,11 @@ def gadm(
         update,
         out_logging,
     )
-    country_shapes.to_file(snakemake.output.country_shapes)
 
     offshore_shapes = eez(
         countries_list, geo_crs, country_shapes, EEZ_gpkg, out_logging
     )
 
-    offshore_shapes.reset_index().to_file(snakemake.output.offshore_shapes)
-
-    africa_shape = gpd.GeoDataFrame(
-        geometry=[country_cover(country_shapes, offshore_shapes.geometry)]
-    )
-    africa_shape.reset_index().to_file(snakemake.output.africa_shape)
-
     gadm_shapes = gadm(
         worldpop_method,
         gdp_method,
@@ -1366,4 +1436,35 @@ def gadm(
         year,
         nprocesses=nprocesses,
     )
+
+    if snakemake.params.get("subregion", False):
+        subregion_config = snakemake.params.subregion["define_subregion_gadm"]
+        country_shapes = crop_country(country_shapes, gadm_shapes, subregion_config)
+        offshore_shapes = crop_offshore(offshore_shapes, gadm_shapes, subregion_config)
+
+        for sub_region in subregion_config:
+            region_GADM = [
+                country + "." + str(region) + "_1"
+                for country in subregion_config[sub_region]
+                for region in subregion_config[sub_region][country]
+            ]
+
+            gadm_shapes.loc[region_GADM, "country"] = sub_region
+
+    if snakemake.params.get("subregion", {}).get("drop_subregion", False):
+        drop_subregion = snakemake.params.subregion["drop_subregion"]
+        country_shapes = country_shapes.query("index != @drop_subregion")
+        offshore_shapes = offshore_shapes.query("index != @drop_subregion")
+        gadm_shapes = gadm_shapes.query("country != @drop_subregion")
+
+    africa_shape = gpd.GeoDataFrame(
+        geometry=[country_cover(country_shapes.geometry, offshore_shapes.geometry)]
+    )
+
+    country_shapes.to_file(snakemake.output.country_shapes)
+
+    offshore_shapes.reset_index().to_file(snakemake.output.offshore_shapes)
+
     save_to_geojson(gadm_shapes, out.gadm_shapes)
+
+    africa_shape.reset_index().to_file(snakemake.output.africa_shape)