Various bugfixes to make multiyear branch run

rules/build_electricity.smk

@@ -202,12 +202,26 @@ if config["enable"].get("build_cutout", False):
             "../scripts/build_cutout.py"
 
 
+# TODO: consider if it's fine to just use the default cutout instead, i.e.:
+# cutout="cutouts/" + CDIR + config["atlite"]["default_cutout"] + ".nc",
+def cutouts_for_extent_input(wildcards):
+    # We only need these cutouts in order to determine geographic
+    # boundaries, so just pick the first weather year available
+    weather_year = config["scenario"]["weather_year"][0]
+    cutouts = []
+    for k in config["electricity"]["renewable_carriers"]:
+        cutout = config["renewable"][k]["cutout"]
+        cutout = cutout.replace("{weather_year}", weather_year)
+        cutouts.append(f"cutouts/{CDIR}{cutout}.nc")
+    return cutouts
+
+
 if config["enable"].get("build_natura_raster", False):
 
     rule build_natura_raster:
         input:
             natura=ancient("data/bundle/natura/Natura2000_end2015.shp"),
-            cutouts=expand("cutouts/" + CDIR + "{cutouts}.nc", **config["atlite"]),
+            cutouts=cutouts_for_extent_input,
         output:
             RESOURCES + "natura.tiff",
         resources:
@@ -223,13 +237,7 @@ if config["enable"].get("build_natura_raster", False):
 rule build_ship_raster:
     input:
         ship_density="data/shipdensity_global.zip",
-        cutouts=expand(
+        cutouts=cutouts_for_extent_input,
-            "cutouts/" + CDIR + "{cutout}.nc",
-            cutout=[
-                config["renewable"][k]["cutout"]
-                for k in config["electricity"]["renewable_carriers"]
-            ],
-        ),
     output:
         RESOURCES + "shipdensity_raster.tif",
     log:

rules/build_sector.smk

@@ -150,7 +150,7 @@ rule build_heat_demands:
     params:
         snapshots=config["snapshots"],
     input:
-        pop_layout=RESOURCES + "pop_layout{weather_year}_{scope}.nc",
+        pop_layout=RESOURCES + "pop_layout_{scope}{weather_year}.nc",
         regions_onshore=RESOURCES
         + "regions_onshore_elec{weather_year}_s{simpl}_{clusters}.geojson",
         cutout="cutouts/" + CDIR + config["atlite"]["default_cutout"] + ".nc",
@@ -174,7 +174,7 @@ rule build_temperature_profiles:
     params:
         snapshots=config["snapshots"],
     input:
-        pop_layout=RESOURCES + "pop_layout{weather_year}_{scope}.nc",
+        pop_layout=RESOURCES + "pop_layout_{scope}{weather_year}.nc",
         regions_onshore=RESOURCES
         + "regions_onshore_elec{weather_year}_s{simpl}_{clusters}.geojson",
         cutout="cutouts/" + CDIR + config["atlite"]["default_cutout"] + ".nc",
@@ -246,7 +246,7 @@ rule build_solar_thermal_profiles:
         snapshots=config["snapshots"],
         solar_thermal=config["solar_thermal"],
     input:
-        pop_layout=RESOURCES + "pop_layout{weather_year}_{scope}.nc",
+        pop_layout=RESOURCES + "pop_layout_{scope}{weather_year}.nc",
         regions_onshore=RESOURCES
         + "regions_onshore_elec{weather_year}_s{simpl}_{clusters}.geojson",
         cutout="cutouts/" + CDIR + config["atlite"]["default_cutout"] + ".nc",
@@ -684,7 +684,8 @@ if config["sector"]["retrofitting"]["retro_endogen"]:
         input:
             building_stock="data/retro/data_building_stock.csv",
             data_tabula="data/bundle-sector/retro/tabula-calculator-calcsetbuilding.csv",
