Geometry simplification script (GPKG, Snakemake, geopandas, pygeos)

simplify.py

import logging
from pathlib import Path
import sys
import warnings

smk = snakemake # type: ignore

logging.basicConfig(filename=Path(str(smk.log)), level=smk.params.get('logLevel', logging.INFO))
logging.getLogger().addHandler(logging.StreamHandler(sys.stdout)) # Log to stdout as well

# pyright: reportMissingModuleSource=false
# pyright: reportMissingImports=false
import geopandas as gpd
import pygeos

quantisation: int = smk.params.get('quantisation', 100) # units of CRS, e.g. a value of 100 → 100 m using NZTM (EPSG:2193)

logging.info('Opening input')
gdf = gpd.read_file(str(smk.input[0]), engine="pyogrio")

logging.info('Setting precision, and repairing any invalid geometries')
gdf.geometry = pygeos.make_valid(pygeos.simplify(gdf.geometry.values.data, tolerance=quantisation, preserve_topology=True))
# # An alternative to the above:
# gdf.geometry = pygeos.make_valid(pygeos.set_precision(gdf.geometry.values.data, quantisation))

logging.info('Removing records with empty and missing geometries')
logging.info(f'{gdf.size} records prior to filter')
warnings.filterwarnings('ignore', 'GeoSeries.isna', UserWarning)
gdf = gdf[~(gdf.geometry.is_empty | gdf.geometry.isna())]
logging.info(f'{gdf.size} records posterior to filter')

logging.info('Writing output')
gdf.to_file(smk.output[0], engine="pyogrio", driver='GPKG')