Iterate through geodataframes and overlay them to create a single output

Question

I am attempting to iterate through dissolved_result and overlay them using how='intersection' to create a single output file that represents all of the input dissolved_results. I am struggling with how best to do that, this is what I have so far and the error I get.

All inputs are polygons/multipolygons.

Edit: added some links to snaps to show what I'm trying to do. I have these three shapes, I want to intersect them all at the end of this function.(One of them appears to have no features, they're just not visible at this scale.)

def usable(gdf, feature_gdf, zoom_level, shuffle_tiles=False):
feature_gdf = gpd.read_file(inshp)
feature_gdf = feature_gdf.to_crs(epsg=5070)   
intermediate_result = gpd.GeoDataFrame()

for index, row in gdf.iterrows(): 
    geometry = row['geometry']
    tiles = list_tiles(geometry, zoom=zoom_level, shuf=shuffle_tiles)

    for tile in tiles:
        tile_bounds = mercantile.bounds(*tile)
        tile_box = gpd.GeoSeries([box(*tile_bounds)], crs=gdf.crs)
        tile_box = tile_box.to_crs(epsg=5070)       

        if feature_gdf.intersects(tile_box.unary_union).any():
            intersected_features = feature_gdf[feature_gdf.intersects(tile_box.unary_union)]
            buffered_features = intersected_features.copy()
            buffered_features['geometry'] = buffered_features.buffer(buffer_distance)
            tile_minus_buffered = tile_box.difference(buffered_features.unary_union)
            tile_buffer = gpd.GeoDataFrame(geometry=tile_minus_buffered, crs=tile_box.crs)
            intermediate_result = pd.concat([intermediate_result, tile_buffer], ignore_index=True)

        else:
            tile_gdf = gpd.GeoDataFrame(geometry=tile_box, crs=tile_box.crs)
            intermediate_result = pd.concat([intermediate_result, tile_gdf], ignore_index=True)

    # save_file(intermediate_result, 'result')

dissolved_result = intermediate_result.dissolve()

final_result = gpd.GeoDataFrame()

for index, row in dissolved_result.iterrows():
    result = gpd.GeoDataFrame(geometry=[row['geometry']], crs=dissolved_result.crs)   
    final_result = gpd.overlay(final_result, result, how='intersection')
print(final_result.head())



Error (from final for-loop):
AttributeError: The CRS attribute of a GeoDataFrame without an active geometry column is not defined. Use GeoDataFrame.set_geometry to set the active geometry column.

I dont understand what you want to do. Can you add some screenshots or drawings to make it clearer? — BERA, Oct 31 '23 at 06:03
BERA - Certainly, I added some links. I don't have enough reputation yet to embed them inline. I am trying to iteratively intersect the three images, though the final product will have an arbitrary number of shapes to intersect. Does that help? — maharg, Oct 31 '23 at 15:13
Do you have only one dataframe and want to find the common intersection of all features in it? — BERA, Oct 31 '23 at 18:12
No, I have three output dataframes (each one is represented in the pictures above) and I want to intersect all three. I understand it will have to be two at a time, since that's how geopandas.overlay works. I'm just struggling with how best to implement that. — maharg, Oct 31 '23 at 20:43
@BERA - Thanks for your patience. This was very helpful! I had to change the 'how' to 'how='intersection'' and now I get the results I was looking for. Do you have any insight into how to use this process for an arbitrary number of dfs? — maharg, Nov 13 '23 at 15:46

BERA · Accepted Answer · 2023-11-01T18:30:03.557

You can use functools.reduce to union all frames, then select:

import geopandas as gpd
from functools import reduce
df = gpd.read_file(filename=r"C:/GIS/GIStest/10k_buffered_points.geojson")
df["id1"] = 1 #A temp id column
df2 = gpd.read_file(filename=r"C:/GIS/GIStest/swamp.gpkg")
df2["id2"] = 1
df3 = gpd.read_file(filename=r"C:/GIS/GIStest/open_ground.gpkg")
df3["id3"] = 1
unioned = reduce(lambda frame1, frame2: gpd.overlay(df1=frame1, df2=frame2, how="union", keep_geom_type=True),
           [df, df2, df3]) #Union df1+df2, then union the output with df3.
unioned[["id1","id2","id3"]].sample(3)
id1  id2  id3
1.0  NaN  1.0       #df1 and df3 intersecting
1.0  NaN  NaN       #Only features from df1
1.0  1.0  NaN       #df1 and df2 intersecting
intersecting = unioned.loc[unioned[["id1","id2","id3"]].sum(axis=1)>1] #Select rows where two frames intersects (rowwise sum >1)
intersecting.to_file(r"C:/GIS/GIStest/intersecting.gpkg")

You might want to dissolve the results.

Iterate through geodataframes and overlay them to create a single output

1 Answers1

unioned[["id1","id2","id3"]].sample(3)

id1 id2 id3

1.0 NaN 1.0 #df1 and df3 intersecting

1.0 NaN NaN #Only features from df1

1.0 1.0 NaN #df1 and df2 intersecting