pythongeopandasintersectionshapely

Problems obtaining an intersected linestring using GeoPandas


I am having an issue obtaining a linestring from an intersection with a polygon using GeoPandas. The linestring is self-intersecting, which is what is causing my issues.

A line intersecting a polygon:

enter image description here

Given the following code:

import geopandas as gp
from shapely.geometry import LineString, Polygon

# Draw a polygon that is 100 x 100 units, starting at coordinates 0, 0
polygon = Polygon([(50, 0), (50, 100), (100, 100), (100, 0)])

# Convert the polygon to a geodataframe
polygon = gp.GeoDataFrame(index=[0], crs='EPSG:4326', geometry=[polygon])

# Draw a horizontal line that starts at coordinates 50, 0 and is 200 units long
line = LineString([(0, 50), (75, 50), (70, 35), (55, 40), (250, 50)])

# Convert the line to a geodataframe
line = gp.GeoDataFrame(index=[0], crs='EPSG:4326', geometry=[line])

print(line)

# Intersect the line with the polygon
intersection = line.intersection(polygon)

print(intersection)

I have the following results:

0 LINESTRING (0.00000 50.00000, 75.00000 50.0000... 0 MULTILINESTRING ((50.00000 50.00000, 75.00000 ...

After intersect, I am returned a multilinestring instead of a linestring. The line is being split by the polygon (desired) but is also being split into multiple lines where it self-intersects (not desired).

I have tried to re-join the multiline with unary-union without success. The output remains a multilinestring.

I'm unsure what else I can do to keep only the portion of the line contained within the polygon, as a single line. Any thoughts on how I might be able to accomplish this?


Solution

  • Because line_merge apparently doesn't reconstruct the single linestring when it is self-intersecting, the only option I see is to extract the coordinates and create a new linestring, which results in a single linestring.

    In the code sample I use some functions that are only available in shapely 2, so it needs a relatively up-to-date version of shapely.

    Code sample:

    import geopandas as gp
    import shapely
    from shapely.geometry import LineString, Polygon
    
    # Draw a polygon that is 100 x 100 units, starting at coordinates 0, 0
    polygon = Polygon([(50, 0), (50, 100), (100, 100), (100, 0)])
    
    # Convert the polygon to a geodataframe
    polygon = gp.GeoDataFrame(index=[0], crs='EPSG:4326', geometry=[polygon])
    
    # Draw a horizontal line that starts at coordinates 50, 0 and is 200 units long
    line = LineString([(0, 50), (75, 50), (70, 35), (55, 40), (250, 50)])
    
    # Convert the line to a geodataframe
    line = gp.GeoDataFrame(index=[0], crs='EPSG:4326', geometry=[line])
    
    # print(line)
    
    # Intersect the line with the polygon
    intersection = line.intersection(polygon)
    
    print(f"{intersection=}")
    
    # Recreate the intersection line from its coordinates
    intersection_line = shapely.linestrings(shapely.get_coordinates(intersection))
    
    print(f"{intersection_line=}")