Shapely Self-Intersections – Identifying LineString Intersections in Python

intersectionlinepythonself-intersectionshapely

I am looking for an efficient way to identify the self-intersections in a shapely LineString.

I can check if there is an intersection by lineStringName.is_simple function. However, I couldn't find a way to identify:

  1. self-intersection points/ node number
  2. the number of intersections.

Are there any methods apart from making a list of lines inside the LineString and checking if any of them intersects with any of them?
Consider below as the LineString.

   import shapely

   shapely.wkt.loads('LINESTRING (9603380.577551289 2719693.31939431, 9602238.01822002 2719133.882441244, 9601011.900844947 2718804.012436028, 9599670.800095448 2718931.680117098, 9599567.204161201 2717889.384686942, 9600852.184025297 2721120.409265322, 9599710.80929024 2720511.270897166, 9602777.832940497 2718125.875545334)')

enter image description here

Best Answer

You can also polygonize the line to see if polygons are getting created because of self-intersection of line with itself. If returned polygons are equal to or more than 1, you can infer that line is self-intersecting and is forming closed polygons. Here is the code to do that


from shapely.ops import polygonize, unary_union
from shapely import wkt

ls = wkt.loads('LINESTRING (9603380.577551289 2719693.31939431, 9602238.01822002 2719133.882441244, 9601011.900844947 2718804.012436028, 9599670.800095448 2718931.680117098, 9599567.204161201 2717889.384686942, 9600852.184025297 2721120.409265322, 9599710.80929024 2720511.270897166, 9602777.832940497 2718125.875545334)')

polygons = list(polygonize(unary_union(ls)))
if len(polygons)>0:
    print("line is forming the polygons by intersecting itself")
    for p in polygons:
        print(p)

Here is the output of this code for your sample line

line is forming the polygons by intersecting itself
POLYGON ((9601676.146449726 2718982.718574256, 9601011.900844947 2718804.012436028, 9599970.383407902 2718903.160927319, 9600397.520253101 2719977.177480989, 9601676.146449726 2718982.718574256))
POLYGON ((9599970.383407902 2718903.160927319, 9599567.204161201 2717889.384686942, 9599670.800095448 2718931.680117098, 9599970.383407902 2718903.160927319))
POLYGON ((9600397.520253101 2719977.177480989, 9599710.80929024 2720511.270897166, 9600852.184025297 2721120.409265322, 9600397.520253101 2719977.177480989))

When you visualize these polygons they appear like this

enter image description here

Edit 1:

If you dont need the actual polygons you can skip the polygonize part and just do the unary_union of your input linestring. This will give you all non-intersecting polyline segments in your input polyline. You can check for duplicate coordinates in these segments to find the self-intersection points.

from shapely.ops import , unary_union
from shapely import wkt
import collections


ls = wkt.loads('LINESTRING (9603380.577551289 2719693.31939431, 9602238.01822002 2719133.882441244, 9601011.900844947 2718804.012436028, 9599670.800095448 2718931.680117098, 9599567.204161201 2717889.384686942, 9600852.184025297 2721120.409265322, 9599710.80929024 2720511.270897166, 9602777.832940497 2718125.875545334)')

mls = unary_union(ls)

coords_list = []
for non_itersecting_ls in mls:
    coords_list.extend(non_itersecting_ls.coords)

print([item for item, count in collections.Counter(coords_list).items() if count > 1])

It gives all the intersectio points as the output

[(9601676.146449726, 2718982.7185742557), (9599970.383407902, 2718903.160927319), (9600397.520253101, 2719977.1774809887)]

Here is how it looks when you visualize the intersection points.

enter image description here

Related Question