7

How to speed up the "union" of two shapefiles using geopandas?

I have two large shapefiles:

geodataframe_left.shape = (3610, 12) (500 MB) (link)

geodataframe_right.shape = (16396, 3) (200 MB) (link)

One file contains global province boundaries (GADM), The other shapefile contains Hydrobasin level 6. I've tried multiple approaches (PostGIS, Google Bigquery, ArcGIS) but nothing beats geopandas' one line command:

gdf_union = gpd.overlay(gdf_left, gdf_right, how='union')

there is however one problem. The command is extremely slow on larger datasets including mine. I've waited for almost an hour now without success. I'm running my code on a large amazon EC2 instance (32GB RAM).

How can I improve the performance of this operation? Is it possible to use multiprocessing for example or can I split the problem into smaller chunks so I can monitor progress?

UPDATE:

  • I tried updating to geopandas 0.4.0 .. still no result after running for two hours.
  • Used the use_sindex = True flag. Deprecated since 0.4.0 (always uses sindex)
  • I implemented the tile approach and found weird behavior. Details below.

I created a fishnet grid of polygons of 10 x 10 degrees that I use to clip both layers before performing the union. The processing times per cell differ wildly and there are 4 cells that take > 10 minutes to process. One cell takes even 40 minutes.

When the process runs in parallel I get the following error (script continues though):

AttributeError: 'IndexStreamHandle' object has no attribute '_ptr'
Exception ignored in: <bound method Handle.__del__ of <rtree.index.IndexStreamHandle object at 0x7f007cb42e10>>
Traceback (most recent call last):
  File "/opt/anaconda3/envs/python35/lib/python3.5/site-packages/rtree/index.py", line 875, in __del__
    self.destroy()
  File "/opt/anaconda3/envs/python35/lib/python3.5/site-packages/rtree/index.py", line 863, in destroy
    if self._ptr is not None:

enter image description here

What's interesting is that these cells are at the 180 meridian. I've uploaded the GPKG here

enter image description here

Maybe shapely doesn't like this hemisphere crossing stuff.

UPDATE 2:

Running ArcMap locally, the process takes appr. 10 minutes. Instead of using geopackages, I used shapefiles. file1 file2

Executing: Union "hybas_lev06_v1c_merged_fiona_V04 #;gadm36_1 #" C:\Users\Rutger.Hofste\Desktop\werkmap\union_benchmark\output\union_arcgis_global.shp ALL # GAPS
Start Time: Fri Nov 30 10:32:04 2018
Reading Features...
Processing Tiles...
Assembling Tile Features...
Succeeded at Fri Nov 30 10:43:02 2018 (Elapsed Time: 10 minutes 57 seconds)

I also tried "Union" (both default and SAGA) in QGIS but this process is also ridiculously slow. I ran for 30min and the progress bar was at 4%.

UPDATE 3: I'm implementing the tiled approach however for a small number of polygons the result of an intersection in shapely is not a multipologon or polygon but a geometrycollection with multiple geometries including lineStrings or just LineStrings. I would expect the result of an intersection to be polygon or multipolygon.

18
  • This seems to be very relevant: jorisvandenbossche.github.io/blog/2017/09/19/geopandas-cython
    – RutgerH
    Commented Nov 28, 2018 at 16:15
  • 1
    You can clip both layers into Tiles, Union the Tiles, then merge the results. More steps but can be run in-parallel.
    – klewis
    Commented Nov 28, 2018 at 18:00
  • 1
    The .overlay command has an optional 5th parameter, use_sindex, have you tried this?
    – klewis
    Commented Nov 29, 2018 at 16:06
  • 1
    @RutgerH I am afraid that the cythonized version of geopandas will not help much. I did a timing with a small subset of the data (subset in Europe, one of the parts that take long), and ~90% of the time is spent inside GEOS itself, the underlying C library that does the union/intersection/difference operations.
    – joris
    Commented Nov 29, 2018 at 23:10
  • 1
    @klewis use_sindex is deprecated in the 0.4 version of geopandas, as the new implementation of overlay always uses a spatial index.
    – joris
    Commented Nov 29, 2018 at 23:12

1 Answer 1

0

The test data is not available anymore, but the geofileops library specifically tries to speed up processing large files and has a union function.

Based on the benchmarking I did, the speedup factor to be expected is about the number of CPUs/cores you have available for relatively simple polygons. If complex polygons are involved (> 1000 points) the speedup can be significantly more.

Sample script:

import logging
import geofileops as gfo

if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO)
    gfo.union(
        input1_path="input1.gpkg",
        input2_path="input2.gpkg",
        output_path="output.gpkg",
    )

Disclaimer: I'm the developer of geofileops.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.