0

I ran into a problem where shapely.intersection cuts a LineString into a MultiLinestring with lots of short lines instead of a simple LineString. The intersection is made between a circle (Polygon) and a LineString. The python code is simple:

intersection = shapely.intersection(line, circle)

The original LineString contains duplicate points, since it is a GPS track going back and forth the same way: Visualization The black circle intersects with the blue line (GPS track). The red line is the result.

These are the sub Lines of the intersection separated a bit to visualize it better: enter image description here

Here I separated the points of the original GPS track so it is easier to see that it goes back the same way: enter image description here

This is the intersection:

    {
  "type": "Feature",
  "geometry": {
    "type": "MultiLineString",
    "coordinates": [
      [
        [
          8.987437,
          60.710481
        ],
        [
          8.987753,
          60.710386
        ]
      ],
      [
        [
          8.987753,
          60.710386
        ],
        [
          8.987975,
          60.710381
        ],
        [
          8.98801,
          60.710381
        ]
      ],
      [
        [
          8.98801,
          60.710381
        ],
        [
          8.988991,
          60.709666
        ]
      ],
      [
        [
          8.988991,
          60.709666
        ],
        [
          8.989304,
          60.709225
        ]
      ],
      [
        [
          8.98801,
          60.710381
        ],
        [
          8.987753,
          60.710386
        ]
      ]
    ]
  },
  "properties": {
    "stroke": "#ff0000",
    "stroke-width": 2,
    "stroke-opacity": 1
  },
  "id": 2
}

As you can see the last point of a sub LineString and the first point of the next sub LineString are the same. These were duplicate points of the original LineString. The duplicate points belong there because the track goes back the same way as it came. Here is the full geojson file to visualize in geojson.io:

{ "type": "FeatureCollection", "features": [ { "type": "Feature", "geometry": { "type": "LineString", "coordinates": [ [ 8.985738, 60.712983 ], [ 8.985523, 60.712769 ], [ 8.985278, 60.712175 ], [ 8.985441, 60.711763 ], [ 8.985674, 60.711554 ], [ 8.986364, 60.711158 ], [ 8.986783, 60.71076 ], [ 8.986913, 60.710638 ], [ 8.987753, 60.710386 ], [ 8.987975, 60.710381 ], [ 8.98801, 60.710381 ], [ 8.988991, 60.709666 ], [ 8.989304, 60.709225 ], [ 8.988991, 60.709666 ], [ 8.98801, 60.710381 ], [ 8.987753, 60.710386 ], [ 8.986913, 60.710638 ], [ 8.986364, 60.711158 ], [ 8.985674, 60.711554 ], [ 8.985427, 60.711546 ], [ 8.985205, 60.711582 ], [ 8.985061, 60.711618 ], [ 8.984823, 60.711663 ], [ 8.984621, 60.711681 ], [ 8.984332, 60.711681 ], [ 8.98412, 60.71166 ], [ 8.983526, 60.711579 ], [ 8.982879, 60.711582 ] ] }, "properties": { "stroke": "#0000a0", "stroke-width": 2, "stroke-opacity": 1 }, "id": 0 }, { "type": "Feature", "geometry": { "type": "Polygon", "coordinates": [ [ [ 8.993146, 60.708983 ], [ 8.993075, 60.708632 ], [ 8.992866, 60.708294 ], [ 8.992526, 60.707983 ], [ 8.992069, 60.707711 ], [ 8.991512, 60.707487 ], [ 8.990876, 60.707321 ], [ 8.990187, 60.707219 ], [ 8.989469, 60.707184 ], [ 8.988752, 60.707219 ], [ 8.988063, 60.707321 ], [ 8.987427, 60.707487 ], [ 8.98687, 60.707711 ], [ 8.986413, 60.707983 ], [ 8.986073, 60.708294 ], [ 8.985864, 60.708632 ], [ 8.985793, 60.708983 ], [ 8.985864, 60.709334 ], [ 8.986073, 60.709671 ], [ 8.986413, 60.709982 ], [ 8.98687, 60.710255 ], [ 8.987427, 60.710478 ], [ 8.988062, 60.710644 ], [ 8.988752, 60.710747 ], [ 8.989469, 60.710781 ], [ 8.990187, 60.710747 ], [ 8.990876, 60.710644 ], [ 8.991512, 60.710478 ], [ 8.992069, 60.710255 ], [ 8.992526, 60.709982 ], [ 8.992866, 60.709671 ], [ 8.993075, 60.709334 ], [ 8.993146, 60.708983 ] ] ] }, "properties": {}, "id": 1 }, { "type": "Feature", "geometry": { "type": "MultiLineString", "coordinates": [ [ [ 8.987437, 60.710481 ], [ 8.987753, 60.710386 ] ], [ [ 8.987753, 60.710386 ], [ 8.987975, 60.710381 ], [ 8.98801, 60.710381 ] ], [ [ 8.98801, 60.710381 ], [ 8.988991, 60.709666 ] ], [ [ 8.988991, 60.709666 ], [ 8.989304, 60.709225 ] ], [ [ 8.98801, 60.710381 ], [ 8.987753, 60.710386 ] ] ] }, "properties": { "stroke": "#ff0000", "stroke-width": 2, "stroke-opacity": 1 }, "id": 2 } ] }

