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