Skip to content

CleanCoverage snapping snaps inside the polygons... leading to multipolygons #1476

Description

@theroggy

Context: I have a layer of polygons, and I would like to use a snapping distance of e.g. 50 cm in CleanCoverage so polygons would be snapped to a neighbour if there is only a free space between them smaller than this distance. Note that these are not gaps, as gaps needs to be a closed area between two polygons (or within a polygon), which is not the case in this example: the polygons don't touch.

"Problem/observation": By setting the snap distance to 50 cm, some geometries in the layer I'm running this on become multipolygons because apparently the snapping is also applied within the outer ring of the polygon itself. Hence, if there is a location in a polygon where the polygon becomes narrower than the snapping distance, the narrow part of the polygon "disappears" and both wider ends are kept, leading to the polygon becoming a multipolygon.

Question: Having narrower sections in the polygons than the snapping distance might be a bit odd, granted, but nonetheless... is this behaviour to be expected and/or wanted?

EDIT: a reproducer to show the problem. The reproducer is a "coverage" of a single polygon that has a narrow location in it... and running clean_coverage with a snap distance larger than the narrow location removes the narrow piece, resulting that the polygon becomes a multipolygon.

Green: the input polygon, Red hatched: the resulting multipolygon after clean_coverage:

Image
import tempfile
import shutil
from pathlib import Path

import geopandas as gpd
import shapely
from matplotlib import pyplot as plt
from osgeo import gdal

gdal.UseExceptions()

# Prepare test data
tmp_dir = Path(tempfile.gettempdir()) / "geos_test_case"
tmp_dir.mkdir(parents=True, exist_ok=True)
poly = shapely.from_wkt("POLYGON ((0 0, 15 0, 15 10, 10 1, 5 1, 5 10, 0 10, 0 0))")
input_gdf = gpd.GeoDataFrame(geometry=[poly], crs="EPSG:31370")
input_path = tmp_dir / "input.gpkg"
input_gdf.to_file(input_path)

# Apply cleaning coverage
output_path = tmp_dir / "output.gpkg"
gdal.alg.vector.clean_coverage(
    input=input_path, output=output_path, snapping_distance=1.5, overwrite=True
)

# Visualize the result
output_gdf = gpd.read_file(output_path)
ax = input_gdf.plot(facecolor="green", alpha=0.5)
output_gdf.plot(ax=ax, facecolor="none", edgecolor="red", hatch="///")
plt.show()

shutil.rmtree(tmp_dir)

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions