# Fast quasi-recursive polygon intersection algorithm?

I'm trying to make a map which is colored based on the distance to a target point. I've done this by creating a series of concentric rings.

These rings are created by buffering out from the target point a set distance, then removing the previous distance circle, like so:

``````gdf4['geometry'] = gdf4['geometry'].buffer((1609.344 * 2))
gdf4['buffer'] = 2

for i in range(2, 20, 2):
gdfx = gdf.copy()
gdfx[f'geometry_{i}'] = gdfx['geometry'].buffer((1609.344 * i))
gdfx[f'geometry_{i + 2}'] = gdfx['geometry'].buffer((1609.344 * (i + 2)))

for center in list(gdf['name'].unique()):
big  = gdfx.loc[gdfx['name'] == center, f'geometry_{i + 2}'].iloc[0]
small = gdfx.loc[gdfx['name'] == center, f'geometry_{i}'].iloc[0]
gdfx.loc[gdfx['name'] == center, 'geometry'] = big.difference(small)

gdfx['buffer'] = i + 2

gdf4 = pd.concat([gdf4, gdfx])
``````

The rings look like this:

The problem I'm having is how to color the polygons created when the rings intersect. I'd like to color these intersects based on the point they are closer to (ie. if a 2 ring intersects with a 4 ring, I want the intersect to be colored the 2 color). This alone is straightforward. The problem arises when I have intersects of intersects like here around Chicago:

Here's what I have so far:

``````def find_intersection_polys(gdf):
spatial_index = gdf.sindex

processed_polys = set()

while True:
new_intersect_gdf = gpd.overlay(gdf, gdf, how='intersection')

new_intersect_gdf = new_intersect_gdf[~new_intersect_gdf.index.isin(processed_polys)]

if new_intersect_gdf.empty:
break

new_intersect_gdf['buffer_value'] = new_intersect_gdf.apply(lambda row: min(row['buffer_1'], row['buffer_2']), axis=1)

processed_polys.update(new_intersect_gdf.index)

gdf = gpd.GeoDataFrame(pd.concat([gdf, new_intersect_gdf]), crs=gdf.crs)

return gdf

intersect_gdf = find_intersection_polys(gdf4)
``````

Obvious the implementation here (while-break) is far too slow. I tried running this for several hours yesterday and it still didn't finish.

Can anyone suggest something faster/better?

• 614 points total. Jul 21, 2023 at 16:17
• Another approach would be to create an euclidian distance raster then ( if you need vector data) reclass/vectorize it
– J.R
Jul 21, 2023 at 16:59

This is based on Kevin Neufeld answer of a PostGIS problem where polygons overlap.

You create lines of the polygon boundaries, union them into a big multiline, polygonize it and transfer attributes from the original polygon table to the new polygons by using points inside polygons.

It finishes in 16 s. for 1000 multibuffered points.

``````import geopandas as gpd
import shapely
from timeit import default_timer

start = default_timer()

#1. Create lines from the polygon boundaries
lines = df.apply(lambda x: x["geometry"].boundary, axis=1).to_list()

#2. Union all the lines to into one multiline.
noded_lines = shapely.ops.unary_union(lines)

#3. Re-polygonize all lines
noded_lines_singleparts = [x for x in noded_lines.geoms] #Create a list of singlepart lines
new_polys = list(shapely.polygonize(noded_lines_singleparts).geoms) #Create polygons from it

#4. Transfer attributes from original polygons to new polygons, by using points in the polygons
new_polys = gpd.GeoDataFrame(geometry=new_polys, crs=df.crs) #Create a dataframe from the list of new polygons
new_polys_point = new_polys.copy() #Create a copy of it to use to create points in polygons
new_polys_point["geometry"] = new_polys.geometry.representative_point()

#Overlay points with original polygons, transferring the attributes (the ring number)
points_with_attributes = gpd.sjoin(left_df=new_polys_point, right_df=df, how="left", predicate="intersects")

#Join the points to the newly formed polygons, transferring the attributes (ring number):
del(points_with_attributes["index_right"])
new_polys_with_attributes = gpd.sjoin(left_df=new_polys, right_df=points_with_attributes,
how="left", predicate="intersects")
new_polys_with_attributes = new_polys_with_attributes.loc[~new_polys_with_attributes["ringnum"].isna()]

#Drop duplicate geometries by geometry and ring number.
new_polys_with_attributes["wkt"] = new_polys_with_attributes.apply(lambda x: x.geometry.wkt, axis=1)
new_polys_with_attributes = new_polys_with_attributes.sort_values(by=["wkt", "ringnum"], ascending=[True, False])
new_polys_with_attributes["ringnum"] = new_polys_with_attributes["ringnum"].astype(int)
new_polys_with_attributes = new_polys_with_attributes.drop_duplicates(subset="wkt", keep="last")
del(new_polys_with_attributes["wkt"])

#Dissolve boundaries
dissolved = new_polys_with_attributes.dissolve(by="ringnum")

#Dissolve creates multipolygons, convert to singleparts
dissolved["singleparts"] = dissolved.apply(lambda x: [poly for poly in x.geometry.geoms], axis=1) #List all parts
dissolved = dissolved.explode("singleparts").set_geometry("singleparts") #Explode so each part becomes one row
dissolved = dissolved.drop(columns=["geometry"]).rename_geometry("geometry").set_crs(df.crs)
dissolved.to_file(r"/home/bera/Desktop/GIStest/multirings_clean/new_polys_with_attributes_dissolved.shp")
print(f"Finished processing in {round(default_timer()-start)} s.")

#Finished processing in 16 s.
``````

• This totally worked, thanks so much!!! Jul 26, 2023 at 19:08