Geopandas - split lines at intersection with polygons, retain polygon ID in new link attributes

Viewed 2957

I have two geodataframes: one containing many Linestrings, and another containing many Polygons. These lines and polygons intersect each other. The output I am trying to achieve is a new geodataframe containing links which have been split at any location that they intersect a polygon.

Simplified test code is as follows:

import geopandas
from shapely.geometry import Polygon, LineString

#Create geodataframe containing 2 links:
linkID = ['1','2']
link_geom = [LineString([(0, 0), (10, 10)]),LineString([(10, 10), (20, 10)])]
a = {'linkID':linkID,'geometry':link_geom}
gdf_links = geopandas.GeoDataFrame(a)

#Create geodataframe containing 2 polygons:
polyID = ['100','200']
poly_geom = [Polygon([(2, 1), (2, 3), (4, 3), (4, 1)]),Polygon([(15, 7), (15, 13), (18, 13), (18, 7)])]
b = {'polyID':polyID,'geometry':poly_geom}
gdf_poly = geopandas.GeoDataFrame(b)

You can see that the lines intersect the polygons by running the following code:

links = gdf_links.unary_union
polys = gdf_poly.unary_union
geopandas.GeoSeries([links,polys]).plot(cmap='tab10')

Which generates the following figure (annotation added in red/green to show the various elements involved):

enter image description here

Splitting the two lines wherever they intersect a polgyon will results in 6 segments (as labelled in green in the figure above).

Ideally, the output I would like to achieve is something like shown below, with a 'poly_intersect' column listing the ID of the polygon at which the split has occurred. Note that items 2 and 5 list two polygon ID's, because they intersect the polygon twice.

enter image description here

1 Answers

You can get those linestrings by doing an "overlay" of the lines with the polygons (using the "union" operation while overlaying). However, it will not give exactly the expected result as you depicted for the "poly_intersect" column:

>> geopandas.overlay(gdf_links, gdf_poly, how='union').explode().reset_index(drop=True)

  linkID polyID                   geometry
0      1    100      LINESTRING (2 2, 3 3)
1      2    200  LINESTRING (15 10, 18 10)
2      1    NaN      LINESTRING (0 0, 2 2)
3      1    NaN    LINESTRING (3 3, 10 10)
4      2    NaN  LINESTRING (10 10, 15 10)
5      2    NaN  LINESTRING (18 10, 20 10)

Some remarks:

  • I am using explode(), because the overlay() operation returns MultiLineStrings for the part of one linestring before and after intersecting the polygon. The explode will split those "multi" geometries.
  • The "polyID" column indicates which linestrings originated from intersecting with which one the polygons. But it does not have information for the other linestrings of which polygon they touch.

So if you want require that information as well (the full "poly_intersect" column), you could do some follow-up analysis checking which of the linestrings is touching which of the polygons.

Related