-            air_temperature=RESOURCES + "temp_air_total_elec{weather_year}_s{simpl}_{clusters}.nc",
+            air_temperature=RESOURCES
+            + "temp_air_total_elec{weather_year}_s{simpl}_{clusters}.nc",
             u_values_PL="data/retro/u_values_poland.csv",
             tax_w="data/retro/electricity_taxes_eu.csv",
             construction_index="data/retro/comparative_level_investment.csv",
@@ -719,7 +720,7 @@ if not config["sector"]["retrofitting"]["retro_endogen"]:
 
 rule build_population_weighted_energy_totals:
     input:
-        energy_totals=RESOURCES + "{kind}}_totals.csv",
+        energy_totals=RESOURCES + "{kind}_totals.csv",
         clustered_pop_layout=RESOURCES
         + "pop_layout_elec{weather_year}_s{simpl}_{clusters}.csv",
     output:
@@ -743,6 +744,8 @@ rule build_shipping_demand:
         regions=RESOURCES
         + "regions_onshore_elec{weather_year}_s{simpl}_{clusters}.geojson",
         demand=RESOURCES + "energy_totals.csv",
+    params:
+        energy_totals_year=config["energy"]["energy_totals_year"],
     output:
         RESOURCES + "shipping_demand{weather_year}_s{simpl}_{clusters}.csv",
     threads: 1

rules/postprocess.smk

@@ -28,7 +28,7 @@ if config["foresight"] != "perfect":
         benchmark:
             (
                 BENCHMARKS
-                + "plot_network/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}"
+                + "plot_network/elec{weather_year}_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}"
             )
         conda:
             "../envs/environment.yaml"
@@ -181,14 +181,14 @@ rule plot_elec_statistics:
         plotting=config["plotting"],
         barplots=STATISTICS_BARPLOTS,
     input:
-        network=RESULTS + "networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc",
+        network=RESULTS + "networks/elec{weather_year}_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc",
     output:
         **{
             f"{plot}_bar": RESULTS
-            + f"figures/statistics_{plot}_bar_elec_s{{simpl}}_{{clusters}}_ec_l{{ll}}_{{opts}}.pdf"
+            + f"figures/statistics_{plot}_bar_elec{{weather_year}}_s{{simpl}}_{{clusters}}_ec_l{{ll}}_{{opts}}.pdf"
             for plot in STATISTICS_BARPLOTS
         },
         barplots_touch=RESULTS
-        + "figures/.statistics_plots_elec_s{simpl}_{clusters}_ec_l{ll}_{opts}",
+        + "figures/.statistics_plots_elec{weather_year}_s{simpl}_{clusters}_ec_l{ll}_{opts}",
     script:
         "../scripts/plot_statistics.py"

rules/retrieve.smk

@@ -144,9 +144,6 @@ if config["enable"]["retrieve"] and config["enable"].get(
         protected(
             directory("data/bundle-sector/eurostat-energy_balances-june_2016_edition")
         ),
-        protected(
-            directory("data/bundle-sector/eurostat-energy_balances-june_2021_edition")
-        ),
         protected(
             directory("data/bundle-sector/eurostat-energy_balances-may_2018_edition")
         ),

scripts/base_network.py

@@ -558,7 +558,7 @@ def _set_countries_and_substations(n, config, country_shapes, offshore_shapes):
         ~buses["under_construction"]
     )
 
-    c_nan_b = buses.country == "na"
+    c_nan_b = buses.country.fillna("na") == "na"
     if c_nan_b.sum() > 0:
         c_tag = _get_country(buses.loc[c_nan_b])
         c_tag.loc[~c_tag.isin(countries)] = np.nan

scripts/build_energy_totals.py