I couldn't figure out why this happens and how I can avoid it. Simply merging the MultiLineString to a LineString is no option, because it can contain legitimate MultiLineString segments if the track leaves the circle, but reenters later on. Then multiple LineStrings inside of a MultiLineString are desired. But not the MultiLineString due to the duplicate points.

This should be a reproducable example in Python 3.7, Shapely 2.0.2:

import shapely

line = shapely.linestrings([ [ 8.985738, 60.712983 ], [ 8.985523, 60.712769 ], [ 8.985278, 60.712175 ], [ 8.985441, 60.711763 ], [ 8.985674, 60.711554 ], [ 8.986364, 60.711158 ], [ 8.986783, 60.71076 ], [ 8.986913, 60.710638 ], [ 8.987753, 60.710386 ], [ 8.987975, 60.710381 ], [ 8.98801, 60.710381 ], [ 8.988991, 60.709666 ], [ 8.989304, 60.709225 ], [ 8.988991, 60.709666 ], [ 8.98801, 60.710381 ], [ 8.987753, 60.710386 ], [ 8.986913, 60.710638 ], [ 8.986364, 60.711158 ], [ 8.985674, 60.711554 ], [ 8.985427, 60.711546 ], [ 8.985205, 60.711582 ], [ 8.985061, 60.711618 ], [ 8.984823, 60.711663 ], [ 8.984621, 60.711681 ], [ 8.984332, 60.711681 ], [ 8.98412, 60.71166 ], [ 8.983526, 60.711579 ], [ 8.982879, 60.711582 ] ])
circle = shapely.polygons([ [ 8.993146, 60.708983 ], [ 8.993075, 60.708632 ], [ 8.992866, 60.708294 ], [ 8.992526, 60.707983 ], [ 8.992069, 60.707711 ], [ 8.991512, 60.707487 ], [ 8.990876, 60.707321 ], [ 8.990187, 60.707219 ], [ 8.989469, 60.707184 ], [ 8.988752, 60.707219 ], [ 8.988063, 60.707321 ], [ 8.987427, 60.707487 ], [ 8.98687, 60.707711 ], [ 8.986413, 60.707983 ], [ 8.986073, 60.708294 ], [ 8.985864, 60.708632 ], [ 8.985793, 60.708983 ], [ 8.985864, 60.709334 ], [ 8.986073, 60.709671 ], [ 8.986413, 60.709982 ], [ 8.98687, 60.710255 ], [ 8.987427, 60.710478 ], [ 8.988062, 60.710644 ], [ 8.988752, 60.710747 ], [ 8.989469, 60.710781 ], [ 8.990187, 60.710747 ], [ 8.990876, 60.710644 ], [ 8.991512, 60.710478 ], [ 8.992069, 60.710255 ], [ 8.992526, 60.709982 ], [ 8.992866, 60.709671 ], [ 8.993075, 60.709334 ], [ 8.993146, 60.708983 ] ])

intersection = shapely.intersection(line, circle)
4
  • Could you share a reproductible example ? In my case if i use unary_union on a separated MultiLinestring, it stays a MultiLinestring Jan 25 at 18:21
  • @MagicMushroom I added the example the end of my question
    – Andreas
    Jan 25 at 18:34
  • In my opinion, th issue comes from the last part of the MultiLineString. If you separate every part of the MLS and ask of intersects or touches, the last (part 5) part will return a False with the 4th part. And you can figure out that part 4 and part 2 are actually the same. Jan 26 at 9:08
  • Do you have any idea why shapely.intersection produces a MultiLineString at all? The given LineString should not be cut into a MultiLineString by the intersection. I added an image to the question how the given LineString looks like (blue line). I guess the problem has something to do with the line going back and forth the same way and that some points on the way back are exactly the same as on the way forward? When I provide a line to shapely.intersection that has no duplicate points everything is fine, it generates a simple LineString!
    – Andreas
    Jan 26 at 9:15

