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:
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?
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=}")