Splitting self-intersecting polygon only returned one polygon in Shapely

Question

I am using Python 3.5 64 bit in Windows 7 64 bit, shapely version 1.5.13.

I have the following code that returned me a self-intersecting polygon:

import numpy as np
from shapely.geometry import Polygon, MultiPolygon
import matplotlib.pyplot as plt

x = np.array([ 0.38517325,  0.40859912,  0.43296919,  0.4583215 ,  0.4583215 ,
               0.43296919,  0.40859912,  0.38517325,  0.36265506,  0.34100929])
y = np.array([ 62.5       ,  56.17977528,  39.39698492,   0.        ,
               0.        ,  17.34605377,  39.13341671,  60.4180932 ,
               76.02574417,  85.47008547])
polygon = Polygon(np.c_[x, y])
plt.plot(*polygon.exterior.xy)

This is correct. Then I tried to obtain the two individual polygons by using buffer(0):

split_polygon = polygon.buffer(0)
plt.plot(*polygon.exterior.xy)
print(type(split_polygon))
plt.fill(*split_polygon.exterior.xy)

Unfortunately, it only returned of the the two polygons:

Could anyone please help? Thanks!

Slocum

Answer 1

The first step is to close the LineString to make a LinearRing, which is what Polygons are made of.

from shapely.geometry import LineString, MultiPolygon
from shapely.ops import polygonize, unary_union

# original data
ls = LineString(np.c_[x, y])
# closed, non-simple
lr = LineString(ls.coords[:] + ls.coords[0:1])
lr.is_simple  # False

However, note that it is non-simple, since the lines cross to make a bow-tie. (The widely used buffer(0) trick usually does not work for fixing bow-ties in my experience). This is unsuitable for a LinearRing, so it needs further work. Make it simple and MultiLineString with unary_union:

mls = unary_union(lr)
mls.geom_type  # MultiLineString'

Then use polygonize to find the Polygons from the linework:

for polygon in polygonize(mls):
    print(polygon)

Or if you want one MultiPolygon geometry:

mp = MultiPolygon(list(polygonize(mls)))

Volumetric

Thanks for this, seems to work in the majority of cases. I actually had to do a buffer(0) afterwards to get it to clean 100% of my features, but at least it didn't drop big chunks. Minor nitpick: no need to call list() when polygonizing, the MultiPolygon object will take a generator. – Addieaddiego

@Mike T: Thanks a lot! I am using it for a thesis at university and I am very happy I found that. Psi have also got to make some notes about the implementation and I still do not entirely get all your code, would you mind explaining it yet another time? – Toole

using shapely 1.8.0 you can use shapely.validation.make_valid(shapely.geometry.Polygon(np.c_[x, y])) to get the same Multipolygon in one line – Bamboozle

Answer 2

I struggled with this for a while still in 2020, and finally just wrote a method that cleans up self intersections.

This requires Shapely v 1.2.1 explain_validity() method to work.

def clean_bowtie_geom(base_linearring):
    base_polygon = Polygon(base_linearring)

    invalidity = explain_validity(base_polygon)
    invalid_regex = re.compile('^(Self-intersection)[[](.+)\s(.+)[]]$')
    match = invalid_regex.match(invalidity)
    if match:
        groups = match.groups()
        intersect_point = (float(groups[1]), float(groups[2]))
        new_linring_coords1 = []
        new_linring_coords2 = []
        pop_new_linring = False

        for i in range(0, len(base_linearring.coords)):
            if i == len(base_linearring.coords) - 1:
                end_point = base_linearring.coords[0]
            else:
                end_point = base_linearring.coords[i + 1]
            start_point = base_linearring.coords[i]

            if not pop_new_linring:
                if is_point_on_line_and_between(start=start_point, end=end_point, pt=intersect_point):
                    new_linring_coords2.append(intersect_point)
                    new_linring_coords1.append(intersect_point)
                    pop_new_linring = True
                else:
                    new_linring_coords1.append(start_point)

            else:
                new_linring_coords2.append(start_point)
                if is_point_on_line_and_between(start=start_point, end=end_point, pt=intersect_point):
                    new_linring_coords2.append(intersect_point)
                    pop_new_linring = False

        corrected_linear_ring1 = LinearRing(coordinates=new_linring_coords1)
        corrected_linear_ring2 = LinearRing(coordinates=new_linring_coords2)

        polygon1 = Polygon(corrected_linear_ring1)
        polygon2 = Polygon(corrected_linear_ring2)
        
def is_point_on_line_and_between(start, end, pt, tol=0.0005):
    """
    Checks to see if pt is directly in line and between start and end coords
    :param start: list or tuple of x, y coordinates of start point of line
    :param end: list or tuple of x, y coordinates of end point of line
    :param pt: list or tuple of x, y coordinates of point to check if it is on the line
    :param tol: Tolerance for checking if point on line
    :return: True if on the line, False if not on the line
    """
    v1 = (end[0] - start[0], end[1] - start[1])
    v2 = (pt[0] - start[0], pt[1] - start[1])
    cross = cross_product(v1, v2)
    if cross <= tol:
        # The point lays on the line, but need to check if in between
        if ((start[0] <= pt[0] <= end[0]) or (start[0] >= pt[0] >= end[0])) and ((start[1] <= pt[1] <= end[1]) or (start[1] >= pt[1] >= end[1])):
            return True
    return False

This is not the cleanest code, but it gets the job done for me.

Input is a LinearRing with self intersecting geometry (is_simple=False) and output can be either 2 LinearRings, or Two Polygons, whichever you prefer (or have condition to pick one or the other, the world is your oyster, really).

---

EDIT

---

In Shapely 1.8.0, new function added. shapely.validation.make_valid() will take a self intersecting Polygon and return a MultiPolygon with each polygon created by splitting at the self intersection point(s).

Disestablish

I think you meant dot product, not cross product? – Hammock

Answer 3

Found in comments that should be an answer (@RunOrVeith)

Starting from shapely 1.8.0 you should use shapely.validation.make_valid

from shapely.geometry import Polygon 
from shapely.validation import make_valid

# poly with multiple self-intersections (p.buffer doesn't work)
p = Polygon([[0,0], [0,20], [5,20], [-5,10], [5,0]])
print(p.is_valid) # False

valid_p = make_valid(p)
print(valid_p.is_valid) # True

Note: valid_p will be MultipPolygon with 3 Polygon in it (in current scenario)

Veroniqueverras