Skip to content

Instantly share code, notes, and snippets.

@alpha-beta-soup

alpha-beta-soup/simplify.py

Created October 10, 2022 20:45
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save alpha-beta-soup/9fa834293398b8e58721b913769e240a to your computer and use it in GitHub Desktop.
Save alpha-beta-soup/9fa834293398b8e58721b913769e240a to your computer and use it in GitHub Desktop.
Download ZIP
Geometry simplification script (GPKG, Snakemake, geopandas, pygeos)
Raw
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')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment