UA-MD availability matrix: consider WDPA protection areas as substitute for Natura

rules/build_electricity.smk

@@ -210,6 +210,10 @@ rule determine_availability_matrix_MD_UA:
     input:
         copernicus=RESOURCES
         + "Copernicus_LC100_global_v3.0.1_2019-nrt_Discrete-Classification-map_EPSG-4326.tif",
+        wdpa=RESOURCES
+        + f"WDPA_{bYYYY}.gpkg",
+        wdpa_marine=RESOURCES
+        + f"WDPA_WDOECM_{bYYYY}_marine.gpkg",
         gebco=lambda w: (
             "data/bundle/GEBCO_2014_2D.nc"
             if "max_depth" in config["renewable"][w.technology].keys()

scripts/determine_availability_matrix_MD_UA.py

@@ -8,12 +8,18 @@ import logging
 import time
 
 import atlite
+import fiona
 import geopandas as gpd
 import numpy as np
 from _helpers import configure_logging
 
 logger = logging.getLogger(__name__)
 
+def get_wdpa_layer_name(wdpa_fn, layer_substring):
+    """Get layername from file "wdpa_fn" whose name contains "layer_substring"."""
+    l = fiona.listlayers(wdpa_fn)
+    return [_ for _ in l if layer_substring in _][0]
+
 if __name__ == "__main__":
     if "snakemake" not in globals():
         from _helpers import mock_snakemake
@@ -64,6 +70,29 @@ if __name__ == "__main__":
             snakemake.input.copernicus, codes=codes, buffer=buffer, crs="EPSG:4326"
         )
 
+    if config["natura"]:
+        wdpa_fn = snakemake.input.wdpa_marine if "offwind" in snakemake.wildcards.technology else snakemake.input.wdpa
+        layer = get_wdpa_layer_name(wdpa_fn, "polygons")
+        wdpa = gpd.read_file(
+            wdpa_fn,
+            bbox=regions.geometry,
+            layer=layer,
+        ).to_crs(3035)
+        if not wdpa.empty:
+            excluder.add_geometry(wdpa.geometry)
+
+        layer = get_wdpa_layer_name(wdpa_fn, "points")
+        wdpa_pts = gpd.read_file(
+            wdpa_fn,
+            bbox=regions.geometry,
+            layer=layer,
+        ).to_crs(3035)
+        wdpa_pts = wdpa_pts[wdpa_pts["REP_AREA"] > 1]
+        wdpa_pts["buffer_radius"] = np.sqrt(wdpa_pts["REP_AREA"] / np.pi) * 1000
+        wdpa_pts = wdpa_pts.set_geometry(wdpa_pts["geometry"].buffer(wdpa_pts["buffer_radius"]))
+        if not wdpa_pts.empty:
+            excluder.add_geometry(wdpa_pts.geometry)
+
     if "max_depth" in config:
         # lambda not supported for atlite + multiprocessing
         # use named function np.greater with partially frozen argument instead