7

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

0

3 Answers 3

6

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

3
  • Hello @Abhi, I appreciate your response. It is a possible option. But I am actually looking for an efficient solution. In this case, the process will get more intense. Thank you
    – rmj
    Commented Feb 10, 2022 at 9:29
  • You can also skip the polygonize part and just used the unary_union part, it will give you a list of non-intersecting polyline segments bundled together in a MultiLineString. I will update the answer with this info.
    – Abhilshit
    Commented Feb 10, 2022 at 13:34
  • 1
    Thank you @Abhilshit it looks good now.
    – rmj
    Commented Feb 11, 2022 at 15:39
4

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)

enter image description here

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)

enter image description here

1
  • thanks for the response. I have already specified that I am not looking for a solution that breaks the LineString and checks if any of them is crossing any of them as the itertool does. I believe there is no direct solution for this issue. However, itertools could speed up the process and make it efficient which could be better than the for loops I mentioned. Hence, I believe that your solution is the best possible one. Thank you so much
    – rmj
    Commented Feb 10, 2022 at 9:24
4

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:

result

1
  • Thank you for the response @Taras, but it won't give the entire list/ total number of intersections (might be with a for loop though).
    – rmj
    Commented Feb 10, 2022 at 9:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.