Flaechenziel for Onshore Wind

config/config.yaml

@@ -4,7 +4,7 @@
 
 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#run
 run:
-  prefix: 20241121-fix-missing-gas-capa
+  prefix: 20241202-adjust-onwind
 
   name:
   # - CurrentPolicies
@@ -62,7 +62,7 @@ scenario:
   ll:
   - vopt
   clusters:
-  - 27 #current options: 27, 49
+  - 49 #current options: 27, 49
   opts:
   - ''
   sector_opts:
@@ -206,38 +206,38 @@ clustering:
   # print(fw.div(fw.sum()).subtract(5e-5).round(4).to_dict().__repr__().replace(",","\n"))
   focus_weights:
   # 27 nodes: 8 for Germany, 3 for Italy, 2 each for Denmark, UK and Spain, 1 per each of other 10 "Stromnachbarn"
-    'DE': 0.2966
-    'AT': 0.0370
-    'BE': 0.0370
-    'CH': 0.0370
-    'CZ': 0.0370
-    'DK': 0.0741
-    'FR': 0.0370
-    'GB': 0.0741
-    'LU': 0.0370
-    'NL': 0.0370
-    'NO': 0.0370
-    'PL': 0.0370
-    'SE': 0.0370
-    'ES': 0.0741
-    'IT': 0.1111
+    # 'DE': 0.2966
+    # 'AT': 0.0370
+    # 'BE': 0.0370
+    # 'CH': 0.0370
+    # 'CZ': 0.0370
+    # 'DK': 0.0741
+    # 'FR': 0.0370
+    # 'GB': 0.0741
+    # 'LU': 0.0370
+    # 'NL': 0.0370
+    # 'NO': 0.0370
+    # 'PL': 0.0370
+    # 'SE': 0.0370
+    # 'ES': 0.0741
+    # 'IT': 0.1111
   # high spatial resolution: change clusters to 49
   # 49 nodes: 30 for Germany, 3 for Italy, 2 each for Denmark, UK and Spain, 1 per each of other 10 "Stromnachbarn"
-    # 'DE': 0.6124
-    # 'AT': 0.0204
-    # 'BE': 0.0204
-    # 'CH': 0.0204
-    # 'CZ': 0.0204
-    # 'DK': 0.0408
-    # 'FR': 0.0204
-    # 'GB': 0.0408
-    # 'LU': 0.0204
-    # 'NL': 0.0204
-    # 'NO': 0.0204
-    # 'PL': 0.0204
-    # 'SE': 0.0204
-    # 'ES': 0.0408
-    # 'IT': 0.0612
+    'DE': 0.6124
+    'AT': 0.0204
+    'BE': 0.0204
+    'CH': 0.0204
+    'CZ': 0.0204
+    'DK': 0.0408
+    'FR': 0.0204
+    'GB': 0.0408
+    'LU': 0.0204
+    'NL': 0.0204
+    'NO': 0.0204
+    'PL': 0.0204
+    'SE': 0.0204
+    'ES': 0.0408
+    'IT': 0.0612
   temporal:
     resolution_sector: 365H
 
@@ -432,6 +432,12 @@ solving:
         H2 pipeline retrofitted: 0.05
       fractional_last_unit_size: true
   constraints:
+    onwind_flaechenziel:
+      enable: true
+      percentage: # % of total land area
+        2035: 0.5
+        2040: 0.75
+        2045: 1.0
     limits_capacity_max:
       Generator:
         onwind:

workflow/Snakefile

@@ -311,6 +311,7 @@ use rule solve_sector_network_myopic from pypsaeur with:
         network=RESULTS
         + "prenetworks-final/base_s_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc",
         co2_totals_name=resources("co2_totals.csv"),
+        regions_onshore=resources("regions_onshore_base_s_{clusters}.geojson"),
 
 
 rule modify_existing_heating:

workflow/scripts/additional_functionality.py

@@ -2,6 +2,7 @@
 
 import logging
 
+import geopandas as gpd
 import pandas as pd
 from prepare_sector_network import determine_emission_sectors
 from xarray import DataArray
@@ -663,6 +664,109 @@ def adapt_nuclear_output(n):
     )
 
 
