1

In using 'intersection' to crop / clip a single 'braided' line using a rectangle (actually a grown version of that rectangle using 'buffer'), the internal portion of the line gets split on self-intersections. The result of the intersection operation is a MultiLineString.

enter image description here

I would like to avoid that, and just get one internal 'braided' line. Why does the intersection method do this? Is there a way to control it? Is there a quicker way to end up with one single line than joining the LineStrings end-to-end after the intersection operation?

The responses to Welding individual line segments into one LineString using Shapely seem to indicate that shapely.ops.linemerge should be the way to go - but it doesn't seem to do the trick, maybe there's something I'm misunderstanding:

enter image description here

from shapely.geometry import LineString,Polygon,MultiLineString,MultiPolygon,GeometryCollection
from shapely.ops import split,linemerge
...
result=targetGeom&boundaryGeom # could be MultiPolygon or MultiLinestring or GeometryCollection
logging.info('crop result:'+str(result))
result=linemerge(result)
logging.info('after linemerge:'+str(result))

Running the code on the above shapes:

INFO:root:crop: target=c1  boundary=c2
INFO:root:crop result:MULTILINESTRING ((-121.108305669868 39.29086839149927, -121.1084962057916 39.29034474676413, -121.1080157152655 39.29023810613076), (-121.1080157152655 39.29023810613076, -121.1070370840875 39.29002090757433, -121.1083674597589 39.28961403057203, -121.1080157152655 39.29023810613076), (-121.1080157152655 39.29023810613076, -121.1076812116669 39.29083159238037))
INFO:root:after linemerge:MULTILINESTRING ((-121.108305669868 39.29086839149927, -121.1084962057916 39.29034474676413, -121.1080157152655 39.29023810613076), (-121.1080157152655 39.29023810613076, -121.1076812116669 39.29083159238037), (-121.1080157152655 39.29023810613076, -121.1070370840875 39.29002090757433, -121.1083674597589 39.28961403057203, -121.1080157152655 39.29023810613076))

(shorthand translation - there is no rounding error, the endpoints are exact matches):

INFO:root:crop result:MULTILINESTRING ((A, B, C), (C, D, E, C), (C, F))
INFO:root:after linemerge:MULTILINESTRING ((A, B, C), (C, F), (C, D, E, C))

It looks like the last two segments of the multilinestring are swapped, but otherwise it's still the same three unmerged segments:

enter image description here

UPDATE: I guess this is really a question of how to preserve complex / non-simple lines (is_simple=False) during the geometry operations. Found the same question here from several years back, which had no great solution other than brute-force.

5
  • Can you show the wanted result? Are those linestrings valid (you mention self intersections...)? Jul 3, 2021 at 17:15
  • graphically, the result would be the same as the current result, except you'd only have one label since there would only be one line. Shortening the output above (I did this with search-and-replace in vi to avoid typos - also it shows there is no rounding error - the points are exact): INFO:root:crop result:MULTILINESTRING ((A, B, C), (C, D, E, C), (C, F)) INFO:root:after linemerge:MULTILINESTRING ((A, B, C), (C, F), (C, D, E, C))
    – Tom Grundy
    Jul 3, 2021 at 17:58
  • Maybe the reason it's dying is because the crossover point (C in the list above) is a four-way intersection, so there's no one clear winner as to which segment should come 'next'.
    – Tom Grundy
    Jul 3, 2021 at 18:01
  • added the shorthand notation to the original question
    – Tom Grundy
    Jul 3, 2021 at 18:04
  • sorry, to answer your initial questions, it's valid but complex - added this in an update to the original post.
    – Tom Grundy
    Jul 4, 2021 at 0:22

1 Answer 1

0

You might be able to use the complex_split function from here (or here) and test the resulting geometries for which are inside of the rectangle. Like this:

# Split the braid on the rectangle w/o splitting on self-intersections
geom_collection = complex_split(braid, rectangle)

# Filter for the result inside the rectangle
linestrings_within = [ls for ls in geom_collection if rectangle.contains(ls)]

# We are only expecting one inside
assert len(linestrings_within) == 1

clipped_braid = linestrings_within[0]

Your Answer

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

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