I've encountered a puzzling issue which is probably linked to some floating point prec

Edge of intersection not matching the edge of original polygon about geos HOT 2 CLOSED

martinfleis commented on August 25, 2024

Edge of intersection not matching the edge of original polygon

from geos.

Comments (2)

pramsey commented on August 25, 2024

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.

from geos.

theroggy commented on August 25, 2024

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()

from geos.

Edge of intersection not matching the edge of original polygon about geos HOT 2 CLOSED

Comments (2)

Related Issues (20)

Recommend Projects

React

Vue.js

Typescript

TensorFlow

Django

Laravel

D3

Recommend Topics

javascript

web

server

Machine learning

Visualization

Game

Recommend Org

Facebook

Microsoft

Google

Alibaba

D3

Tencent