Identifying self-intersections in LineString using shapely

Question

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

Answer 1

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

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.

Answer 2

The problem with the @Taras solution (.explain_validity()) is that only the first point is given

Another solution is to split the LineString in segments and find the intersections between these segments

l = 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)')
segments = list(map(LineString, zip(l.coords[:-1], l.coords[1:])))
for line in segments:
    print(line)  
LINESTRING (9603380.577551289 2719693.31939431, 9602238.01822002 2719133.882441244)
LINESTRING (9602238.01822002 2719133.882441244, 9601011.900844947 2718804.012436028)
LINESTRING (9601011.900844947 2718804.012436028, 9599670.800095448 2718931.680117098)
LINESTRING (9599670.800095448 2718931.680117098, 9599567.204161201 2717889.384686942)
LINESTRING (9599567.204161201 2717889.384686942, 9600852.184025297 2721120.409265322)
LINESTRING (9600852.184025297 2721120.409265322, 9599710.80929024 2720511.270897166)
LINESTRING (9599710.80929024 2720511.270897166, 9602777.832940497 2718125.875545334)

Now use itertools to find the intersections

for seg1,seg2 in itertools.combinations(segments,2):
   if seg1.crosses(seg2):
        inter = seg1.intersection(seg2)        
        print(inter) 
 POINT (9601676.146449726 2718982.718574256)
 POINT (9599970.383407902 2718903.160927319)
 POINT (9600397.520253101 2719977.177480989)

Answer 3

As you already mentioned, checking whether feature collides with itself can be done by means of the .is_simple

Returns True if the feature does not cross itself.

from shapely.wkt import loads

l = 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)')

print(l.is_simple) # False

To a certain extent the LinearRing constructor can be helpful

The LinearRing constructor takes an ordered sequence of (x, y[, z]) point tuples.

together with .explain_validity()

Returns a string explaining the validity or invalidity of the object.

from shapely.wkt import loads
from shapely.validation import explain_validity
from shapely.geometry.polygon import LinearRing

l = 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)')

ring = LinearRing(list(l.coords))

print(ring.is_simple) # False
print(explain_validity(ring)) # Ring Self-intersection[9600397.5202531 2719977.17748099]

One of the points: