Skip to content

Rule build_natura_raster

Workflow Diagram

See the complete workflow in the repository.

digraph snakemake_dag { graph [bgcolor=white, margin=0, size="8,5" ]; node [fontname=sans, fontsize=10, penwidth=2, shape=box, style=rounded ]; edge [color=grey, penwidth=2 ]; 9 [color="0.22 0.6 0.85", label=build_renewable_profiles]; 12 [color="0.31 0.6 0.85", fillcolor=gray, label=build_natura_raster, style=filled]; 12 -> 9; }

Script Documentation

Builds a rasterized Natura and environmental exclusion mask from vector geometries. The geometries are sourced from Protected Planet Data and filtered for the required regions.

Relevant Settings

countries:

enable:
    build_natura_raster:
    progress_bar:

crs:
    area_crs:

natura:
    natura_size:
    natura_resolution:
    window_size:
    buffer_size:

renewable:
    {technology}:
        cutout:

Inputs

  • data/landcover/world_protected_areas/*.shp: Vectorized shapefiles representing the world protected areas from Protected Planet Data.
  • cutouts/{CDIR}/{cutout}.nc: Atlite cutout specified in renewable: {technology}: cutout for each technology. The cutout(s) are used to determine geographic coverage and resolution alignment.
  • resources/{RDIR}/shapes/country_shapes.geojson: Onshore country geometries used to determine the spatial extent of the rasterization.
  • resources/{RDIR}/shapes/offshore_shapes.geojson: Offshore regions associated with the selected countries.

Outputs

  • resources/{RDIR}/natura.tiff: Rasterized version of the world protected areas.

Description

The rule build_natura_raster converts large collections of vector-based environmental exclusion geometries into a rasterized binary mask containing:

  • 1 for raster cells intersecting exclusion geometries, and
  • 0 elsewhere.

First, the rule determines the spatial extent of the rasterization which is configured through natura: natura_size and can be based on:

  • the full global extent,
  • the combined extent of all configured renewable technology cutouts (renewable: {technology}: cutout), or
  • the combined extent of all selected country (countries) and their offshore regions.

In case of the countries extent, the country and offshore geometries are reprojected into (crs: area_crs). Additionally, a buffer (natura: buffer_size) is applied to the geometries to ensure that all regions of interest are fully included in the rasterization.

The rasterization is performed as follows:

  • an empty GeoTIFF is created using the extent determined above and a grid resolution defined by natura: natura_resolution,
  • the extent is split into smaller windows with a maximum width and height of natura: window_size,
  • each input .shp file is processed independently and rasterized window by window.

This windowed processing strategy reduces memory usage and allows the rule to run efficiently on standard compute environments.

Notes

This script only runs in the current PyPSA-Earth workflow if enable: build_natura_raster is true. Otherwise, the workflow will use the general data/natura/natura.tiff file, which is copied using the rule copy_defaultnatura_tiff.

There are two main differences between the two options, the data source and the license:

get_relevant_regions(country_shapes, offshore_shapes, natura_crs, buffer)

Load and merge country and offshore regions into a unified geometry.

The resulting geometry is buffered to ensure all relevant nearby regions are included.

Parameters

country_shapes : str Path to the vector file containing the country geometries. offshore_shapes : str Path to the vector file containing the offshore geometries. natura_crs : str Coordinate reference system used for all geometries. buffer : float Buffer distance applied to both country and offshore geometries. Units are determined by natura_crs.

Returns

gpd.GeoDataFrame GeoDataFrame containing a single merged geometry representing the buffered union of all country and offshore regions.

Notes

Both input datasets are reprojected to natura_crs before merging.

Examples

regions = get_relevant_regions(
    "resources/shapes/country_shapes.geojson",
    "resources/shapes/offshore_shapes.geojson",
    "ESRI:54009",
    10000,
)
len(regions)
# 1
regions.crs.to_string()
# "ESRI:54009"

get_fileshapes(list_paths, accepted_formats=('.shp',))

Function to parse the list of paths and identify the ones with one of the accepted file formats.

Parameters

list_paths : list[str] List of paths to check. accepted_formats : str | tuple[str, ...], optional File format(s) to accepted. If a string is provided, only that extension is accepted. If a tuple is provided, any of the extensions are accepted. Default is (".shp",).

Returns

list[str] List of all paths which are of one of the accepted formats.

Examples

paths = ["shape_file.shp", "shape_index.shi"]
get_fileshapes(paths)
# ["shape_file.shp"]
get_fileshapes(paths, ".shi")
# ["shape_index.shi"]

determine_region_xXyY(cutout_name, regions, natura_size, out_logging)

Determine the bounds of the analyzed regions.

Parameters

cutout_name : str Path of the cutout file. regions : gpd.GeoDataFrame | None GeoDataFrame containing the analyzed regions. Required if natura_size is "countries", otherwise ignored. natura_size : str Flag to determine which region should be used. "global" includes the entire world, "cutout" the extent of the cutout, and "countries" only includes the bounds of the requested countries and their offshore regions. out_logging : bool If True, emits progress information via the module logger.

Returns

list[float] Bounding box of the region in the format: [min_lon, max_lon, min_lat, max_lat].

Examples

cutout_path = "cutouts/cutout-2013-era5.nc"
regions = None

# Global extent:
determine_region_xXyY(cutout_path, regions, "global", False)
# [-180, 180, -90, 90]

# Cutout extent (Africa):
determine_region_xXyY(cutout_path, regions, "cutout", False)
# [-19.8, 67.8, -37.8, 39.6]

# Countries extent:
import geopandas as gpd
from shapely.geometry import Polygon

regions = gpd.GeoDataFrame(
    geometry=[Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])],
    crs="EPSG:4326",
)
determine_region_xXyY(cutout_path, regions, "countries", False)
# [0.0, 10.0, 0.0, 10.0]

get_transform_and_shape(bounds, res, out_logging)

Compute an affine transform and raster shape from spatial bounds and resolution.

Parameters

bounds : list[float] Bounding box in the format: [min_lon, min_lat, max_lon, max_lat]. res : float Spatial resolution of the grid (in coordinate units, e.g. degrees). out_logging : bool If True, emits progress information via the module logger.

Returns

transform : rio.Affine Affine transform mapping raster indices (row, col) to spatial coordinates. shape : tuple[int, int] Raster shape as (n_lat, n_lon), corresponding to (rows, cols).

Examples

import rasterio as rio

# [min_lon, min_lat, max_lon, max_lat]
bounds = [0.0, 50.0, 3.0, 52.0]
res = 1.0

transform, shape = get_transform_and_shape(bounds, res, False)
shape
# (2, 3)
transform
# Affine(1.0, 0.0, 0.0,
# 0.0, -1.0, 52.0)

decide_bigtiff_flag(out_shape, dtype='uint8', safety_factor=1.1)

Decide whether a raster requires BIGTIFF storage based on raster shape.

BIGTIFF is required for GeoTIFF files larger than approximately 4 GB.

Parameters

out_shape : tuple[int, int] Raster shape as (n_rows, n_cols). dtype : str, optional Data type of the raster values. Default is "uint8". safety_factor : float, optional Multiplicative buffer applied to the estimated size to account for metadata, compression inefficiencies, and file overhead. Default is 1.1.

Returns

str "YES" if the estimated size is larger than the BIGTIFF threshold, otherwise "NO".

Notes

The size estimate is computed as:

n_rows * n_cols * bytes_per_pixel * safety_factor

The BIGTIFF threshold is fixed at 4,000,000,000 bytes (~4 GB).

Examples

decide_bigtiff_flag((100, 100), dtype="uint8")
# "NO"