@@ -356,7 +356,7 @@ def build_idees(countries, year=None):
     func = partial(idees_per_country, year=year, base_dir=snakemake.input.idees)
     tqdm_kwargs = dict(
         ascii=False,
-        unit="" country",
+        unit=" country",
         total=len(countries),
         desc="Build from IDEES database",
         disable=disable_progress
@@ -382,9 +382,6 @@ def build_idees(countries, year=None):
     total_heat = df[subset].sum(axis=1)
     df["district heat share"] = district_heat.div(total_heat)
 
-    if year:
-        df = df.xs(int(year), level='year')
-
     df.columns.name = 'item'
 
     return df
@@ -741,12 +738,13 @@ if __name__ == "__main__":
 
     countries = snakemake.params.countries
     idees_countries = pd.Index(countries).intersection(eu28)
+    countries_without_ch = pd.Index(countries).difference(["CH"])
 
     data_year = params["energy_totals_year"]
     report_year = snakemake.params.energy["eurostat_report_year"]
     input_eurostat = snakemake.input.eurostat
-    eurostat = build_eurostat(input_eurostat, countries, report_year, data_year)
+    eurostat = build_eurostat(countries_without_ch)
-    swiss = build_swiss(data_year)
+    swiss = build_swiss()
     idees = build_idees(idees_countries, data_year)
 
     energy = build_energy_totals(countries, eurostat, swiss, idees)

scripts/build_heat_totals.py

@@ -29,6 +29,14 @@ def approximate_heat_demand(energy_totals, hdd):
             Y = demand[c].dropna()
             X = hdd.loc[Y.index, c]
 
+            # Sometimes (looking at you, Switzerland) we only have
+            # _one_ year of heating data to base the prediction on. In
+            # this case we add a point at 0, 0 to make a "polynomial"
+            # fit work.
+            if len(X) == len(Y) == 1:
+                X.loc[-1] = 0
+                Y.loc[-1] = 0
+
             to_predict = hdd.index.difference(Y.index)
             X_pred = hdd.loc[to_predict, c]
 

scripts/build_hydro_profile.py

@@ -184,6 +184,8 @@ if __name__ == "__main__":
     fn = snakemake.input.eia_hydro_generation
     eia_stats = get_eia_annual_hydro_generation(fn, countries)
 
+    config_hydro = snakemake.config["renewable"]["hydro"]
+
     if config_hydro.get("eia_correct_by_capacity"):
         fn = snakemake.input.eia_hydro_capacity
         correct_eia_stats_by_capacity(eia_stats, fn, countries)

scripts/build_population_weighted_energy_totals.py

@@ -26,7 +26,7 @@ if __name__ == "__main__":
 
     pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0)
 
-    totals = pd.read_csv(snakemake.input.totals, index_col=[0, 1])
+    totals = pd.read_csv(snakemake.input.energy_totals, index_col=[0, 1])
     totals = totals.xs(data_year, level="year")
 
     nodal_totals = totals.loc[pop_layout.ct].fillna(0.0)

scripts/build_renewable_profiles.py

@@ -216,7 +216,15 @@ if __name__ == "__main__":
     else:
         client = None
 
-    sns = pd.date_range(freq="h", **snakemake.config["snapshots"])
+    weather_year = snakemake.wildcards.weather_year
+    if weather_year:
+        snapshots = dict(
+            start=weather_year, end=str(int(weather_year) + 1), inclusive="left"
+        )
+    else:
+        snapshots = snakemake.config["snapshots"]
+    sns = pd.date_range(freq="h", **snapshots)
+
     cutout = atlite.Cutout(snakemake.input.cutout).sel(time=sns)
     regions = gpd.read_file(snakemake.input.regions)
     assert not regions.empty, (

scripts/build_shipping_demand.py

@@ -24,9 +24,10 @@ if __name__ == "__main__":
 
     scope = gpd.read_file(snakemake.input.scope).geometry[0]
     regions = gpd.read_file(snakemake.input.regions).set_index("name")
-    demand = pd.read_csv(snakemake.input.demand, index_col=0)[
+    demand = pd.read_csv(snakemake.input.demand, index_col=[0, 1])[
         "total international navigation"
     ]
+    demand = demand.xs(snakemake.params.energy_totals_year, level=1)
 
     # read port data into GeoDataFrame
     with open(snakemake.input.ports, "r", encoding="latin_1") as f: