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?
- 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:
What's interesting is that these cells are at the 180 meridian. I've uploaded the GPKG here
Maybe shapely doesn't like this hemisphere crossing stuff.
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.