Fix workflow for running Africa

Snakefile

@@ -1,8 +1,13 @@
+import sys
+
+sys.path.append("pypsa-earth/scripts")
+
 from os.path import exists
 from shutil import copyfile, move
 from scripts.helpers import get_last_commit_message
 
 from snakemake.remote.HTTP import RemoteProvider as HTTPRemoteProvider
+from _helpers import create_country_list
 
 HTTP = HTTPRemoteProvider()
 
@@ -13,7 +18,10 @@ if not exists("config.yaml"):
 configfile: "config.yaml"
 
 
-PYPSAEARTH_FOLDER = "./pypsa-earth"
+PYPSAEARTH_FOLDER = "pypsa-earth"
+
+# convert country list according to the desired region
+config["countries"] = create_country_list(config["countries"])
 
 
 SDIR = config["summary_dir"] + config["run"]
@@ -41,13 +49,21 @@ wildcard_constraints:
     h2export="[0-9]+m?|all",
 
 
-subworkflow pypsaearth:
-    workdir:
-        PYPSAEARTH_FOLDER
-    snakefile:
-        PYPSAEARTH_FOLDER + "/Snakefile"
-    configfile:
-        "./config.pypsa-earth.yaml"
+if not config.get("disable_subworkflow", False):
+
+    subworkflow pypsaearth:
+        workdir:
+            PYPSAEARTH_FOLDER
+        snakefile:
+            PYPSAEARTH_FOLDER + "/Snakefile"
+        configfile:
+            "./config.pypsa-earth.yaml"
+
+
+if config.get("disable_subworkflow", False):
+
+    def pypsaearth(path):
+        return PYPSAEARTH_FOLDER + "/" + path
 
 
 if config["enable"].get("retrieve_cost_data", True):

config.default.yaml

@@ -9,6 +9,9 @@ run: test_run
 
 foresight: overnight
 
+# option to disable the subworkflow to ease the analyses
+disable_subworkflow: false
+
 scenario:
   simpl: # only relevant for PyPSA-Eur
   - ""

config.pypsa-earth.yaml

@@ -2,7 +2,7 @@
 #
 # SPDX-License-Identifier: CC0-1.0
 
-version: 0.2.3
+version: 0.3.0
 tutorial: false
 
 logging:
@@ -16,7 +16,7 @@ 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
   download_osm_data: true  # If 'true', OpenStreetMap data will be downloaded for the above given countries
-  build_natura_raster: false  # If True, than an exclusion raster will be build
+  build_natura_raster: false  # If True, then an exclusion raster will be build
   build_cutout: false
   # If "build_cutout" : true, then environmental data is extracted according to `snapshots` date range and `countries`
   # requires cds API key https://cds.climate.copernicus.eu/api-how-to
@@ -62,9 +62,10 @@ cluster_options:
     feature: solar+onwind-time # only for hac. choose from: [solar+onwind-time, solar+onwind-cap, solar-time, solar-cap, solar+offwind-cap] etc.
     exclude_carriers: []
     remove_stubs: true
-    remove_stubs_across_borders: true
+    remove_stubs_across_borders: false
     p_threshold_drop_isolated: 20 # [MW] isolated buses are being discarded if bus mean power is below the specified threshold
-    p_threshold_merge_isolated: 300 # [MW] isolated buses are being merged into a single isolated bus if bus mean power is below the specified threshold
+    p_threshold_merge_isolated: 300 # [MW] isolated buses are being merged into a single isolated bus if a bus mean power is below the specified threshold
+    s_threshold_fetch_isolated: 0.05 # [-] a share of the national load for merging an isolated network into a backbone network
   cluster_network:
     algorithm: kmeans
     feature: solar+onwind-time
@@ -109,7 +110,7 @@ clean_osm_data_options:  # osm = OpenStreetMap
 
 build_osm_network:  # Options of the build_osm_network script; osm = OpenStreetMap
   group_close_buses: true  # When "True", close buses are merged and guarantee the voltage matching among line endings
-  group_tolerance_buses: 5000  # [m] (default 5000) Tolerance in meters of the close buses to merge
+  group_tolerance_buses: 10000  # [m] (default 5000) Tolerance in meters of the close buses to merge
   split_overpassing_lines: true  # When True, lines overpassing buses are splitted and connected to the bueses
   overpassing_lines_tolerance: 1  # [m] (default 1) Tolerance to identify lines overpassing buses
   force_ac: false  # When true, it forces all components (lines and substation) to be AC-only. To be used if DC assets create problem.
