Fast selection of intersecting geometries in GPKG file

Question

I have a GPKG file, example.gpkg, with 2 layers, base and circles. There is only 1 geometry in base, and there is a large number of geometries in circles. I would like to select all the geometries in circles which intersect with the single geometry in base. Below I present the two ways to achieve this.

I would like to learn how to make this process faster when circles has a large number of geometries.

---

Simple example

Let base have a polygon near origin, and let circles be 500 circles near origin. I generate such GeoSeries using Python:

import shapely.geometry
import geopandas as gpd
import numpy as np
np.random.seed(42)

base=gpd.GeoSeries(shapely.geometry.Point(0,0).buffer(1,resolution=3))

circles = gpd.GeoSeries(
    [
        shapely.geometry.Point(e)
        for e in np.random.uniform(low=-10, high=10, size=[500, 2])
    ]
).buffer(distance=0.1)

Visualize result:

import matplotlib.pyplot as plt
fig, ax = plt.subplots()

base.plot(ax=ax, edgecolor="red",facecolor="red")
circles.plot(ax=ax, edgecolor="black",facecolor="none")

Write result to GPKG file:

base.to_file("example.gpkg",layer="base")
circles.to_file("example.gpkg",layer="circles")

We can use ogr2ogr and ST_INTERSECTS to write file intersecting_circles.gpkg with the circles geometries intersecting base:

ogr2ogr -sql "SELECT circles.geom FROM circles, base WHERE ST_INTERSECTS(circles.geom, base.geom)" intersecting_circles.gpkg example.gpkg
 

It might also be possible to use spatial indexes to make this query faster. There is an rtree_circles_geom table in the GPKG file (this thread explains how we can know that). rtree_circles_geom contains the bounding boxes of each geometry in circles. Let's manually filter the geometries so that we only need to run ST_INTERSECTS on those circles which have bounding boxes overlapping with base's bounding box.

To achieve this, I create intersecting_circles_with_rtree.gpkg via:

extent=$(ogrinfo -sql "SELECT ST_Envelope(geom) FROM base" example.gpkg -json | jq -r '.layers[].geometryFields[].extent' |  tr -d '[] ')

l=(${extent//,/ })

base_minx=${l[0]}
base_miny=${l[1]}
base_maxx=${l[2]}
base_maxy=${l[3]}


ogr2ogr -sql \
"SELECT geom_inrange FROM ( \
    SELECT fid, geom as geom_inrange FROM circles \
        WHERE fid IN ( \
            SELECT id FROM rtree_circles_geom \
            WHERE maxx>=${base_minx} \
            AND minx<=${base_maxx} \
            AND maxy>=${base_miny} \
            AND miny<=${base_maxy} \
        )\
    ), base WHERE ST_INTERSECTS(geom_inrange,base.geom)" \
intersecting_circles_with_rtree.gpkg example.gpkg

I open both intersecting_circles.gpkg and intersecting_circles_with_rtree.gpkg in QGIS. They look the same: they both contain the circles which intersect the red polygon (see visualization above). Two different methods, same result: below I investigate which method selects the right geometries faster.

---

Benchmarking

I would like to know how the time needed to select the intersecting geometries depends on the number of geometries in circles. Since I care about selecting the right polygons, not writing a new GPKG file, I am going to use ogrinfo instead of ogr2ogr (and I'll suppress the non-error outputs by using > /dev/null, as described here).

I create query.sh:

extent=$(ogrinfo -sql "SELECT ST_Envelope(geom) FROM base" example.gpkg -json | jq -r '.layers[].geometryFields[].extent' |  tr -d '[] ')

l=(${extent//,/ })

base_minx=${l[0]}
base_miny=${l[1]}
base_maxx=${l[2]}
base_maxy=${l[3]}


ogrinfo -sql \
"SELECT geom_inrange FROM ( \
    SELECT fid, geom as geom_inrange FROM circles \
        WHERE fid IN ( \
            SELECT id FROM rtree_circles_geom \
            WHERE maxx>=${base_minx} \
            AND minx<=${base_maxx} \
            AND maxy>=${base_miny} \
            AND miny<=${base_maxy} \
        )\
    ), base WHERE ST_INTERSECTS(geom_inrange,base.geom)" \
example.gpkg > /dev/null

I am going to call this script above, and ogrinfo -sql "SELECT circles.geom FROM circles, base WHERE ST_INTERSECTS(circles.geom, base.geom)" example.gpkg > /dev/null, measuring the time needed for their completion for different number of geometries in circles. I work in a Jupyter-Lab notebook, where lines starting with ! are interpreted as bash calls. I do:

import shapely.geometry
import geopandas as gpd
import numpy as np
import time
np.random.seed(42)

elapsed_with_rtree = []
elapsed_without_rtree = []

base=gpd.GeoSeries(shapely.geometry.Point(0,0).buffer(1,resolution=3))

circle_counts = [2000,5000,10000,15000,20000,40000,80000,200000]

for circle_count in circle_counts:

    circles = gpd.GeoSeries(
        [
            shapely.geometry.Point(e)
            for e in np.random.uniform(low=-10, high=10, size=[circle_count, 2])
        ]
    ).buffer(distance=0.1)

    !rm -f example.gpkg
    base.to_file("example.gpkg",layer="base")
    circles.to_file("example.gpkg",layer="circles")
    
    t0 = time.time()
    !ogrinfo -sql "SELECT circles.geom FROM circles, base WHERE ST_INTERSECTS(circles.geom, base.geom)" example.gpkg > /dev/null    
    t1 = time.time()
    !bash query.sh
    t2 = time.time()
    
    elapsed_without_rtree.append(t1-t0)
    elapsed_with_rtree.append(t2-t1)

I plot the results:

import matplotlib.pyplot as plt

plt.figure(figsize=(8,4))

plt.scatter(circle_counts,elapsed_without_rtree,c='k')
plt.plot(circle_counts,elapsed_without_rtree,c='r',label="without rtree")

plt.scatter(circle_counts,elapsed_with_rtree,c='k')
plt.plot(circle_counts,elapsed_with_rtree,c='b',label="with rtree")

plt.xlabel("Number of intersections checked")
plt.ylabel("Time [sec]")
plt.legend()

It is clear: when circles have a large number of geometries, the solution utilising manual queries in rtree_circles_geom performs better.

---

Fast execution without manual rtree processing

I am looking for a way to achieve the similarly fast runtime without manually having to write up the WHERE clause on rtree_circles_geom. Something like:

ogr2ogr -sql "SELECT circles.geom FROM circles, base WHERE ST_INTERSECTS(circles.geom, base.geom)" \
output.gpkg input.gpkg \
-SOME_OPTION Optimize_for_large_input_Use_rtree_tables

in other, more general words:

Given a GPKG file with two layers, first_layer with a single geometry, second_layer with many of geometries, how can I easily select those geometries from second_layer which intersect the the single geometry in first_layer, in a way which scales for high number of geometries in the second_layer?

Answer 1

Depending on your exact needs there are different ways to do this efficiently.

For the specific case presented here I would use the mask parameter in geopandas.read_file. Under the hood, using GDAL, the spatial index will be used as first pass filter before checking the actual intersection, so it should be fast for small and large datasets being filtered.

Timings

For 2000 circles, took 0.06476800000018557 seconds
For 200000 circles, took 0.3340591000014683 seconds

Script

from pathlib import Path
import time
import geopandas as gpd
import numpy as np
import shapely.geometry

np.random.seed(42)
base=gpd.GeoSeries(shapely.geometry.Point(0,0).buffer(1,resolution=3))

circle_counts = [2000, 200000]
for circle_count in circle_counts:
    circles = gpd.GeoSeries(
        [
            shapely.geometry.Point(e)
            for e in np.random.uniform(low=-10, high=10, size=[circle_count, 2])
        ]
    ).buffer(distance=0.1)

    Path("example.gpkg").unlink(missing_ok=True)
    base.to_file("example.gpkg", layer="base")
    circles.to_file("example.gpkg", layer="circles")

    start = time.perf_counter()
    base = gpd.read_file("example.gpkg", layer="base")
    circles_filtered = gpd.read_file("example.gpkg", layer="circles", mask=base)
    print(f"For {circle_count} circles, took {time.perf_counter() - start} seconds")