pythonshapely

Splitting polygon by line using intersection returns original polygon in Shapely

I reduced my problem to some MWE so that it can be reproduced. I'm using Shapely 2.0 I want to split a polygon P in 2 parts by a line L.

POLYGON ((292718.0381447676 6638193.414029885, 292694.50537013356 6637994.803334004, 292718.0381447676 6638193.414029885, 292718.9708331647 6638193.303518486, 292722.0992155936 6638192.97038651, 292725.48975053953 6638192.648911643, 292729.16966360115 6638192.34135049, 292733.16283984896 6638192.050457474, 292737.48817684915 6638191.779469689, 292742.1600771896 6638191.532042792, 292747.1894678763 6638191.312155023, 292752.58358723175 6638191.124185786, 292758.34681013395 6638190.973003162, 292753.9574101781 6637991.021175833, 292758.34681013395 6638190.973003162, 292753.9574101781 6637991.021175833, 292694.50537013356 6637994.803334004, 292718.0381447676 6638193.414029885))

First, I get the intersection between P and L. The result is a LineString L2 as expected : LINESTRING (292756.2221414001 6638094.18724655, 292743.0611 6638095.284, 292709.9881 6638095.284, 292706.4410368586 6638095.537357549)

Second, I split P by L2, and I expect it to return 2 polygons P1 and P2.

from shapely.ops import split

box = wkt.loads(my_polygon_as_str)
line = wt.loads(my_line_as_str)
result = split(box, line) # should return 2 polygons

Yet, the split() method returns only one polygon for some reason, which is the original P shape. POLYGON ((292694.50537013356 6637994.803334004, 292718.0381447676 6638193.414029885, 292718.9708331647 6638193.303518486, 292722.0992155936 6638192.97038651, 292725.48975053953 6638192.648911643, 292729.16966360115 6638192.34135049, 292733.16283984896 6638192.050457474, 292737.48817684915 6638191.779469689, 292742.1600771896 6638191.532042792, 292747.1894678763 6638191.312155023, 292752.58358723175 6638191.124185786, 292758.34681013395 6638190.973003162, 292753.9574101781 6637991.021175833, 292694.50537013356 6637994.803334004))

What is happening here ? Is it a matter of floating points with Shapely considering L2 boundaries not intersecting exactly on P sides ? If so, how can I achieve such operation with real world data ?

P.S: I want to achieve this for hundreds of polygons, and the splitting randomly works for a few of them (about 10%) which is really confusing.

Solution

  • As mentioned already in the comments, you should split with the entire, original line rather than first intersecting it with the polygon.

    The reason for this is that by intersecting first you can run into "robustness issues". When calculating an intersection point between two lines, the exact point where they cross will sometimes need more digits to be "exact" than the precision used by a computer. Hence the intersection point calculated between the line and the polygon will be slightly (up to some nanometers) off. Depending on where the point ends up exactly, this can lead to the line not crossing the polygon anymore and hence the split will fail.

    Reference: GEOS FAQ