@@ -126,9 +127,9 @@ load_options:
 
 electricity:
   base_voltage: 380.
-  voltages: [220., 300., 380.]
-  co2limit: 72.0e+6  # European default, 0.05 * 3.1e9*0.5, needs to be adjusted for Africa
-  co2base: 72.0e+6  # European default, adjustment to Africa necessary
+  voltages: [132., 220., 300., 380., 500., 750.]
+  co2limit: 7.75e+7  # European default, 0.05 * 3.1e9*0.5, needs to be adjusted for Africa
+  co2base: 1.487e+9  # European default, adjustment to Africa necessary
   agg_p_nom_limits: data/agg_p_nom_minmax.csv
   hvdc_as_lines: false  # should HVDC lines be modeled as `Line` or as `Link` component?
   automatic_emission: false
@@ -168,11 +169,15 @@ electricity:
       PV: [solar]
 
 lines:
-  types:
+  ac_types:
+    132.: "243-AL1/39-ST1A 20.0"
     220.: "Al/St 240/40 2-bundle 220.0"
     300.: "Al/St 240/40 3-bundle 300.0"
     380.: "Al/St 240/40 4-bundle 380.0"
-  dc_type: "HVDC XLPE 1000"
+    500.: "Al/St 240/40 4-bundle 380.0"
+    750.: "Al/St 560/50 4-bundle 750.0"
+  dc_types:
+    500.: "HVDC XLPE 1000"
   s_max_pu: 0.7
   s_nom_max: .inf
   length_factor: 1.25

scripts/build_base_industry_totals.py

@@ -26,7 +26,7 @@ def create_industry_base_totals(df):
     # Converting values of mass (ktons) to energy (TWh)
     index_mass = df.loc[df["Unit"] == "Metric tons,  thousand"].index
     df.loc[index_mass, "Quantity_TWh"] = df.loc[index_mass].apply(
-        lambda x: x["Quantity"] * fuels_conv_toTWh.get(x["Commodity"], "not found"),
+        lambda x: x["Quantity"] * fuels_conv_toTWh.get(x["Commodity"], float("nan")),
         axis=1,
     )
 

scripts/build_clustered_population_layouts.py

@@ -6,6 +6,7 @@
 import geopandas as gpd
 import pandas as pd
 import xarray as xr
+from helpers import read_csv_nafix, to_csv_nafix
 
 if __name__ == "__main__":
     if "snakemake" not in globals():
@@ -43,15 +44,15 @@
 
     pop["ct"] = gpd.read_file(snakemake.input.regions_onshore).set_index("name").country
     country_population = pop.total.groupby(pop.ct).sum()
-    pop["fraction"] = pop.total / pop.ct.map(country_population)
+    pop["fraction"] = (pop.total / pop.ct.map(country_population)).fillna(0.0)
 
-    pop.to_csv(snakemake.output.clustered_pop_layout)
+    to_csv_nafix(pop, snakemake.output.clustered_pop_layout)
 
     gdp_layout = xr.open_dataarray(snakemake.input["gdp_layout"])
     gdp = I.dot(gdp_layout.stack(spatial=("y", "x")))
     gdp = pd.DataFrame(gdp, index=clustered_regions.index, columns=["total"])
 
     gdp["ct"] = gpd.read_file(snakemake.input.regions_onshore).set_index("name").country
     country_gdp = gdp.total.groupby(gdp.ct).sum()
-    gdp["fraction"] = gdp.total / gdp.ct.map(country_gdp)
-    gdp.to_csv(snakemake.output.clustered_gdp_layout)
+    gdp["fraction"] = (gdp.total / gdp.ct.map(country_gdp)).fillna(0.0)
+    to_csv_nafix(gdp, snakemake.output.clustered_gdp_layout)

scripts/build_industry_demand.py

@@ -108,7 +108,6 @@ def country_to_nodal(industrial_production, keys):
         cagr = read_csv_nafix(snakemake.input.industry_growth_cagr, index_col=0)
 
         # Building nodal industry production growth
-
         for country in countries:
             if country not in cagr.index:
                 cagr.loc[country] = cagr.loc["DEFAULT"]
@@ -129,7 +128,8 @@ def country_to_nodal(industrial_production, keys):
         production_base = cagr.applymap(lambda x: 1)
         production_tom = production_base * growth_factors
 
-        industry_totals = (production_tom * industry_base_totals).fillna(0)
+        # non-used line; commented out
+        # industry_totals = (production_tom * industry_base_totals).fillna(0)
 
         industry_util_factor = snakemake.config["sector"]["industry_util_factor"]
 
