Edge of intersection not matching the edge of original polygon

[martinfleis]

I've encountered a puzzling issue which is probably linked to some floating point precision (though setting a precision grid does not help), when the result of intersection between geometries does not have a boundary matching the expected boundary of the smallest of the intersecting polygons. You can see the original use case in pysal/momepy#598.

Here's a small example derived from the original issue using shapely.

import shapely

geom = shapely.from_wkt(
    "POLYGON ((1603268.502117987 6464060.781328565, 1603296.8217964454 6464047.851641227, 1603349.1085612718 6464035.338499875, 1603363.557831175 6464031.88480676, 1603361.0308787343 6464021.107210826, 1603317.7832565615 6463836.796863219, 1603217.2506244255 6463859.844110653, 1603202.3783404578 6463872.287568242, 1603157.794884393 6463914.581579376, 1603146.6963311615 6463924.630126579, 1603157.0490438067 6463936.205929175, 1603258.3275165292 6464049.413615812, 1603268.502117987 6464060.781328565))"
)
points = shapely.from_wkt(
    [
        "POINT (1603284.828798125 6464019.5353853395)",
        "POINT (1603323.0633351507 6464000.228348275)",
        "POINT (1603278.2522754562 6464030.13553264)",
        "POINT (1603304.896049631 6463958.075765142)",
        "POINT (1603264.5876622554 6463903.927670779)",
    ]
)
voronoi = shapely.get_parts(shapely.voronoi_polygons(shapely.union_all(points)))
clipped = shapely.intersection(voronoi, geom)
symmdiff = shapely.union_all(clipped).symmetric_difference(geom)
symmdiff.area
# 2.1684507522934382e-08

Visually:

geom:

output of voronoi:

clipped:

symmdiff:

My expectation is that the boundary in this case matches, hence the symmetric difference is an empty polygon. In the use case linked above, we do tessellation within many geometries in a coverage but due to this issue, the result is not a valid coverage and the spatial contiguity is broken.

Using GEOS 3.12.1 via shapely 2.0.3.

[pramsey]

No, you've introduced new vertices where the voronoi edges hit the boundary edges, so the old and new boundaries won't be the same. New vertices have to fit into the precision grid of the underlying numeric format (IEEE doubles in this case), they do not fall into the precise "rational number line" location you are thinking of.

[theroggy]

As I started typing this already yesterday, I can as well post it as an elaboration on the answer by @pramsey.

During the clip operation, new points will need to be calculated at the crossings between the voronoi boundaries and the original geom, as marked with red circles in the image below. Due to non-infinite precision, in many cases this point will not be "exactly" on the crossing... so if you overlay with the original polygon in this case you'll get gaps/slivers. In my experience these slivers are always narrower than 1e-8, and this is also the case in this example.

2 options:

• if you run set_precision on the result of the overlay those slivers will be cleaned up. This is what I always do in practice.
• make the data topologically sound, so make sure that such calculated points are also added to neighbouring polygons. If the points are introduced in the neighbours, the edges fit perfectly and there will no gaps or slivers. Hopefully GEOSCoverageMakeValid is on the way to support this ;-)...

Script used

from matplotlib import pyplot as plt
import shapely
import shapely.plotting as plot

geom = shapely.from_wkt(
    "POLYGON ((1603268.502117987 6464060.781328565, 1603296.8217964454 6464047.851641227, 1603349.1085612718 6464035.338499875, 1603363.557831175 6464031.88480676, 1603361.0308787343 6464021.107210826, 1603317.7832565615 6463836.796863219, 1603217.2506244255 6463859.844110653, 1603202.3783404578 6463872.287568242, 1603157.794884393 6463914.581579376, 1603146.6963311615 6463924.630126579, 1603157.0490438067 6463936.205929175, 1603258.3275165292 6464049.413615812, 1603268.502117987 6464060.781328565))"
)
points = shapely.from_wkt(
    [
        "POINT (1603284.828798125 6464019.5353853395)",
        "POINT (1603323.0633351507 6464000.228348275)",
        "POINT (1603278.2522754562 6464030.13553264)",
        "POINT (1603304.896049631 6463958.075765142)",
        "POINT (1603264.5876622554 6463903.927670779)",
    ]
)
voronoi = shapely.get_parts(shapely.voronoi_polygons(shapely.union_all(points)))
clipped = shapely.intersection(voronoi, geom)
symmdiff = shapely.union_all(clipped).symmetric_difference(geom)
print(f"{symmdiff=}")
# output: symmdiff=<MULTIPOLYGON (...)>
print(f"{shapely.set_precision(symmdiff, 1e-8)=}")
# output: shapely.set_precision(symmdiff, 1e-8)=<MULTIPOLYGON EMPTY>

figure, ax = plt.subplots()
for voronoi_geom in voronoi:
    plot.plot_polygon(voronoi_geom, ax=ax, color="blue")
plot.plot_polygon(geom, ax=ax, color="green")
plt.show()