2 Answers 2

0

The issue comes from your input Line wich is not a real Linestring because it is made up of superfluous points.

Your are in a particular case, if your take a closer look at your geometries there is a separation between 2 parts.

import shapely

line = shapely.geometry.LineString([ [ 8.985738, 60.712983 ], [ 
8.985523, 60.712769 ], [ 8.985278, 60.712175 ], [ 8.985441, 
60.711763 ], [ 8.985674, 60.711554 ], [ 8.986364, 60.711158 ], [ 
8.986783, 60.71076 ], [ 8.986913, 60.710638 ], [ 8.987753, 
60.710386 ], [ 8.987975, 60.710381 ], [ 8.98801, 60.710381 ], [ 
8.988991, 60.709666 ], [ 8.989304, 60.709225 ], [ 8.988991, 
60.709666 ], [ 8.98801, 60.710381 ], [ 8.987753, 60.710386 ], [ 
8.986913, 60.710638 ], [ 8.986364, 60.711158 ], [ 8.985674, 
60.711554 ], [ 8.985427, 60.711546 ], [ 8.985205, 60.711582 ], [ 
8.985061, 60.711618 ], [ 8.984823, 60.711663 ], [ 8.984621, 
60.711681 ], [ 8.984332, 60.711681 ], [ 8.98412, 60.71166 ], [ 
8.983526, 60.711579 ], [ 8.982879, 60.711582 ] ])
circle = shapely.geometry.Polygon([ [ 8.993146, 60.708983 ], [ 
8.993075, 60.708632 ], [ 8.992866, 60.708294 ], [ 8.992526, 
60.707983 ], [ 8.992069, 60.707711 ], [ 8.991512, 60.707487 ], [ 
8.990876, 60.707321 ], [ 8.990187, 60.707219 ], [ 8.989469, 
60.707184 ], [ 8.988752, 60.707219 ], [ 8.988063, 60.707321 ], [ 
8.987427, 60.707487 ], [ 8.98687, 60.707711 ], [ 8.986413, 
60.707983 ], [ 8.986073, 60.708294 ], [ 8.985864, 60.708632 ], [ 
8.985793, 60.708983 ], [ 8.985864, 60.709334 ], [ 8.986073, 
60.709671 ], [ 8.986413, 60.709982 ], [ 8.98687, 60.710255 ], [ 
8.987427, 60.710478 ], [ 8.988062, 60.710644 ], [ 8.988752, 
60.710747 ], [ 8.989469, 60.710781 ], [ 8.990187, 60.710747 ], [ 
8.990876, 60.710644 ], [ 8.991512, 60.710478 ], [ 8.992069, 
60.710255 ], [ 8.992526, 60.709982 ], [ 8.992866, 60.709671 ], [ 
8.993075, 60.709334 ], [ 8.993146, 60.708983 ] ])
intersection = shapely.intersection(line, circle)

g1 = [*intersection.geoms][0]
g2 = [*intersection.geoms][1]
g3 = [*intersection.geoms][2]
g4 = [*intersection.geoms][3]
g5 = [*intersection.geoms][4]

print(g1.touches(g2))
print(g2.touches(g3))
print(g3.touches(g4))
# issue with part 5
print(g5.touches(g5))
print(g5.touches(g3))
print(g5.touches(g2))
print(g5.touches(g1))

enter image description here

If you really want it to be a LineString you can concatenate coords

coords = [np.array(i.coords) for i in [g1, g2, g5, g3, g4]]
shapely.geometry.LineString(np.concatenate(coords))

But be careful with the order of the geometries

0

The problem is that there are duplicate nodes (coordinates) and thus, duplicate segments, in your LineString (points in red)

from collections import Counter
duplicates = [node for node, count in Counter(list(line.coords)).items() if count > 1]
print(duplicates)
[(8.985674, 60.711554), (8.986364, 60.711158), (8.986913, 60.710638), (8.987753, 60.710386), (8.98801, 60.710381), (8.988991, 60.709666)]

Duplicate nodes

enter image description here

Duplicate segments

enter image description here

Therefore is it normal that the result of the intersection is a MultiLineString

print(type(line.intersection(circle)))
<class 'shapely.geometry.linestring.LineString'>

enter image description here

If you eliminate the duplicates (segments) the result is a LineString

print(type(line2.intersection(circle)))
<class 'shapely.geometry.multilinestring.MultiLineString'>

enter image description here

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.