From 3c53747f4cdd82a9aa2af2924785056edb53f20a Mon Sep 17 00:00:00 2001 From: Ivan Ruiz Manuel <72193617+irm-codebase@users.noreply.github.com> Date: Wed, 12 Aug 2026 21:49:06 +0200 Subject: [PATCH 1/5] Update clip wrapper and replace clip script --- README.md | 10 ++-- config/config.yaml | 6 +-- tests/integration/test_config.yaml | 4 +- workflow/internal/settings.yaml | 8 +-- workflow/rules/automatic.smk | 85 +++++++++++------------------- workflow/rules/process.smk | 8 +-- workflow/scripts/clip_raster.py | 42 --------------- workflow/scripts/report.py | 3 ++ workflow/scripts/resample.py | 4 +- workflow/scripts/tif_to_png.py | 3 ++ 10 files changed, 58 insertions(+), 115 deletions(-) delete mode 100644 workflow/scripts/clip_raster.py diff --git a/README.md b/README.md index 63f2b79..416a7cd 100644 --- a/README.md +++ b/README.md @@ -70,12 +70,12 @@ Here is a `wind_offshore` example. We start with the `pixel_area`, the total sur wind_onshore: initial_area: pixel_area continuous_layers: - slope: - min: 0 - max: 20 + slope_deg: + min: 0 + max: 20 settlement_share: - min: 0 - max: 0.01 + min: 0 + max: 0.01 binary_layers: regions_maritime: 0 regions_land: 1 diff --git a/config/config.yaml b/config/config.yaml index 947d727..9f5807a 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -2,7 +2,7 @@ # - "UTM": project each shape to the UTM zone of its centroid for buffering # - "epsg:xxxx": use the specified CRS for all buffering # A good option is "epsg:8857" (WGS 84 / Equal Earth Greenwich) for global coverage -buffer_crs: "epsg:8857" +buffer_crs: "epsg:3035" split_by: country_id # likely country_id or shape_id @@ -30,7 +30,7 @@ techs: pv_open_field: initial_area: pixel_area continuous_layers: - slope: + slope_deg: min: 0 max: 3 settlement_share: @@ -49,7 +49,7 @@ techs: wind_onshore: initial_area: pixel_area continuous_layers: - slope: + slope_deg: min: 0 max: 20 settlement_share: diff --git a/tests/integration/test_config.yaml b/tests/integration/test_config.yaml index 3a579f4..b541d41 100644 --- a/tests/integration/test_config.yaml +++ b/tests/integration/test_config.yaml @@ -24,7 +24,7 @@ module_area_potentials: pv_open_field: initial_area: pixel_area continuous_layers: - slope: + slope_deg: min: 0 max: 3 settlement_share: @@ -43,7 +43,7 @@ module_area_potentials: wind_onshore: initial_area: pixel_area continuous_layers: - slope: + slope_deg: min: 0 max: 20 settlement_share: diff --git a/workflow/internal/settings.yaml b/workflow/internal/settings.yaml index 531ffd5..c48d532 100644 --- a/workflow/internal/settings.yaml +++ b/workflow/internal/settings.yaml @@ -12,7 +12,7 @@ resources: ghsl_tif: "GHS_BUILT_S_E2025_GLOBE_R2023A_4326_30ss_V1_0.tif" land_cover_types: - POST_FLOODING: FARM + POST_FLOODING_CROPLANDS: NOT_SUITABLE RAINFED_CROPLANDS: FARM MOSAIC_CROPLAND: FARM MOSAIC_VEGETATION: FARM @@ -23,13 +23,13 @@ land_cover_types: OPEN_NEEDLELEAVED_FOREST: FOREST CLOSED_TO_OPEN_MIXED_FOREST: FOREST MOSAIC_FOREST: FOREST - CLOSED_TO_OPEN_REGULARLY_FLOODED_FOREST: FOREST - CLOSED_REGULARLY_FLOODED_FOREST: FOREST + CLOSED_TO_OPEN_REGULARLY_FLOODED_FOREST: NOT_SUITABLE + CLOSED_REGULARLY_FLOODED_FOREST: NOT_SUITABLE MOSAIC_GRASSLAND: OTHER CLOSED_TO_OPEN_SHRUBLAND: OTHER CLOSED_TO_OPEN_HERBS: OTHER SPARSE_VEGETATION: OTHER - CLOSED_TO_OPEN_REGULARLY_FLOODED_GRASSLAND: OTHER + CLOSED_TO_OPEN_REGULARLY_FLOODED_GRASSLAND: NOT_SUITABLE BARE_AREAS: OTHER ARTIFICIAL_SURFACES_AND_URBAN_AREAS: URBAN WATER_BODIES: WATER diff --git a/workflow/rules/automatic.smk b/workflow/rules/automatic.smk index ea63978..0789271 100644 --- a/workflow/rules/automatic.smk +++ b/workflow/rules/automatic.smk @@ -8,7 +8,7 @@ if config.get("tiny_files", False): rule clip_slope: input: - vector="", + like_vector="", output: path="/automatic/cutout/{shape}/slope.tif", log: @@ -19,11 +19,11 @@ if config.get("tiny_files", False): message: "Download slope data covering the bounds of the input shapefile." wrapper: - "v7.2.0/geo/rasterio/clip-geotiff" + "v9.14.0/geo/rasterio/clip" rule clip_bathymetry: input: - vector="", + like_vector="", output: path="/automatic/cutout/{shape}/bathymetry.tif", log: @@ -34,7 +34,7 @@ if config.get("tiny_files", False): message: "Download bathymetry data covering the bounds of the input shapefile." wrapper: - "v7.2.0/geo/rasterio/clip-geotiff" + "v9.14.0/geo/rasterio/clip" else: @@ -56,7 +56,7 @@ else: "Download global slope data." shell: """ - curl -sSLo {output:q} {params.url:q} + curl -sSLo {output:q} {params.url:q} > {log:q} 2>&1 """ rule download_bathymetry: @@ -73,45 +73,34 @@ else: "Download global bathymetry data." shell: """ - curl -sSLo {output:q} {params.url:q} + curl -sSLo {output:q} {params.url:q} > {log:q} 2>&1 """ rule clip_slope: input: - script=workflow.source_path("../scripts/clip_raster.py"), - shapes="", - slope=rules.download_slope.output, + like_vector="", + raster=rules.download_slope.output[0], output: - "/automatic/cutout/{shape}/slope.tif", + path="/automatic/cutout/{shape}/slope.tif", log: "/{shape}/clip_slope.log", - conda: - "../envs/module.yaml" message: "Cut slope data to the bounds of the input shapefile." - shell: - """ - python {input.script:q} {input.slope:q} {input.shapes:q} {output:q} 2>{log:q} - """ + wrapper: + "v9.14.0/geo/rasterio/clip" rule clip_bathymetry: input: - script=workflow.source_path("../scripts/clip_raster.py"), - shapes="", - bathymetry=rules.download_bathymetry.output, + like_vector="", + raster=rules.download_bathymetry.output[0], output: - "/automatic/cutout/{shape}/bathymetry.tif", + path="/automatic/cutout/{shape}/bathymetry.tif", log: "/{shape}/clip_bathymetry.log", - conda: - "../envs/module.yaml" message: "Cut bathymetry data to the bounds of the input shapefile." - shell: - """ - python {input.script:q} {input.bathymetry:q} {input.shapes:q} {output:q} 2>{log:q} - """ - + wrapper: + "v9.14.0/geo/rasterio/clip" ## # Globcover @@ -132,7 +121,7 @@ rule download_globcover: "Download the GlobCover land cover data (~380 MB)." shell: """ - curl -sSLo {output:q} {params.url:q} + curl -sSLo {output:q} {params.url:q} > {log:q} 2>&1 """ @@ -152,27 +141,22 @@ rule unzip_globcover: "Unzip the relevant TIF files from the GlobCover zip file." shell: """ - python {input.script:q} {input.zipfile:q} -f {params.target_file:q} -o {output:q} 2>{log:q} + python {input.script:q} {input.zipfile:q} -f {params.target_file:q} -o {output:q} > {log:q} 2>&1 """ rule clip_landcover: input: - script=workflow.source_path("../scripts/clip_raster.py"), - shapes="", - landcover=rules.unzip_globcover.output, + like_vector="", + raster=rules.unzip_globcover.output[0], output: - "/automatic/cutout/{shape}/landcover.tif", + path="/automatic/cutout/{shape}/landcover.tif", log: "/{shape}/clip_landcover.log", - conda: - "../envs/module.yaml" message: "Cut land cover data to the bounds of the input shapefile." - shell: - """ - python {input.script:q} {input.landcover:q} {input.shapes:q} {output:q} 2>{log:q} - """ + wrapper: + "v9.14.0/geo/rasterio/clip" ## @@ -194,7 +178,7 @@ rule download_ghsl: "Download the GHSL (Global Human Settlement Layer) built-up surface data." shell: """ - curl -sSLo {output:q} {params.url:q} + curl -sSLo {output:q} {params.url:q} > {log:q} 2>&1 """ @@ -214,27 +198,22 @@ rule unzip_ghsl: "Unzip the relevant TIF file from the GHSL data." shell: """ - python {input.script:q} {input.zipfile:q} -f {params.target_file:q} -o {output:q} 2>{log:q} + python {input.script:q} {input.zipfile:q} -f {params.target_file:q} -o {output:q} > {log:q} 2>&1 """ rule clip_settlement: input: - script=workflow.source_path("../scripts/clip_raster.py"), - shapes="", - settlement=rules.unzip_ghsl.output, + like_vector="", + raster=rules.unzip_ghsl.output[0], output: - "/automatic/cutout/{shape}/settlement.tif", + path="/automatic/cutout/{shape}/settlement.tif", log: "/{shape}/clip_settlement.log", - conda: - "../envs/module.yaml" message: "Cut settlement data to the bounds of the input shapefile." - shell: - """ - python {input.script:q} {input.settlement:q} {input.shapes:q} {output:q} 2>{log:q} - """ + wrapper: + "v9.14.0/geo/rasterio/clip" ## @@ -246,7 +225,7 @@ rule rasterise_clip_wdpa: input: script=workflow.source_path("../scripts/clip_and_rasterise_polys.py"), shapes="", - reference_raster=rules.clip_landcover.output, + reference_raster=rules.clip_landcover.output[0], protected_areas="", output: "/automatic/cutout/{shape}/wdpa.tif", @@ -258,5 +237,5 @@ rule rasterise_clip_wdpa: "Rasterise and cut WDPA data to the bounds of the input shapefile, using the landcover raster as reference for the rasterisation." shell: """ - python {input.script:q} {input.shapes:q} {input.reference_raster:q} {input.protected_areas:q} {output:q} 2>{log:q} + python {input.script:q} {input.shapes:q} {input.reference_raster:q} {input.protected_areas:q} {output:q} > {log:q} 2>&1 """ diff --git a/workflow/rules/process.smk b/workflow/rules/process.smk index 7e1f7b1..9c7be82 100644 --- a/workflow/rules/process.smk +++ b/workflow/rules/process.smk @@ -14,7 +14,7 @@ checkpoint breakup_shape: "Break up {wildcards.shape} into the configured subunits." shell: """ - python {input.script:q} {input.shapes:q} {params.split_by:q} {output:q} 2>{log:q} + python {input.script:q} {input.shapes:q} {params.split_by:q} {output:q} > {log:q} 2>&1 """ @@ -49,7 +49,7 @@ 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} 2>{log:q} + {output.resampled_input:q} {output.plot:q} > {log:q} 2>&1 """ @@ -78,7 +78,7 @@ rule area_potential: "Compute area potential for the tech {wildcards.tech} and {wildcards.subunit} in {wildcards.shape}." shell: """ - python {input.script:q} "{input.shapes}/{wildcards.subunit}.parquet" {input.resampled_path:q} {params.config:q} {params.buffer_crs:q} {output.area_potential:q} {output.plot:q} --override_config={params.subunit_override_config:q} 2>{log:q} + python {input.script:q} "{input.shapes}/{wildcards.subunit}.parquet" {input.resampled_path:q} {params.config:q} {params.buffer_crs:q} {output.area_potential:q} {output.plot:q} --override_config={params.subunit_override_config:q} > {log:q} 2>&1 """ @@ -95,7 +95,7 @@ rule aggregate_area_potential: "Aggregate area potential for the tech {wildcards.tech} in {wildcards.shape}." shell: """ - gdalwarp --config GDAL_CACHEMAX 3000 -wm 3000 -of GTiff -co COMPRESS=LZW {input} {output.aggregated_area_potential:q} + gdalwarp --config GDAL_CACHEMAX 3000 -wm 3000 -of GTiff -co COMPRESS=LZW {input} {output.aggregated_area_potential:q} > {log:q} 2>&1 """ diff --git a/workflow/scripts/clip_raster.py b/workflow/scripts/clip_raster.py deleted file mode 100644 index 73b424c..0000000 --- a/workflow/scripts/clip_raster.py +++ /dev/null @@ -1,42 +0,0 @@ -"""Clip raster files based on the bounding box from a parquet shapefile.""" - -import subprocess - -import click - - -@click.command() -@click.argument("input_tif", type=click.Path(exists=True)) -@click.argument("input_parquet", type=click.Path(exists=True)) -@click.argument("output_tif", type=click.Path()) -def clip_raster(input_tif, input_parquet, output_tif): - """Clip INPUT_TIF using the bounding box from INPUT_PARQUET and save as OUTPUT_TIF. - - This script calls 'fio' and 'rio' directly, assuming they are installed. - - """ - try: - # Step 1: Get bounds from the input_parquet using fio - fio_cmd = ["fio", "info", input_parquet, "--bounds"] - result = subprocess.run(fio_cmd, capture_output=True, text=True, check=True) - bounds = result.stdout.strip() - - # Step 2: Run rio clip with the bounds obtained - rio_cmd = [ - "rio", - "clip", - "--overwrite", - input_tif, - output_tif, - "--bounds", - bounds, - ] - subprocess.run(rio_cmd, check=True) - - except subprocess.CalledProcessError as e: - click.echo(f"Error running command: {e.cmd}", err=True) - click.echo(e.stderr, err=True) - - -if __name__ == "__main__": - clip_raster() diff --git a/workflow/scripts/report.py b/workflow/scripts/report.py index 5deb5ff..f551b3c 100644 --- a/workflow/scripts/report.py +++ b/workflow/scripts/report.py @@ -1,5 +1,7 @@ """This script generates a report summarizing area potentials for different technologies.""" +import sys + import geopandas as gpd import pandas as pd import rioxarray as rxr @@ -58,6 +60,7 @@ def report(shapes, area_potentials, csv_path, html_path): if __name__ == "__main__": + sys.stderr = open(snakemake.log[0], "w", buffering=1) report( snakemake.input.shapes, snakemake.input.area_potentials, diff --git a/workflow/scripts/resample.py b/workflow/scripts/resample.py index 986778a..7f29cc9 100644 --- a/workflow/scripts/resample.py +++ b/workflow/scripts/resample.py @@ -18,7 +18,7 @@ GLOBCOVER_TYPES = { - 11: "POST_FLOODING", + 11: "POST_FLOODING_CROPLANDS", 14: "RAINFED_CROPLANDS", 20: "MOSAIC_CROPLAND", 30: "MOSAIC_VEGETATION", @@ -224,7 +224,7 @@ def resample_inputs( ## da_slope = rxr.open_rasterio(slope_path, masked=True) / 100 print(f"Slope resolution: {da_slope.rio.resolution()}") - resampled["slope"] = da_slope.rio.reproject_match( + resampled["slope_slope"] = da_slope.rio.reproject_match( reference_raster, resampling=Resampling.average ) del da_slope diff --git a/workflow/scripts/tif_to_png.py b/workflow/scripts/tif_to_png.py index 993be40..645ee9e 100644 --- a/workflow/scripts/tif_to_png.py +++ b/workflow/scripts/tif_to_png.py @@ -1,5 +1,7 @@ """This script plots a TIF file to PNG format.""" +import sys + import rioxarray as rxr from _script_utils import plot_all_dataset_variables @@ -15,4 +17,5 @@ def tif_to_png(tif_file_in, png_file_out): if __name__ == "__main__": + sys.stderr = open(snakemake.log[0], "w", buffering=1) tif_to_png(snakemake.input[0], snakemake.output[0]) From 502e44c8bdd66ae2cfd4aced9cc05579a06966ca Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 12 Aug 2026 19:53:48 +0000 Subject: [PATCH 2/5] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- workflow/rules/automatic.smk | 15 ++++++++------- workflow/rules/process.smk | 8 ++++---- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/workflow/rules/automatic.smk b/workflow/rules/automatic.smk index 0789271..0608660 100644 --- a/workflow/rules/automatic.smk +++ b/workflow/rules/automatic.smk @@ -56,7 +56,7 @@ else: "Download global slope data." shell: """ - curl -sSLo {output:q} {params.url:q} > {log:q} 2>&1 + curl -sSLo {output:q} {params.url:q} >{log:q} 2>&1 """ rule download_bathymetry: @@ -73,7 +73,7 @@ else: "Download global bathymetry data." shell: """ - curl -sSLo {output:q} {params.url:q} > {log:q} 2>&1 + curl -sSLo {output:q} {params.url:q} >{log:q} 2>&1 """ rule clip_slope: @@ -102,6 +102,7 @@ else: wrapper: "v9.14.0/geo/rasterio/clip" + ## # Globcover ## @@ -121,7 +122,7 @@ rule download_globcover: "Download the GlobCover land cover data (~380 MB)." shell: """ - curl -sSLo {output:q} {params.url:q} > {log:q} 2>&1 + curl -sSLo {output:q} {params.url:q} >{log:q} 2>&1 """ @@ -141,7 +142,7 @@ rule unzip_globcover: "Unzip the relevant TIF files from the GlobCover zip file." shell: """ - python {input.script:q} {input.zipfile:q} -f {params.target_file:q} -o {output:q} > {log:q} 2>&1 + python {input.script:q} {input.zipfile:q} -f {params.target_file:q} -o {output:q} >{log:q} 2>&1 """ @@ -178,7 +179,7 @@ rule download_ghsl: "Download the GHSL (Global Human Settlement Layer) built-up surface data." shell: """ - curl -sSLo {output:q} {params.url:q} > {log:q} 2>&1 + curl -sSLo {output:q} {params.url:q} >{log:q} 2>&1 """ @@ -198,7 +199,7 @@ rule unzip_ghsl: "Unzip the relevant TIF file from the GHSL data." shell: """ - python {input.script:q} {input.zipfile:q} -f {params.target_file:q} -o {output:q} > {log:q} 2>&1 + python {input.script:q} {input.zipfile:q} -f {params.target_file:q} -o {output:q} >{log:q} 2>&1 """ @@ -237,5 +238,5 @@ rule rasterise_clip_wdpa: "Rasterise and cut WDPA data to the bounds of the input shapefile, using the landcover raster as reference for the rasterisation." shell: """ - python {input.script:q} {input.shapes:q} {input.reference_raster:q} {input.protected_areas:q} {output:q} > {log:q} 2>&1 + python {input.script:q} {input.shapes:q} {input.reference_raster:q} {input.protected_areas:q} {output:q} >{log:q} 2>&1 """ diff --git a/workflow/rules/process.smk b/workflow/rules/process.smk index 9c7be82..805b340 100644 --- a/workflow/rules/process.smk +++ b/workflow/rules/process.smk @@ -14,7 +14,7 @@ checkpoint breakup_shape: "Break up {wildcards.shape} into the configured subunits." shell: """ - python {input.script:q} {input.shapes:q} {params.split_by:q} {output:q} > {log:q} 2>&1 + python {input.script:q} {input.shapes:q} {params.split_by:q} {output:q} >{log:q} 2>&1 """ @@ -49,7 +49,7 @@ 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} >{log:q} 2>&1 """ @@ -78,7 +78,7 @@ rule area_potential: "Compute area potential for the tech {wildcards.tech} and {wildcards.subunit} in {wildcards.shape}." shell: """ - python {input.script:q} "{input.shapes}/{wildcards.subunit}.parquet" {input.resampled_path:q} {params.config:q} {params.buffer_crs:q} {output.area_potential:q} {output.plot:q} --override_config={params.subunit_override_config:q} > {log:q} 2>&1 + python {input.script:q} "{input.shapes}/{wildcards.subunit}.parquet" {input.resampled_path:q} {params.config:q} {params.buffer_crs:q} {output.area_potential:q} {output.plot:q} --override_config={params.subunit_override_config:q} >{log:q} 2>&1 """ @@ -95,7 +95,7 @@ rule aggregate_area_potential: "Aggregate area potential for the tech {wildcards.tech} in {wildcards.shape}." shell: """ - gdalwarp --config GDAL_CACHEMAX 3000 -wm 3000 -of GTiff -co COMPRESS=LZW {input} {output.aggregated_area_potential:q} > {log:q} 2>&1 + gdalwarp --config GDAL_CACHEMAX 3000 -wm 3000 -of GTiff -co COMPRESS=LZW {input} {output.aggregated_area_potential:q} >{log:q} 2>&1 """ From 59c9c3a16b1fb43e170c89b7bd6aca70b1aec6be Mon Sep 17 00:00:00 2001 From: Ivan Ruiz Manuel <72193617+irm-codebase@users.noreply.github.com> Date: Thu, 13 Aug 2026 14:05:15 +0200 Subject: [PATCH 3/5] standardise processing to output EPSG:4326 --- INTERFACE.yaml | 5 +-- config/config.yaml | 2 +- workflow/Snakefile | 1 + workflow/internal/settings.yaml | 3 ++ workflow/rules/automatic.smk | 14 ++++---- workflow/rules/preprocess.smk | 37 ++++++++++++++++++++ workflow/rules/process.smk | 21 +---------- workflow/scripts/clip_and_rasterise_polys.py | 1 + workflow/scripts/normalise_shapes.py | 27 ++++++++++++++ 9 files changed, 81 insertions(+), 30 deletions(-) create mode 100644 workflow/rules/preprocess.smk create mode 100644 workflow/scripts/normalise_shapes.py diff --git a/INTERFACE.yaml b/INTERFACE.yaml index 20e10bf..1cfd1ca 100644 --- a/INTERFACE.yaml +++ b/INTERFACE.yaml @@ -16,7 +16,8 @@ pathvars: shapes: default: "/user/shapes/{shape}.parquet" description: > - Region geometries in parquet format. + Region geometries in GeoParquet format. + The module normalises these geometries to EPSG:4326 before processing. These should conform to the schema defined in https://github.com/calliope-project/module_geo_boundaries/blob/main/workflow/internal/shape.schema.yaml. wdpa: default: "/user/wdpa.gdb" @@ -25,7 +26,7 @@ pathvars: results: area_potential: default: "/{shape}/area_potential_{tech}.tif" - description: "Area potential GeoTIFF raster for the specified technology across all subunits." + description: "Area potential GeoTIFF raster in EPSG:4326 for the specified technology across all subunits." wildcards: shape: "Name of the shape to be processed, e.g., 'world', 'europe', 'MEX'." tech: "Name of the technology, e.g., 'pv_rooftop' or 'wind_offshore'. Available technologies are defined in the module configuration." diff --git a/config/config.yaml b/config/config.yaml index 9f5807a..d5af9a2 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -2,7 +2,7 @@ # - "UTM": project each shape to the UTM zone of its centroid for buffering # - "epsg:xxxx": use the specified CRS for all buffering # A good option is "epsg:8857" (WGS 84 / Equal Earth Greenwich) for global coverage -buffer_crs: "epsg:3035" +buffer_crs: "epsg:8857" split_by: country_id # likely country_id or shape_id diff --git a/workflow/Snakefile b/workflow/Snakefile index fcc85e9..82d2113 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -41,6 +41,7 @@ wildcard_constraints: # Add all your includes here. include: "rules/functions.smk" +include: "rules/preprocess.smk" include: "rules/automatic.smk" include: "rules/process.smk" diff --git a/workflow/internal/settings.yaml b/workflow/internal/settings.yaml index c48d532..57305df 100644 --- a/workflow/internal/settings.yaml +++ b/workflow/internal/settings.yaml @@ -11,6 +11,9 @@ resources: 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" +# NOTE: only works because GlobCover is already in this CRS! +target_crs: "EPSG:4326" + land_cover_types: POST_FLOODING_CROPLANDS: NOT_SUITABLE RAINFED_CROPLANDS: FARM diff --git a/workflow/rules/automatic.smk b/workflow/rules/automatic.smk index 0608660..3d7c184 100644 --- a/workflow/rules/automatic.smk +++ b/workflow/rules/automatic.smk @@ -8,7 +8,7 @@ if config.get("tiny_files", False): rule clip_slope: input: - like_vector="", + like_vector=rules.normalise_shapes.output.shapes, output: path="/automatic/cutout/{shape}/slope.tif", log: @@ -23,7 +23,7 @@ if config.get("tiny_files", False): rule clip_bathymetry: input: - like_vector="", + like_vector=rules.normalise_shapes.output.shapes, output: path="/automatic/cutout/{shape}/bathymetry.tif", log: @@ -78,7 +78,7 @@ else: rule clip_slope: input: - like_vector="", + like_vector=rules.normalise_shapes.output.shapes, raster=rules.download_slope.output[0], output: path="/automatic/cutout/{shape}/slope.tif", @@ -91,7 +91,7 @@ else: rule clip_bathymetry: input: - like_vector="", + like_vector=rules.normalise_shapes.output.shapes, raster=rules.download_bathymetry.output[0], output: path="/automatic/cutout/{shape}/bathymetry.tif", @@ -148,7 +148,7 @@ rule unzip_globcover: rule clip_landcover: input: - like_vector="", + like_vector=rules.normalise_shapes.output.shapes, raster=rules.unzip_globcover.output[0], output: path="/automatic/cutout/{shape}/landcover.tif", @@ -205,7 +205,7 @@ rule unzip_ghsl: rule clip_settlement: input: - like_vector="", + like_vector=rules.normalise_shapes.output.shapes, raster=rules.unzip_ghsl.output[0], output: path="/automatic/cutout/{shape}/settlement.tif", @@ -225,7 +225,7 @@ rule clip_settlement: rule rasterise_clip_wdpa: input: script=workflow.source_path("../scripts/clip_and_rasterise_polys.py"), - shapes="", + shapes=rules.normalise_shapes.output.shapes, reference_raster=rules.clip_landcover.output[0], protected_areas="", output: diff --git a/workflow/rules/preprocess.smk b/workflow/rules/preprocess.smk new file mode 100644 index 0000000..5e3c757 --- /dev/null +++ b/workflow/rules/preprocess.smk @@ -0,0 +1,37 @@ +"""Rules related to data harmonisation.""" + +rule normalise_shapes: + input: + shapes="", + output: + shapes="/automatic/normalized_shapes/{shape}.parquet", + log: + "/{shape}/normalise_shapes.log", + conda: + "../envs/module.yaml" + params: + crs=internal["target_crs"] + message: + "Validate and normalize {wildcards.shape} to {params.crs}." + script: + "../scripts/normalise_shapes.py" + + +checkpoint breakup_shape: + input: + script=workflow.source_path("../scripts/breakup_shape.py"), + shapes=rules.normalise_shapes.output.shapes, + output: + directory("/automatic/shapes/{shape}"), + log: + "/{shape}/breakup_shape.log", + conda: + "../envs/module.yaml" + params: + split_by=config["split_by"], + message: + "Break up {wildcards.shape} into the configured subunits." + shell: + """ + python {input.script:q} {input.shapes:q} {params.split_by:q} {output:q} >{log:q} 2>&1 + """ diff --git a/workflow/rules/process.smk b/workflow/rules/process.smk index 805b340..c5e3b38 100644 --- a/workflow/rules/process.smk +++ b/workflow/rules/process.smk @@ -1,22 +1,3 @@ -checkpoint breakup_shape: - input: - script=workflow.source_path("../scripts/breakup_shape.py"), - shapes="", - output: - directory("/automatic/shapes/{shape}"), - log: - "/{shape}/breakup_shape.log", - conda: - "../envs/module.yaml" - params: - split_by=config["split_by"], - message: - "Break up {wildcards.shape} into the configured subunits." - shell: - """ - python {input.script:q} {input.shapes:q} {params.split_by:q} {output:q} >{log:q} 2>&1 - """ - rule prepare_resampled_inputs: input: @@ -119,7 +100,7 @@ rule plot_aggregated_area_potential: rule area_potential_report: input: - shapes="", + shapes=rules.normalise_shapes.output.shapes, area_potentials=expand( workflow.pathvars.apply(""), tech=config["techs"].keys(), diff --git a/workflow/scripts/clip_and_rasterise_polys.py b/workflow/scripts/clip_and_rasterise_polys.py index 0da3bf4..448516e 100644 --- a/workflow/scripts/clip_and_rasterise_polys.py +++ b/workflow/scripts/clip_and_rasterise_polys.py @@ -21,6 +21,7 @@ def clip_and_rasterise_polys( xmin, ymin, xmax, ymax = shapes.total_bounds protected_areas = gpd.read_file(protected_area_path) print(f"Protected areas: {len(protected_areas)}") + protected_areas = protected_areas.to_crs(shapes.crs) protected_areas = protected_areas.cx[xmin:xmax, ymin:ymax] print(f"Protected areas after applying total_bounds: {len(protected_areas)}") diff --git a/workflow/scripts/normalise_shapes.py b/workflow/scripts/normalise_shapes.py new file mode 100644 index 0000000..a937681 --- /dev/null +++ b/workflow/scripts/normalise_shapes.py @@ -0,0 +1,27 @@ +"""Input shape normalisation.""" + +import sys +from typing import TYPE_CHECKING, Any + +import geopandas as gpd +from _schemas import ShapesSchema + +if TYPE_CHECKING: + snakemake: Any + + +def main() -> None: + """Standardise the provided shapes to an internally standardised CRS.""" + shapes = gpd.read_parquet(snakemake.input.shapes) + shapes = ShapesSchema.validate(shapes) + + target_crs = snakemake.params.crs + if not shapes.crs.equals(target_crs): + shapes = shapes.to_crs(target_crs) + + shapes.to_parquet(snakemake.output.shapes) + + +if __name__ == "__main__": + sys.stderr = open(snakemake.log[0], "w", buffering=1) + main() From 58d841dee8cf528a29132fed8dd4968eb128b580 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 13 Aug 2026 12:05:28 +0000 Subject: [PATCH 4/5] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- workflow/rules/preprocess.smk | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/workflow/rules/preprocess.smk b/workflow/rules/preprocess.smk index 5e3c757..0cd3099 100644 --- a/workflow/rules/preprocess.smk +++ b/workflow/rules/preprocess.smk @@ -1,5 +1,6 @@ """Rules related to data harmonisation.""" + rule normalise_shapes: input: shapes="", @@ -10,7 +11,7 @@ rule normalise_shapes: conda: "../envs/module.yaml" params: - crs=internal["target_crs"] + crs=internal["target_crs"], message: "Validate and normalize {wildcards.shape} to {params.crs}." script: From 8f69a3da7442022742aed2ba658566d017750c02 Mon Sep 17 00:00:00 2001 From: Ivan Ruiz Manuel <72193617+irm-codebase@users.noreply.github.com> Date: Thu, 13 Aug 2026 18:21:47 +0200 Subject: [PATCH 5/5] Run shape repair for special cases --- workflow/scripts/_schemas.py | 5 +++++ workflow/scripts/normalise_shapes.py | 1 + 2 files changed, 6 insertions(+) diff --git a/workflow/scripts/_schemas.py b/workflow/scripts/_schemas.py index ff17826..f40df5d 100644 --- a/workflow/scripts/_schemas.py +++ b/workflow/scripts/_schemas.py @@ -19,6 +19,11 @@ class ShapesSchema(pa.DataFrameModel): parent_name: Series[str] | None "Human-readable name in the parent dataset." + @pa.parser("geometry") + def make_geometries_valid(cls, geometries): + """Repair invalid geometries before validating them.""" + return geometries.make_valid().rename(geometries.name) + @pa.check("geometry", element_wise=True) def check_geometries(cls, geom): return (geom is not None) and (not geom.is_empty) and geom.is_valid diff --git a/workflow/scripts/normalise_shapes.py b/workflow/scripts/normalise_shapes.py index a937681..ab16f1d 100644 --- a/workflow/scripts/normalise_shapes.py +++ b/workflow/scripts/normalise_shapes.py @@ -19,6 +19,7 @@ def main() -> None: if not shapes.crs.equals(target_crs): shapes = shapes.to_crs(target_crs) + shapes = ShapesSchema.validate(shapes) shapes.to_parquet(snakemake.output.shapes)