Ways to speed up geopandas overlay when applying it on large datasets

Question

I have a large grid of 1km square cells covering an the area around the Atlantic Ocean around the Atlantic Provinces of Canada from Labrador to Greenland south to the US Canada border. In all it is about 4.3 million cells. I have 400k tracks representing shipping; about 10Gb of data. I am trying to overlay the grid onto the tracks to generate statistics for each cell. I have been passing it to the multiprocessing objects with various results using methods such as using numpy array_split to break the tracks and the grid into smaller chunks which worked for a small sample set but didn't for the whole dataset. I even tried doing it to one cell at a time which was working but was taking too long. I started looking into other methods like...

def cell_calc(df1, df2):
    applied_df1 = df1.apply(lambda row: gpoverlay(df1, df2).groupby(by='grid_id').size().reset_index(name='ALL').fillna(0), gpoverlay(df1, df2, how='intersection') 
   [gpoverlay(df1, df2) 
   ['type']=='CARGO'].groupby(by='grid_id').size().reset_index(name='CARGO').fillna(0))
   return concat(df1, applied_df1)

partialTask = ft.partial(cell_calc, df1=GRIDGDF, df2=TRACKSGDF)
    with mp.Pool() as pool:
        count = 0
        stime = ttime()
        results_gdf = pool.apply(partialTask)
        pool.close()

    print(results_gdf)

Which returned TypeError: 'DataFrame' objects are mutable, thus they cannot be hashed

and...

partialTask = ft.partial(gp.overlay, df1=GRIDGDF, df2=TRACKSGDF, how='intersection')
with Pool() as pool:
    results = pool.apply(partialTask)
    pool.close()
print(results)

Of which works on a test sample but when applied to the whole dataset it runs, then all the processes go to 0% excepts for one which sits at 6% to 7% and it never seems to come out of it.

Simply put, what methods work to speed up the application of geopandas overlay on large datasets and what are some working examples for applying them?

Answer 1

I ended up using multiprocessor and numpy array_split to pass sections of the data to an overlay. In order to get by the 2Gb limit I did all the grid math and returned a final dataframe to be concatinated into one dataframe and saved.

from geopandas import read_file, overlay
from pandas import merge, concat
from multiprocessing import Pool
from functools inport partial
from numpy import array_split

def func(df1, df2):
    df = overlay(df2.to_crs(df1.crs), df1, how='intersection')
    df = df.groupby(by='id').size().reset_index(name='Count').fillna(0)
    return(df1, df, on='id', how='left').fillna(0)[['cell_id', 'Count', 'geometry']]

if __name__ == '__main__':
    grid = #Grid Data has unique id as 'cell_id'
    tracks = #Track Data
    n = len(grid)//10000
    grid_results = []

    if n > 8:
        x = 8
    else:
        x = n
    pfunc = partial(func, df2=tracks)
    with Pool(x) as pool:
        for results in pool.map(pfunc, array_split(grid, n)):
        grid_results.append(results)
    pool.close()
    final_grid = concat(grid_results)

It takes about 15min to do the overlay and another 15 to 20 to do the cell calculations (I only included 1 count function, I have 78). In all there are 430 sets of 10000 grid squares and when using 3 processors I end up with an estimated time to completion of (((430*30)/3)/60)/24 = 2.9 days. I am looking at applying pandarallel's parallel_apply to speed up the cell math. If anyone has a solution using dask I would like to try it.

Answer 2

After some more testing it actually works better if you array_split the target df and not the mask df.

def mp_overlay(df2, df1):
    func = partial(overlay, df2=df2)
    overlays = []
    n = 16 #Number of workers you would like to use...
    with Pool(n) as pool:
        for overlay in pool.map(func, array_split(df1, n)):
            overlays.append(overlay)
    return concat(overlays)

if __name__ == '__main__':
    from pandas import concat, merge
    from geopandas import, read_file, overlay
    from functools import partial
    from multiprocessing import Pool
    from numpy import array_split

    df1 = read_file(target_date)
    df2 = read_file(mask_data)
    x = 8 #number of splits

    foverlays = []
    for grid in array_split(df2, x):
        df = mp_overlay(grid, df1.to_crs(grid.crs))
        foverlays.append(df)
        del df
    overlay = concat(foverlays)

My first answer took about 16h per array split and I had 430 of them with 10000 cells each. (len(grid)//10000 equals roughly 430 grids; if in every 16h I processed 8 it would have taken 35 days. The reason why it was faster in my testing is because I was overlaying a target that was clipped to my mask. As soon as the target had a large number of features outside the mask performance dropped. I then tested an overlay with a mask of 100 cells, 1000 cells, and 10000 cells and found that more cells didn't have a massive effect on the time to perform an overlay but more target features did so I swapped the target features to the process pool.

This also works if your PC resources cannot handle the results of the overlay, like mine, so each return could be an export to file and then each file can be worked with independently. In my case, I did some grid math using the overlay results which were much smaller than the results of the overlay so I could pass them back to the main worker, concatenate them, and export them to a file.

Answer 3

Disclaimer: I am the developer of the Geofileops library.

You could try using geofileops. This is a library that aims to efficiently process larger geo files. It will apply multiprocessing under the hood. For some operations geopandas is used, but not for overlay operations like intersection. I mainly work with polygon data... but line data should normally work as well.

The input data should be in Geopackage format.

import logging
import geofileops as gfo

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

    # Convert input data to Geopackage
    gfo.copy_layer(target_data, target_data_gpkg)
    gfo.copy_layer(mask_data, mask_data_gpkg)

    # Calculate intersection
    gfo.intersection(
        input1_path=target_data_gpkg,
        input2_path=mask_data_gpkg,
        output_path=output_path,
    )