+def add_onwind_constraint(n, investment_year, snakemake, onwind_constraint):
+
+    regions_onshore = gpd.read_file(snakemake.input.regions_onshore)
+    regions_onshore.set_index("name", inplace=True)
+    de_regions_onshore = regions_onshore[regions_onshore.index.str.startswith("DE")]
+    de_regions_onshore["centroid_lat"] = de_regions_onshore.geometry.centroid.y
+    de_regions_onshore["centroid_lon"] = de_regions_onshore.geometry.centroid.x
+
+    # get index that centroid_lat is smaller < 50
+    bins = [[], [], [], [], []]
+    # south by
+    bins[0] = de_regions_onshore[
+        (de_regions_onshore["centroid_lat"] < 50)
+        & (de_regions_onshore["centroid_lon"] > 10)
+    ].index.to_list()
+    # bw
+    bins[1] = de_regions_onshore[
+        (de_regions_onshore["centroid_lat"] < 49)
+        & (de_regions_onshore["centroid_lon"] <= 10)
+    ].index.to_list()
+    # sl, he, rp, nrw
+    bins[2] = de_regions_onshore[
+        (de_regions_onshore["centroid_lat"] >= 49)
+        & (de_regions_onshore["centroid_lat"] < 52)
+        & (de_regions_onshore["centroid_lon"] <= 10)
+    ].index.to_list()
+    # north by, th, sa
+    bins[3] = de_regions_onshore[
+        (de_regions_onshore["centroid_lat"] >= 50)
+        & (de_regions_onshore["centroid_lat"] < 52)
+        & (de_regions_onshore["centroid_lon"] > 10)
+    ].index.to_list()
+    bins[4] = de_regions_onshore[
+        (de_regions_onshore["centroid_lat"] >= 52)
+    ].index.to_list()
+
+    regions_onshore = regions_onshore.to_crs(epsg=3035)
+    # area in sqkm
+    area = (
+        regions_onshore[regions_onshore.index.str.startswith("DE")].geometry.area / 1e6
+    )
+
+    mw_per_sqkm = 30  # source: https://www.ffe.de/wp-content/uploads/2022/02/FfE-Discussion-Paper-2-der-Landesflaeche-fuer-Windenergie-ein-geeignetes-Mass.pdf page 7
+    # % into fraction
+    flaechenziel = onwind_constraint[investment_year] / 100
+
+    for entry in bins:
+        if not entry:
+            continue
+        area_region = area[entry]
+        limit = flaechenziel * area_region.sum() * mw_per_sqkm
+
+        valid_components = (n.generators.bus.isin(entry)) & (
+            n.generators.carrier == "onwind"
+        )
+
+        existing_index = n.generators.index[
+            valid_components & ~n.generators.p_nom_extendable
+        ]
+        extendable_index = n.generators.index[
+            valid_components & n.generators.p_nom_extendable
+        ]
+
+        if n.generators.loc[extendable_index, "p_nom_max"].sum() < 10:
+            logger.warning(
+                f"Skipping onwind constraint for region consisting of {entry} as no extendable capacity exists"
+            )
+            continue
+        elif existing_index.empty:
+            logger.info(f"No existing onwind capacity in region consisting of {entry}")
+            existing_cap = 0
+        else:
+            existing_cap = n.generators.loc[existing_index, "p_nom"].sum()
+
+        lhs = n.model["Generator-p_nom"].loc[extendable_index].sum()
+        rhs = limit - existing_cap
+        if rhs <= 0:
+            logger.info(
+                f"Existing capacity {existing_cap/1000} GW fullfills the Flaechenziel of {limit/1000} GW in {entry}."
+            )
+            continue
+
+        cname = f"flaechenziel_onwind-{entry}"
+
+        logger.info(f"Adding Flaechenziel for {entry} of {rhs/1000} GW.")
+        n.model.add_constraints(lhs >= rhs, name=f"GlobalConstraint-{cname}")
+
+        if cname in n.global_constraints.index:
+            logger.warning(
+                f"Global constraint {cname} already exists. Dropping and adding it again."
+            )
+            n.global_constraints.drop(cname, inplace=True)
+
+        n.add(
+            "GlobalConstraint",
+            cname,
+            constant=rhs,
+            sense=">=",
+            type="",
+            carrier_attribute="",
+        )
+
+
 def additional_functionality(n, snapshots, snakemake):
 
     logger.info("Adding Ariadne-specific functionality")
@@ -712,3 +816,14 @@ def additional_functionality(n, snapshots, snakemake):
 
     if investment_year == 2020:
         adapt_nuclear_output(n)
+
+    if (
+        constraints["onwind_flaechenziel"]["enable"]
+        and investment_year in constraints["onwind_flaechenziel"]["percentage"].keys()
+    ):
+        add_onwind_constraint(
+            n,
+            investment_year,
+            snakemake,
+            constraints["onwind_flaechenziel"]["percentage"],
+        )

workflow/scripts/build_scenarios.py

@@ -178,8 +178,9 @@ def write_to_scenario_yaml(input, output, scenarios, df):
             df.loc[:, fallback_reference_scenario, :], planning_horizons
         )
 
-
-        if reference_scenario.startswith("KN2045plus"): # Still waiting for REMIND uploads
+        if reference_scenario.startswith(
+            "KN2045plus"
+        ):  # Still waiting for REMIND uploads
             fallback_reference_scenario = reference_scenario
 
         co2_budget_source = config[scenario]["co2_budget_DE_source"]