Splitting many polygons by many lines

Question

I have many polygons that need to be cut by lines we have set up to split those polygons. Some polygons may not have any lines to split them, others will have only 1 line, while others may have several line splits.

Here are a couple example screenshots. First one, ploy 406 has no lines to split it, poly 25 has 3 lines that will split it, an finally in poly 27 has only 1 line to split it. All lines have been checked to make sure they completely cross all polys.

Im using python with geopandas and shapely libraries to handle this.

So far I have come up with a solution that appears to identify where each line intersects with a polygon and splits it with shapely.ops.split. But once it is complete the geometries are separated from the original geodataframe and Id like to maintain the attributes from the polygons.

# iterate over each cutline and find intersecting polygon
split_lots = []
for index, crow in df_cutlines.iterrows():
    for index, prow in df_parkinglots.iterrows():
        if crow['geom'].intersects(prow['geom']):
            geometry = split(prow['geom'],crow['geom']) # split poly with line
            split_lots.append(geometry)

There must be a better way to iterate through each line to find the intersecting poly, cut it and keep the attributes.

I thought about using the apply function to go through each row and add the new geom as another column to the original polygons (this would replace my last 2 lines above).

df_parkinglots['split_geom'] = df_parkinglots.apply(lambda prow: split(prow['geom'],crow['geom']), axis=1)

But shapely doesnt seem to like this.

For now I dont really care if the cut polys end up being multipart or single.

Answer 1

You can list all intersecting lines, split the polygon with each line, extract their boundaries, union, and polygonize. Then explode the dataframe to create a row for each polygon part:

import geopandas as gpd
from shapely.ops import split, unary_union
from shapely import polygonize

db = r"/home/bera/Desktop/GIStest/split.gpkg"
poly = gpd.read_file(filename=db, layer="polygon")
poly["polyid"] = range(poly.shape[0]) #Create a unique id column

line = gpd.read_file(filename=db, layer="line")
line["linegeom"] = line.geometry

#Create a column of all lines intersecting each polygon
sj = gpd.sjoin(left_df=poly, right_df=line[["linegeom","geometry"]], how="left")
poly["splitlines"] = sj.groupby("polyid")["linegeom"].apply(list)


def splitpoly(polygeom, linelist):
"""A function to split an input polygon geometry by a list of linestrings"""
    a = [list(split(polygeom, line).geoms) for line in linelist] #List of lists of split polygons
    b = [item for sublist in a for item in sublist] #Flattened into a list of polygons
    boundaries = [x.boundary for x in b]+[polygeom.boundary]
    uu = unary_union(boundaries)
    split_polygon = list(polygonize(list(uu.geoms)).geoms) #A list of polygons. Possible to multipolygonize them he
    return split_polygon

#Create a column with all the parts of the split geometry
poly["polylist"] = poly.apply(lambda x: splitpoly(x["geometry"], x["splitlines"]), axis=1)

result = poly.explode(column="polylist")
result = result.set_geometry(col="polylist", crs=poly.crs).drop(columns="geometry").rename_geometry("geometry")
del(result["splitlines"])

result = result[result.geometry.area>1] #There are some slivers, drop the polygons with an area < 1

result["newid"] = range(result.shape[0])
result.to_file(r"/home/bera/Desktop/GIStest/split_polygons_result.shp")