# Fast selection of intersecting geometries in GPKG file

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?

• GDAL knows how to use the RTree in simple use case that is a selection by envelope/bbox. Therefore `ogr2ogr -spat <xmin> <ymin> <xmax> <ymax>...` where min/max values present the envelope of your single geometry should be fast. From that filtered selection (bounding boxes intersect) you can continue with ST_Intersects and find out the features which really intersect. Maybe you can find a way for saving the interim resultset into memory for making the query fast. Commented Dec 5, 2023 at 15:57
• Try using DWithin rather than buffers, in postgis it would definitely be faster not sure for sqlite Commented Dec 5, 2023 at 16:25

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)