diff --git a/workflow/internal/config.schema.yaml b/workflow/internal/config.schema.yaml index 8dfa47f..b8edb2a 100644 --- a/workflow/internal/config.schema.yaml +++ b/workflow/internal/config.schema.yaml @@ -88,6 +88,19 @@ $defs: default: 1 minimum: 0 + gdalSettings: + description: Optional GDAL runtime configuration. + type: object + additionalProperties: false + properties: + geojson_max_obj_size_mb: + description: > + Maximum accepted size in megabytes for a single GeoJSON object read by + GDAL. Set to 0 to disable GDAL's single-object limit. + type: integer + default: 1024 + minimum: 0 + scenario: type: object required: @@ -222,6 +235,8 @@ properties: $ref: "#/$defs/crsSettings" voronoi_eez: $ref: "#/$defs/partialVoronoiEEZSettings" + gdal: + $ref: "#/$defs/gdalSettings" overture_release: description: | Overture data release to use. Two options are possible: diff --git a/workflow/internal/settings.yaml b/workflow/internal/settings.yaml index 2315670..b56599c 100644 --- a/workflow/internal/settings.yaml +++ b/workflow/internal/settings.yaml @@ -12,4 +12,6 @@ timeouts: connect_seconds: 15 read_seconds: 30 initial_retry_seconds: 5 +gdal: + geojson_max_obj_size_mb: 1024 overture_release: "latest" diff --git a/workflow/rules/_utils.smk b/workflow/rules/_utils.smk index 6db3347..46c6036 100644 --- a/workflow/rules/_utils.smk +++ b/workflow/rules/_utils.smk @@ -35,6 +35,12 @@ def get_voronoi_eez_config(scenario: str) -> dict: return internal["voronoi_eez"] | user_overrides +def get_gdal_config() -> dict: + """Get GDAL configuration with user overrides.""" + gdal = internal["gdal"] | config.get("gdal", {}) + return gdal + + def get_eez_file(scenario: str, country: str) -> str: """Build an EEZ filename reusable by scenarios with matching EEZ requests.""" extra_eez = config["scenarios"][scenario]["countries"][country].get("extra_eez", []) diff --git a/workflow/rules/download.smk b/workflow/rules/download.smk index 6d28fd8..b4b3a54 100644 --- a/workflow/rules/download.smk +++ b/workflow/rules/download.smk @@ -26,6 +26,7 @@ rule download_geoboundaries: "../envs/shape.yaml" params: timeouts=internal["timeouts"], + geojson_max_obj_size_mb=get_gdal_config()["geojson_max_obj_size_mb"], message: "Downloading '{wildcards.country}_{wildcards.subtype}_{wildcards.release_type}' dataset from geoBoundaries." script: @@ -42,6 +43,7 @@ rule download_gadm: "../envs/shape.yaml" params: timeouts=internal["timeouts"], + geojson_max_obj_size_mb=get_gdal_config()["geojson_max_obj_size_mb"], message: "Download '{wildcards.country}_{wildcards.subtype}' dataset from GADM." script: diff --git a/workflow/scripts/_geo.py b/workflow/scripts/_geo.py index 0e48a82..054c789 100644 --- a/workflow/scripts/_geo.py +++ b/workflow/scripts/_geo.py @@ -8,7 +8,7 @@ import geopandas as gpd from pyproj import CRS from rasterio import warp -from shapely import MultiPolygon, Polygon +from shapely import MultiPolygon, Polygon, make_valid from shapely.geometry import GeometryCollection, mapping, shape from shapely.geometry.base import BaseGeometry @@ -48,6 +48,36 @@ def extract_polygonal_geometry( return result +def make_geometry_valid(geom: BaseGeometry | None) -> Polygon | MultiPolygon | None: + """Return a valid polygonal geometry, dropping empty or non-polygonal parts.""" + result = None + + if geom is not None and not geom.is_empty: + if isinstance(geom, (Polygon, MultiPolygon)) and geom.is_valid: + result = geom + else: + polygonal = extract_polygonal_geometry(make_valid(geom)) + if ( + polygonal is not None + and not polygonal.is_empty + and not polygonal.is_valid + ): + polygonal = extract_polygonal_geometry(make_valid(polygonal)) + + if polygonal is not None and not polygonal.is_empty and polygonal.is_valid: + result = polygonal + + return result + + +def make_geometries_valid(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame: + """Apply polygonal geometry repair to a GeoDataFrame and drop empty results.""" + result = gdf.copy() + result["geometry"] = result["geometry"].map(make_geometry_valid) + result = result.loc[result["geometry"].notna()] + return result + + def _rasterio_to_crs(gdf: gpd.GeoDataFrame, to_crs: CRS) -> gpd.GeoDataFrame: """CRS conversion using rasterio's more powerful toolset. diff --git a/workflow/scripts/_schemas.py b/workflow/scripts/_schemas.py index 1e878c3..3a066f1 100644 --- a/workflow/scripts/_schemas.py +++ b/workflow/scripts/_schemas.py @@ -1,7 +1,7 @@ """Reusable schemas.""" import geopandas as gpd -from _geo import extract_polygonal_geometry +from _geo import make_geometries_valid from pandera import pandas as pa from pandera.typing.geopandas import GeoSeries from pandera.typing.pandas import Series @@ -37,14 +37,7 @@ class Config: @pa.dataframe_parser def fix_geometries(cls, gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame: # type: ignore[misc] """Attempt to correct empty, malformed, or non-polygonal geometries.""" - mask = gdf["geometry"].notna() & ~gdf["geometry"].is_empty - gdf = gdf.loc[mask].copy() - - invalid = ~gdf.geometry.is_valid - gdf.loc[invalid, "geometry"] = gdf.loc[invalid, "geometry"].make_valid() - - gdf["geometry"] = gdf["geometry"].apply(extract_polygonal_geometry) - return gdf.loc[gdf["geometry"].notna()] + return make_geometries_valid(gdf) @pa.check("geometry", element_wise=True) def check_geometries(cls, geom) -> bool: diff --git a/workflow/scripts/_utils.py b/workflow/scripts/_utils.py index 5cbf4f4..54d3905 100644 --- a/workflow/scripts/_utils.py +++ b/workflow/scripts/_utils.py @@ -5,6 +5,7 @@ from pathlib import Path import geopandas as gpd +import pyogrio import requests from matplotlib import pyplot as plt from matplotlib.axes import Axes @@ -75,3 +76,11 @@ def download_file(url: str, path: Path, timeouts: DownloadTimeouts) -> None: raise RuntimeError( f"Maximum retries {max_retries} reached without a valid response" ) + + +def read_geojson_file(path: Path, geojson_max_obj_size_mb: int) -> gpd.GeoDataFrame: + """Read a GeoJSON file with the configured GDAL object-size limit.""" + pyogrio.set_gdal_config_options( + {"OGR_GEOJSON_MAX_OBJ_SIZE": str(geojson_max_obj_size_mb)} + ) + return gpd.read_file(path) diff --git a/workflow/scripts/build_combined_area.py b/workflow/scripts/build_combined_area.py index 6a99b9c..ca533bd 100644 --- a/workflow/scripts/build_combined_area.py +++ b/workflow/scripts/build_combined_area.py @@ -23,6 +23,7 @@ def remove_overlaps(gdf: gpd.GeoDataFrame, crs: dict[str, CRS]) -> gpd.GeoDataFr from the later geometry. """ projected = _geo.to_projected_crs(gdf, crs["projected"]) + projected = _geo.make_geometries_valid(projected) left_idx, right_idx = projected.sindex.query( projected.geometry, predicate="intersects" @@ -51,13 +52,10 @@ def remove_overlaps(gdf: gpd.GeoDataFrame, crs: dict[str, CRS]) -> gpd.GeoDataFr continue geom = geom.difference(shapely.unary_union(cutters)) - if geom is not None and not geom.is_empty and not geom.is_valid: - geom = shapely.make_valid(geom) - - geoms[right] = geom + geoms[right] = _geo.make_geometry_valid(geom) projected["geometry"] = geoms - projected = projected.loc[projected.geometry.notna() & ~projected.geometry.is_empty] + projected = projected.loc[projected.geometry.notna()] result = _geo.to_geographic_crs(projected, crs["geographic"]) return result diff --git a/workflow/scripts/build_country.py b/workflow/scripts/build_country.py index 25a7dd0..add471b 100644 --- a/workflow/scripts/build_country.py +++ b/workflow/scripts/build_country.py @@ -12,7 +12,7 @@ import geopandas as gpd import pandas as pd from pyproj import CRS -from shapely import make_valid, voronoi_polygons +from shapely import voronoi_polygons from shapely.geometry import ( GeometryCollection, LineString, @@ -233,7 +233,7 @@ def _split_one_maritime( cells = gpd.GeoDataFrame(geometry=[], crs=crs) else: # Run Voronoi - maritime_geometry = make_valid(maritime_row.geometry) + maritime_geometry = maritime_row.geometry cells = gpd.GeoDataFrame( geometry=list( voronoi_polygons( @@ -243,7 +243,7 @@ def _split_one_maritime( ), crs=crs, ) - cells["geometry"] = cells.geometry.map(make_valid) + cells = _geo.make_geometries_valid(cells) # Identify where each voronoi cell belongs cells = gpd.sjoin( @@ -259,15 +259,18 @@ def _split_one_maritime( # Keep only cells on the maritime area and test for completeness cells["geometry"] = cells.geometry.intersection(maritime_geometry) - cells = cells.loc[~cells.geometry.is_empty & (cells.geometry.area > 0)] + cells = _geo.make_geometries_valid(cells) + cells = cells.loc[cells.geometry.area > 0] pieces = ( cells[["assigned_shape_id", "geometry"]] .dissolve(by="assigned_shape_id", as_index=False) .rename(columns={"assigned_shape_id": "_assigned_shape_id"}) ) - # FIXME: check if rerunning `warp` is necessary to fix things near the antimeridian - pieces = _geo.to_projected_crs(pieces, crs) - uncovered_area = maritime_geometry.difference(pieces.geometry.union_all()).area + pieces = _geo.make_geometries_valid(pieces) + uncovered = _geo.make_geometry_valid( + maritime_geometry.difference(pieces.geometry.union_all()) + ) + uncovered_area = 0 if uncovered is None else uncovered.area if uncovered_area > voronoi_config.uncovered_area_tolerance: raise RuntimeError( f"Voronoi split did not fully cover maritime shape {maritime_row.shape_id!r}. " @@ -297,6 +300,7 @@ def split_maritime_by_shoreline_voronoi( return shapes, cells shapes_proj = _geo.to_projected_crs(shapes, crs["projected"]) + shapes_proj = _geo.make_geometries_valid(shapes_proj) land = shapes_proj.loc[shapes_proj["shape_class"].eq("land")].copy() maritime = shapes_proj.loc[shapes_proj["shape_class"].eq("maritime")].copy() if maritime.empty: @@ -326,7 +330,6 @@ def split_maritime_by_shoreline_voronoi( geometry="geometry", crs=crs["projected"], ) - result.geometry = result.geometry.make_valid() result = _geo.to_geographic_crs(result, crs["geographic"]) return result, cells @@ -341,14 +344,16 @@ def build_country( # Reproject p_land = _geo.to_projected_crs(land, crs["projected"]) p_marine = _geo.to_projected_crs(maritime, crs["projected"]) + p_land = _geo.make_geometries_valid(p_land) + p_marine = _geo.make_geometries_valid(p_marine) # Remove contested zones p_marine = p_marine[p_marine["contested"].eq(False)].drop(columns="contested") # Give priority to EEZs p_land.geometry = p_land.geometry.difference(p_marine.geometry.union_all()) + p_land = _geo.make_geometries_valid(p_land) combined = gpd.GeoDataFrame( pd.concat([p_land, p_marine], ignore_index=True), crs=crs["projected"] ) - combined.geometry = combined.geometry.make_valid() return _geo.to_geographic_crs(combined, crs["geographic"]) diff --git a/workflow/scripts/download_gadm.py b/workflow/scripts/download_gadm.py index e437a64..5ff2f94 100644 --- a/workflow/scripts/download_gadm.py +++ b/workflow/scripts/download_gadm.py @@ -10,7 +10,7 @@ from typing import TYPE_CHECKING, Any import geopandas as gpd -from _utils import DownloadTimeouts, download_file +from _utils import DownloadTimeouts, download_file, read_geojson_file if TYPE_CHECKING: snakemake: Any @@ -23,7 +23,7 @@ def download_country_gadm( - country: str, subtype: str, timeouts: DownloadTimeouts + country: str, subtype: str, timeouts: DownloadTimeouts, geojson_max_obj_size_mb: int ) -> gpd.GeoDataFrame: """Attempts to download country GADM data in .json or zipped json.""" last_error: Exception | None = None @@ -35,7 +35,7 @@ def download_country_gadm( tmp_path = Path(tmp_dir) / f"download.json{zip_ext}" download_file(url, tmp_path, timeouts) - gdf = gpd.read_file(tmp_path) + gdf = read_geojson_file(tmp_path, geojson_max_obj_size_mb) if gdf.empty: raise RuntimeError(f"Downloaded empty GADM file from {url!r}.") return gdf.to_crs(GADM_CRS) @@ -51,7 +51,10 @@ def main(): """Main snakemake process.""" timeouts = DownloadTimeouts(**snakemake.params.timeouts) country = download_country_gadm( - snakemake.wildcards.country, snakemake.wildcards.subtype, timeouts + snakemake.wildcards.country, + snakemake.wildcards.subtype, + timeouts, + snakemake.params.geojson_max_obj_size_mb, ) country.to_parquet(snakemake.output.path) diff --git a/workflow/scripts/download_geoboundaries.py b/workflow/scripts/download_geoboundaries.py index 526c80e..68ab2b7 100644 --- a/workflow/scripts/download_geoboundaries.py +++ b/workflow/scripts/download_geoboundaries.py @@ -10,7 +10,7 @@ from typing import TYPE_CHECKING, Any import geopandas as gpd -from _utils import DownloadTimeouts, download_file +from _utils import DownloadTimeouts, download_file, read_geojson_file if TYPE_CHECKING: snakemake: Any @@ -23,7 +23,11 @@ def download_country_geoboundaries( - country: str, subtype: str, release_type: str, timeouts: DownloadTimeouts + country: str, + subtype: str, + release_type: str, + timeouts: DownloadTimeouts, + geojson_max_obj_size_mb: int, ) -> gpd.GeoDataFrame: """Download country data from geoBoundaries. @@ -46,7 +50,7 @@ def download_country_geoboundaries( geojson_path = tmp_path / "download.geojson" download_file(geojson_url, geojson_path, timeouts) - gdf = gpd.read_file(geojson_path) + gdf = read_geojson_file(geojson_path, geojson_max_obj_size_mb) if gdf.empty: raise RuntimeError( f"Downloaded empty geoBoundaries file from {geojson_url!r}." @@ -64,6 +68,7 @@ def main(): snakemake.wildcards.subtype, snakemake.wildcards.release_type, timeouts, + snakemake.params.geojson_max_obj_size_mb, ) country.to_parquet(snakemake.output.path)