Merge pull request #401 from p-glaum/shipping_raster

Consider shipping routes for offshore potential calculation

Snakefile

@@ -185,6 +185,22 @@ if config['enable'].get('retrieve_natura_raster', True):
         run: move(input[0], output[0])
 
 
+rule retrieve_ship_raster:
+    input: HTTP.remote("https://zenodo.org/record/6953563/files/shipdensity_global.zip", keep_local=True, static=True)
+    output: "data/shipdensity_global.zip"
+    run: move(input[0], output[0])
+
+
+rule build_ship_raster:
+    input:
+        ship_density="data/shipdensity_global.zip",
+        cutouts=expand("cutouts/{cutouts}.nc", **config['atlite'])
+    output: "resources/shipdensity_raster.nc"
+    log: "logs/build_ship_raster.log"
+    benchmark: "benchmarks/build_ship_raster"
+    script: "scripts/build_ship_raster.py"
+
+
 rule build_renewable_profiles:
     input:
         base_network="networks/base.nc",
@@ -195,6 +211,9 @@ rule build_renewable_profiles:
         gebco=lambda w: ("data/bundle/GEBCO_2014_2D.nc"
                          if "max_depth" in config["renewable"][w.technology].keys()
                          else []),
+        ship_density= lambda w: ("resources/shipdensity_raster.nc"
+                         if "ship_threshold" in config["renewable"][w.technology].keys()
+                         else []),
         country_shapes='resources/country_shapes.geojson',
         offshore_shapes='resources/offshore_shapes.geojson',
         regions=lambda w: ("resources/regions_onshore.geojson"

config.default.yaml

@@ -135,6 +135,7 @@ renewable:
     # until done more rigorously in #153
     corine: [44, 255]
     natura: true
+    ship_threshold: 400
     max_depth: 50
     max_shore_distance: 30000
     excluder_resolution: 200
@@ -153,6 +154,7 @@ renewable:
     # until done more rigorously in #153
     corine: [44, 255]
     natura: true
+    ship_threshold: 400
     max_depth: 50
     min_shore_distance: 30000
     excluder_resolution: 200

config.tutorial.yaml

@@ -88,6 +88,7 @@ renewable:
     # correction_factor: 0.93
     corine: [44, 255]
     natura: true
+    ship_threshold: 400
     max_shore_distance: 30000
     excluder_resolution: 200
     potential: simple # or conservative
@@ -102,6 +103,7 @@ renewable:
     # correction_factor: 0.93
     corine: [44, 255]
     natura: true
+    ship_threshold: 400
     min_shore_distance: 30000
     excluder_resolution: 200
     potential: simple # or conservative

doc/configtables/offwind-ac.csv

@@ -6,6 +6,7 @@ resource,,,
 capacity_per_sqkm,:math:`MW/km^2`,float,"Allowable density of wind turbine placement."
 corine,--,"Any *realistic* subset of the `CORINE Land Cover code list <http://www.eea.europa.eu/data-and-maps/data/corine-land-cover-2006-raster-1/corine-land-cover-classes-and/clc_legend.csv/at_download/file>`_","Specifies areas according to CORINE Land Cover codes which are generally eligible for AC-connected offshore wind turbine placement."
 natura,bool,"{true, false}","Switch to exclude `Natura 2000 <https://en.wikipedia.org/wiki/Natura_2000>`_ natural protection areas. Area is excluded if ``true``."
+ship_threshold,--,float,"Ship density threshold from which areas are excluded."
 max_depth,m,float,"Maximum sea water depth at which wind turbines can be build. Maritime areas with deeper waters are excluded in the process of calculating the AC-connected offshore wind potential."
 min_shore_distance,m,float,"Minimum distance to the shore below which wind turbines cannot be build. Such areas close to the shore are excluded in the process of calculating the AC-connected offshore wind potential."
 potential,--,"One of {'simple', 'conservative'}","Method to compute the maximal installable potential for a node; confer :ref:`renewableprofiles`"

doc/configtables/offwind-dc.csv

@@ -6,6 +6,7 @@ resource,,,
 capacity_per_sqkm,:math:`MW/km^2`,float,"Allowable density of wind turbine placement."
 corine,--,"Any *realistic* subset of the `CORINE Land Cover code list <http://www.eea.europa.eu/data-and-maps/data/corine-land-cover-2006-raster-1/corine-land-cover-classes-and/clc_legend.csv/at_download/file>`_","Specifies areas according to CORINE Land Cover codes which are generally eligible for AC-connected offshore wind turbine placement."
 natura,bool,"{true, false}","Switch to exclude `Natura 2000 <https://en.wikipedia.org/wiki/Natura_2000>`_ natural protection areas. Area is excluded if ``true``."
+ship_threshold,--,float,"Ship density threshold from which areas are excluded."
 max_depth,m,float,"Maximum sea water depth at which wind turbines can be build. Maritime areas with deeper waters are excluded in the process of calculating the AC-connected offshore wind potential."
 min_shore_distance,m,float,"Minimum distance to the shore below which wind turbines cannot be build. Such areas close to the shore are excluded in the process of calculating the AC-connected offshore wind potential."
 potential,--,"One of {'simple', 'conservative'}","Method to compute the maximal installable potential for a node; confer :ref:`renewableprofiles`"

doc/release_notes.rst

@@ -10,7 +10,7 @@ Release Notes
 Upcoming Release
 ================
 
-* new feature
+* Add functionality to consider shipping routes when calculating the available area for offshore technologies. Data for the shipping density comes from the `Global Shipping Traffic Density dataset <https://datacatalog.worldbank.org/search/dataset/0037580/Global-Shipping-Traffic-Density>`
 
 PyPSA-Eur 0.5.0 (27th July 2022)
 =====================================

scripts/build_renewable_profiles.py

@@ -243,6 +243,11 @@ if __name__ == '__main__':
         buffer = corine["distance"]
         excluder.add_raster(snakemake.input.corine, codes=codes, buffer=buffer, crs=3035)
     
+    if "ship_threshold" in config:
+        shipping_threshold=config["ship_threshold"] * 8760 * 6 # approximation because 6 years of data which is hourly collected
+        func = functools.partial(np.less, shipping_threshold)
+        excluder.add_raster(snakemake.input.ship_density, codes=func, crs=4326, allow_no_overlap=True)
+
     if "max_depth" in config:
         # lambda not supported for atlite + multiprocessing
         # use named function np.greater with partially frozen argument instead

scripts/build_ship_raster.py

@@ -0,0 +1,63 @@
+# SPDX-FileCopyrightText: : 2022 The PyPSA-Eur Authors
+#
+# SPDX-License-Identifier: MIT
+
+"""
+Transforms the global ship density data from https://datacatalog.worldbank.org/search/dataset/0037580/Global-Shipping-Traffic-Density to the size of the considered cutout. The global ship density raster is later used for the exclusion when calculating the offshore potentials.
+
+Relevant Settings
+-----------------
+
+.. code:: yaml
+
+    renewable:
+        {technology}:
+            cutout:
+
+.. seealso::
+    Documentation of the configuration file ``config.yaml`` at
+    :ref:`renewable_cf`
+
+Inputs
+------
+
+- ``data/bundle/shipdensity/shipdensity_global.zip``: `Global ship density from <https://datacatalog.worldbank.org/search/dataset/0037580/Global-Shipping-Traffic-Density>`.
+
+Outputs
+-------
+
+- ``resources/europe_shipdensity_raster.nc``: Reduced version of `Global ship density from <https://datacatalog.worldbank.org/search/dataset/0037580/` to reduce computation time.
+
+Description
+-----------
+
+"""
+
+import logging
+from _helpers import configure_logging
+from build_natura_raster import determine_cutout_xXyY
+
+import zipfile
+import xarray as xr
+import os
+
+logger = logging.getLogger(__name__)
+
+if __name__ == "__main__":
+    if 'snakemake' not in globals():
+        from _helpers import mock_snakemake
+        snakemake = mock_snakemake('build_ship_raster')
+    configure_logging(snakemake)
+
+    cutouts = snakemake.input.cutouts
+    xs, Xs, ys, Ys = zip(*(determine_cutout_xXyY(cutout) for cutout in cutouts))
+
+    with zipfile.ZipFile(snakemake.input.ship_density) as zip_f:
+        zip_f.extract("shipdensity_global.tif")
+        ship_density = xr.open_rasterio("shipdensity_global.tif")
+        os.remove("shipdensity_global.tif") 
+
+    ship_density = ship_density.drop(["band"]).sel(x=slice(min(xs),max(Xs)), y=slice(max(Ys),min(ys)))
+
+    ship_density.to_netcdf(snakemake.output[0])
+

test/config.test1.yaml

@@ -87,6 +87,7 @@ renewable:
     # correction_factor: 0.93
     corine: [44, 255]
     natura: true
+    ship_threshold: 400
     max_shore_distance: 30000
     excluder_resolution: 200
     potential: simple # or conservative
@@ -101,6 +102,7 @@ renewable:
     # correction_factor: 0.93
     corine: [44, 255]
     natura: true
+    ship_threshold: 400
     min_shore_distance: 30000
     excluder_resolution: 200
     potential: simple # or conservative