@@ -175,6 +175,11 @@ def country_to_nodal(industrial_production, keys):
             "other": 0,
         }
 
+        # fill industry_base_totals
+        level_2nd = industry_base_totals.index.get_level_values(1).unique()
+        mlv_index = pd.MultiIndex.from_product([countries, level_2nd])
+        industry_base_totals = industry_base_totals.reindex(mlv_index, fill_value=0)
+
         geo_locs = pd.read_csv(
             snakemake.input.industrial_database,
             sep=",",
@@ -199,12 +204,15 @@ def match_technology(df):
 
         # Calculating emissions
 
-        geo_locs = match_technology(geo_locs).loc[countries]
+        # get the subset of countries that al
+        countries_geo = geo_locs.index.unique().intersection(countries)
+        geo_locs = match_technology(geo_locs).loc[countries_geo]
 
+        aluminium_year = snakemake.config["demand_data"]["aluminium_year"]
         AL = read_csv_nafix("data/AL_production.csv", index_col=0)
-        AL_prod_tom = AL[AL.Year == snakemake.config["demand_data"]["aluminium_year"]][
+        AL_prod_tom = AL.query("Year == @aluminium_year and index in @countries_geo")[
             "production[ktons/a]"
-        ].loc[countries]
+        ].reindex(countries_geo, fill_value=0.0)
         AL_emissions = AL_prod_tom * emission_factors["non-ferrous metals"]
 
         Steel_emissions = (
@@ -277,12 +285,20 @@ def match_technology(df):
 
         # Fill missing carriers with 0s
         for country in countries:
-            carriers_present = industry_base_totals.xs(country, level="country").index
+            carriers_present = industry_base_totals.xs(country, level=0).index
             missing_carriers = set(all_carriers) - set(carriers_present)
             for carrier in missing_carriers:
                 # Add the missing carrier with a value of 0
                 industry_base_totals.loc[(country, carrier), :] = 0
 
+        # temporary fix: merge other manufacturing, construction and non-fuel into other and drop the column
+        other_cols = list(set(industry_base_totals.columns) - set(clean_industry_list))
+        if len(other_cols) > 0:
+            industry_base_totals["other"] += industry_base_totals[other_cols].sum(
+                axis=1
+            )
+            industry_base_totals.drop(columns=other_cols, inplace=True)
+
         nodal_df = pd.DataFrame()
 
         for country in countries:

scripts/build_population_layouts.py

@@ -65,6 +65,7 @@
 
     # Squeeze into a Series
     urban_fraction = urban_percent_df.squeeze() / 100.0
+    urban_fraction = urban_fraction.groupby(urban_fraction.index).sum()
 
     # population in each grid cell
     pop_cells = pd.Series(I.dot(nuts3["pop"]))

scripts/helpers.py

@@ -568,7 +568,9 @@ def locate_bus(
     except ValueError:
         return gdf_co[gdf_co.geometry == min(gdf_co.geometry, key=(point.distance))][
             col
-        ].item()  # looks for closest one shape=node
+        ].iloc[
+            0
+        ]  # looks for closest one shape=node
 
 
 def get_conv_factors(sector):
@@ -748,6 +750,18 @@ def read_csv_nafix(file, **kwargs):
         return pd.DataFrame()
 
 
+def to_csv_nafix(df, path, **kwargs):
+    "Function to export a pandas object into a csv and standardize the na value"
+    if "na_rep" in kwargs:
+        del kwargs["na_rep"]
+    # if len(df) > 0:
+    if not df.empty or not df.columns.empty:
+        return df.to_csv(path, **kwargs, na_rep=NA_VALUES[0])
+    else:
+        with open(path, "w") as fp:
+            pass
+
+
 def safe_divide(numerator, denominator, default_value=np.nan):
     """
     Safe division function that returns NaN when the denominator is zero

scripts/prepare_gas_network.py

@@ -460,7 +460,10 @@ def get_GADM_layer(
 
         # create a subindex column that is useful
         # in the GADM processing of sub-national zones
-        geodf_temp["GADM_ID"] = geodf_temp[f"GID_{cur_layer_id}"]
+        # Fix issues with missing "." in selected cases
+        geodf_temp["GADM_ID"] = geodf_temp[f"GID_{cur_layer_id}"].apply(
+            lambda x: x if x[3] == "." else x[:3] + "." + x[3:]
+        )
 
         # append geodataframes
         geodf_list.append(geodf_temp)