Intersection in geopandas takes long time to execute

Question

I got two dataframes in GeoPandas Python. One table contains around 2.2 million rows of data (geometry type is Multipolygon) and other table contains around 3800 rows (geometry type is Multiplygon). I am trying to calculate how many polygons from bigger table are either completely 'within' smaller table's polygon or if they intersect with each other how much area does it overlap with other table's polygon. Following is the code I have written:

import geopandas as gpd
import pandas as pd

with_in = gpd.sjoin(parcels_gdf, coverage_df, how='inner', predicate='within')
with_in['Full_covered'] = 100
remaining_parcels  = parcels_gdf.drop(with_in.index)
intersections = remaining_parcels.intersection(coverage_df.unary_union)
intersection_areas = intersections.area
total_intersection_area = intersection_areas.sum()

parcels_gdf is the the table that contains 2.2 million rows. coverage_df contains 3800 rows. remaining_parcels contains around 1.5 million rows. The issue I have, program is taking very long (more than 12 hours as I write and still running) when it execute intersections = remaining_parcels.intersection(coverage_df.unary_union). I am not sure how long further it takes to compelte the execution. I got a laptop with core i7 with 16 GB. Is there any better way to program it faster?

Answer 1

This: intersections = remaining_parcels.intersection(coverage_df.unary_union) is very slow,

because each polygon in remaining_parcels is intersected with a huge dissolved/unioned multipolygon.

Try this:

import geopandas as gpd

bigdf = gpd.read_file("/path/to/file")
bigdf["bigid"] = range(bigdf.shape[0])
smalldf = gpd.read_file("/path/to/file2")

within = gpd.sjoin(left_df=bigdf, right_df=smalldf, predicate="within")
within["full_coverage"] = 100

#Intersect the polygons in bigdf which are not within, with the smalldf.
inter = gpd.overlay(df1=bigdf.loc[~bigdf.index.isin(within.index)],
            df2=smalldf, how="intersection", keep_geom_type="True")
inter["area"] = inter.geometry.area

inter.groupby("bigid")["area"].sum() #If you want each bigdf's intersected polygons area
inter.area.sum() #Or the total

Answer 2

There are several options that will/should be faster.

geopandas.clip

geopandas.intersection doesn't use a spatial index to check first if it is useful to calculate the intersection before calculating it, which can make a big difference when a significant number of the remaining_parcels don't overlap. The geopandas.clip tool will optimize this case. I don't know you data, but I would first look at the next options because I think they will be a better fit.

To try this, replace

intersections = remaining_parcels.intersection(coverage_df.unary_union)

with

intersections = remaining_parcels.clip(coverage_df)

geopandas.overlay

Depending on your specific needs another option might be to use geopandas.overlay to calculate the intersections between the two dataframes. Because in this case coverage_df doesn't need to be unioned first it will be able to use the spatial index efficiently and the intersections to be calculated will be a lot easier. Disadvantage is that if an input feature intersects with multiple rows in coverage_df, this will result in the intersections being divided in multiple rows. If you only need the total intersection area like in your sample above though this is no problem.

intersections = remaining_parcels.overlay(coverage_df, how='intersection')

geofileops.clip

Disclaimer: I am the developer of the Geofileops library.

Finally, you could try using geofileops. This library also has a clip functionality but it has some extra optimizations and will use multiprocessing to speed it up. Possible disadvantage is that the input data should be in a GeoPackage file.

import logging
import geofileops as gfo

if __name__ == "__main__":
    # Init logging so progress printed by gfo is shown
    logging.basicConfig(level=logging.INFO)

    # Write input to file if it only exists in memory
    parcels_path = "/path/you/like/parcels.gpkg"
    coverage_path = "/path/you/like/coverage.gpkg"
    gfo.to_file(parcels_gdf, parcels_path)
    gfo.to_file(coverage_df, coverage_path)

    # Calculate clip
    output_path = "/path/you/like/clip.gpkg"
    gfo.clip(input_path=parcels_path, clip_path=coverage_path, output_path=output_path)