Add processing for shipping routes

README.md

@@ -149,6 +149,8 @@ This module is based on the following research and datasets:
     * License: "The GEBCO Grid is placed in the public domain and may be used free of charge. [...] Users must: Acknowledge the source of The GEBCO Grid. A suitable form of attribution is given in the documentation that accompanies The GEBCO Grid."
 * [GHSL (Global Human Settlement Layer)](https://human-settlement.emergency.copernicus.eu/download.php) built-up surface data (R2023, GHS-BUILT-S, 100m resolution)
     * License: "The GHSL has been produced by the EC JRC as open and free data. Reuse is authorised, provided the source is acknowledged."
+* [World Bank Global Shipping Traffic Density](https://datacatalog.worldbank.org/search/dataset/0037580/global-shipping-traffic-density), the number of AIS-recorded vessel positions per square kilometre accumulated over 2015–2021
+    * License: Creative Commons Attribution 4.0 International (CC BY 4.0)
 * [WDPA (World Database on Protected Areas)](https://www.protectedplanet.net/)
     * License: Non-commercial allowed. Citation: "UNEP-WCMC and IUCN (2025), Protected Planet: The World Database on Protected Areas (WDPA) and World Database on Other Effective Area-based Conservation Measures (WD-OECM) [Online], June 2025, Cambridge, UK: UNEP-WCMC and IUCN. Available at: www.protectedplanet.net."
 

config/config.yaml

@@ -72,6 +72,12 @@ techs:
         min: -50
         max: 0
         share: 0.8
+      # AIS vessel-position counts per km2 accumulated over 2015–2021.
+      ship_travel:
+        # PyPSA-Eur's default for high-density ship routes
+        # 400 AIS * 8760 hours * 6 years
+        min: 0
+        max: 21024000
     binary_layers:
       regions_land: 0
       regions_maritime: 1

workflow/internal/settings.yaml

@@ -10,6 +10,8 @@ resources:
     globcover_landcover_tif: "GLOBCOVER_L4_200901_200912_V2.3.tif"
     ghsl: "https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/GHSL/GHS_BUILT_S_GLOBE_R2023A/GHS_BUILT_S_E2025_GLOBE_R2023A_4326_30ss/V1-0/GHS_BUILT_S_E2025_GLOBE_R2023A_4326_30ss_V1_0.zip"
     ghsl_tif: "GHS_BUILT_S_E2025_GLOBE_R2023A_4326_30ss_V1_0.tif"
+    ship_travel: "https://datacatalogfiles.worldbank.org/ddh-published/0037580/5/DR0045406/shipdensity_global.zip"
+    ship_travel_tif: "shipdensity_global.tif"
 
 # NOTE: only works because GlobCover is already in this CRS!
 target_crs: "EPSG:4326"

workflow/rules/automatic.smk

@@ -160,6 +160,58 @@ rule clip_landcover:
         "v9.14.0/geo/rasterio/clip"
 
 
+##
+# Global Ship Traffic Density
+##
+
+
+rule download_ship_travel:
+    output:
+        "<resources>/automatic/global/ship_travel_density.zip",
+    log:
+        "<logs>/download_ship_travel.log",
+    localrule: True
+    conda:
+        "../envs/module.yaml"
+    params:
+        url=internal["resources"]["automatic"]["ship_travel"],
+    message:
+        "Download Global Ship Density for all vessel types."
+    shell:
+        """
+        curl -sSLo {output:q} {params.url:q} >{log:q} 2>&1
+        """
+
+
+rule unzip_ship_travel:
+    input:
+        rules.download_ship_travel.output[0],
+    output:
+        temp("<resources>/automatic/global/ship_travel.tif"),
+    log:
+        "<logs>/unzip_ship_travel.log",
+    params:
+        internal_paths=internal["resources"]["automatic"]["ship_travel_tif"],
+    message:
+        "Unzip the relevant TIF file from the ship travel density data."
+    wrapper:
+        "v9.8.0/utils/libarchive/extract"
+
+
+rule clip_ship_travel:
+    input:
+        like_vector=rules.normalise_shapes.output.shapes,
+        raster=rules.unzip_ship_travel.output[0],
+    output:
+        path="<resources>/automatic/cutout/{shape}/ship_travel.tif",
+    log:
+        "<logs>/{shape}/clip_ship_travel.log",
+    message:
+        "Cut ship travel data to the bounds of the input shapefile."
+    wrapper:
+        "v9.14.0/geo/rasterio/clip"
+
+
 ##
 # Global Human Settlement Layer (GHSL)
 ##

workflow/rules/functions.smk

@@ -6,3 +6,13 @@ def get_subunits(wildcards):
             os.path.join(checkpoint_output, "{subunit}.parquet")
         ).subunit,
     )
+
+
+def uses_ship_travel(wildcards):
+    """Return whether a base tech or this subunit's overrides use ship travel."""
+    tech_configs = list(config.get("techs", {}).values())
+    tech_configs.extend(config.get("overrides", {}).get(wildcards.subunit, {}).values())
+    return any(
+        "ship_travel" in tech_config.get("continuous_layers", {})
+        for tech_config in tech_configs
+    )

workflow/rules/process.smk

@@ -1,3 +1,5 @@
+import shlex
+
 
 rule prepare_resampled_inputs:
     input:
@@ -8,6 +10,11 @@ rule prepare_resampled_inputs:
         settlement_path=rules.clip_settlement.output,
         bathymetry_path=rules.clip_bathymetry.output,
         protected_area_path=rules.rasterise_clip_wdpa.output,
+        ship_travel_path=branch(
+            condition=uses_ship_travel,
+            then=rules.clip_ship_travel.output,
+            otherwise=[],
+        ),
     output:
         resampled_input="<resources>/automatic/resampled_inputs/{shape}/{subunit}.nc",
         plot=report(
@@ -22,6 +29,11 @@ rule prepare_resampled_inputs:
         # Use internal defaults if not overridden
         land_cover_types_yaml_string=internal["land_cover_types"]
         | config.get("land_cover_types", {}),
+        ship_travel_arg=lambda wildcards, input: (
+            f"--ship-travel-path {shlex.quote(str(input.ship_travel_path))}"
+            if input.ship_travel_path
+            else ""
+        ),
     message:
         "Resample inputs for {wildcards.subunit} in {wildcards.shape} to the projection and resolution of the land cover data, while aggregating land cover types."
     shell:
@@ -30,7 +42,8 @@ rule prepare_resampled_inputs:
             "{input.shapes}/{wildcards.subunit}.parquet" \
             {input.land_cover_path:q} {input.slope_path:q} {input.settlement_path:q} {input.bathymetry_path:q} {input.protected_area_path:q} \
             {params.land_cover_types_yaml_string:q} \
-            {output.resampled_input:q} {output.plot:q} >{log:q} 2>&1
+            {output.resampled_input:q} {output.plot:q} \
+            {params.ship_travel_arg} >{log:q} 2>&1
         """
 
 

workflow/scripts/resample.py

@@ -139,6 +139,7 @@ def _rasterize_regions(shapes, reference_raster):
 @click.argument("land_cover_configuration_yaml_string", type=str)
 @click.argument("output_path", type=str)
 @click.argument("plot_path", type=str)
+@click.option("--ship-travel-path", type=str)
 def resample_inputs(
     shapes_path,
     land_cover_path,
@@ -149,6 +150,7 @@ def resample_inputs(
     land_cover_configuration_yaml_string,
     output_path,
     plot_path,
+    ship_travel_path,
 ):
     """Resample various geospatial datasets to a common shape and resolution.
 
@@ -274,6 +276,17 @@ def resample_inputs(
     )
     del protected_areas
 
+    ##
+    # Global shipping traffic density (AIS position counts per square kilometre)
+    ##
+    if ship_travel_path:
+        ship_travel = rxr.open_rasterio(ship_travel_path, masked=True)
+        print(f"Ship travel resolution: {ship_travel.rio.resolution()}")
+        resampled["ship_travel"] = ship_travel.rio.reproject_match(
+            reference_raster, resampling=Resampling.average
+        )
+        del ship_travel
+
     netcdf4_encoding = {
         var: {"zlib": True, "complevel": 1}
         for var in resampled.data_vars