Fast subtraction of geometries in GPKG file

I have a GPKG file, with 2 layers, p and q. I would like to create a third layer, r. I would like r to contain those parts of p which do not intersect with q. We could phrase this as r=p-q.

Below is a way to achieve this, using `ogr2ogr`'s `-sql` option. I would like to learn how to make the presented method faster.

Simple example

Let p be a chessboard-like grid of squares, and q be circles distributed randomly. I generate such GeoSeries using Python:

``````import geopandas as gpd
import shapely.geometry
import numpy as np

np.random.seed(42)

n = 8 # <- will try multiple different ns below

p = gpd.GeoSeries(
[shapely.geometry.box(j, i, j + 1, i + 1) for i in range(n) for j in range(n)]
)

q = gpd.GeoSeries(
[
shapely.geometry.Point(e)
for e in np.random.uniform(low=0, high=n, size=[n * n * 10, 2])
]
).buffer(distance=0.1)
``````

Visualization of the result:

``````import matplotlib.pyplot as plt
np.random.seed(42)

fig, ax = plt.subplots()

p.plot(ax=ax, edgecolor="red",facecolor="none")
q.plot(ax=ax, edgecolor="black",facecolor="none")

ax.axis("off");
``````

Write GeoSeries to GPKG:

``````!rm -f united.gpkg
p.to_file("united.gpkg",layer="p")
q.to_file("united.gpkg",layer="q")
``````

(I work in a Jupyter-Lab notebook where lines starting with `!` are interpreted as bash calls.)

I use `ogr2ogr` to compute r:

``````ogr2ogr \
-update \
-sql "SELECT ST_Difference(p.geom, (SELECT ST_UNION(geom) from q WHERE ST_INTERSECTS(p.geom,geom))) AS geom FROM p" \
-nln r \
-nlt "MULTIPOLYGON" \
united.gpkg united.gpkg
``````

Inspect results in QGIS:

The query worked as expected.

Benchmarking

Above, the r=p-q operation was performed on relatively few p & q geometries. I would like to time how fast r is calculated for different numbers of geometries within p & q. In other words: how fast does the SQL SELECT statement run, if I vary `n`? (There are `n^2` geometries in layer p.)

I specifically care about how fast the `SQL` query is run, not how fast the GPKG file is written. For this reason, I am going to switch to `ogrinfo` (and I'll suppress the non-error outputs by using `> /dev/null`, as described here).

Execution time measuring script:

``````import geopandas as gpd
import shapely.geometry
import numpy as np
import time

np.random.seed(42)

elapsedList = []
ns = [8,16,24,32,40,48,56,64]

for n in ns:

p = gpd.GeoSeries(
[shapely.geometry.box(j, i, j + 1, i + 1) for i in range(n) for j in range(n)]
)

q = gpd.GeoSeries(
[
shapely.geometry.Point(e)
for e in np.random.uniform(low=0, high=n, size=[n * n * 10, 2])
]
).buffer(distance=0.1)

!rm -f united.gpkg
p.to_file("united.gpkg",layer="p")
q.to_file("united.gpkg",layer="q")

t0 = time.time()
!ogrinfo -q -sql "SELECT ST_Difference(p.geom, (SELECT ST_UNION(geom) from q WHERE ST_INTERSECTS(p.geom,geom))) AS geom FROM p" united.gpkg > /dev/null
elapsedList.append(time.time()-t0)
``````

Visualize execution time:

``````import matplotlib.pyplot as plt

plt.figure(figsize=(8,4))
plt.scatter([e**2 for e in ns],elapsedList)
plt.plot([e**2 for e in ns],elapsedList,c='r')
plt.xlabel("Number of geometries in p")
plt.ylabel("ogrinfo execution time")
``````

It seems that the total time needed for the `SELECT` statement scales worse than linearly with the number of geometries in p.

Fast execution in the large p, q case

My aim is to make the above process faster for large p and q.

A possible approach would be: instead of straight running `ogrinfo` with the `SELECT ST_Difference ... geom FROM p` command as above, one could use `ogrinfo` to create spatial indexes in the GPKG file, so that the `SELECT ST_Difference ... geom FROM p` executes faster, when p and q are large. I don't know much about spatial indexes, but based on `ogrinfo -sql "SELECT HasSpatialIndex('q','geom')" united.gpkg` and `ogrinfo -sql "SELECT HasSpatialIndex('p','geom')" united.gpkg`, it seems these layers already have a spatial index. Output for both:

``````INFO: Open of `united.gpkg'
using driver `GPKG' successful.

Layer name: SELECT
Geometry: Unknown (any)
Feature Count: 1
Layer SRS WKT:
(unknown)
HasSpatialIndex: Integer (0.0)
OGRFeature(SELECT):0
HasSpatialIndex (Integer) = 1
``````

Another approach would be to improve the query itself, ie

``````SELECT ST_Difference(p.geom, (SELECT ST_UNION(geom) FROM q WHERE ST_INTERSECTS(p.geom,geom))) AS geom FROM p
``````

I don't yet know the exact way.

How to make the above r=p-q procedure faster for large p and q?

• Hi. ogrinfo does not create indexes. Also, spatial indexes are created by default when write the geopackage. Maybe you wanted to ask how to apply a spatial filter (based on the spatial index) to the select statement? Dec 4, 2023 at 2:33
• Of course ogrinfo creates indexes if user wants so. It is just SQL and usage is like `ogrinfo my.gpkg -sql "create index my_index on my_table (my_column)"`. Dec 4, 2023 at 6:54
• "patial indexes are created by default when write the geopackage" - that's useful to know, thanks. Does this mean I there are no indexes which could be added to significantly improve performance? Dec 4, 2023 at 9:21
• "Maybe you wanted to ask how to apply a spatial filter (based on the spatial index) to the select statement?" - to rephrase this: should I ask about how to improve the query `SELECT ST_Difference(p.geom, (SELECT ST_UNION(geom) from q WHERE ST_INTERSECTS(p.geom,geom))) AS geom FROM p`? Could you give me a pointer on what should I look up to understand how that query is not nearly as good as it could be @GabrielDeLuca? Dec 4, 2023 at 9:25
• @user30184, thank you for that info, i thought ogrinfo could not write the index to the database, but ` SELECT HasSpatialIndex('table_name','geom_col_name')` can be useful. Dec 4, 2023 at 11:29

The following sql query uses the spatial index:

``````SELECT ST_Difference(
layer1.geom,
(SELECT ST_UNION(layer2.geom)
FROM q layer2
JOIN rtree_q_geom layer2tree ON layer2.rowid = layer2tree.id
WHERE ST_MinX(layer1.geom) <= layer2tree.maxx
AND ST_MaxX(layer1.geom) >= layer2tree.minx
AND ST_MinY(layer1.geom) <= layer2tree.maxy
AND ST_MaxY(layer1.geom) >= layer2tree.miny
AND ST_INTERSECTS(layer1.geom, layer2.geom) = 1
)
) AS geom
FROM p layer1
``````

However, you should be aware that ST_difference(geom, NULL) returns NULL. So if a rectangle doesn't intersect with any circle, NULL will be returned for it! So, unless you are really sure this situation will never occur, you'd rather want to use something like this to get correct results. Performance is the same, at least for this case.

``````SELECT (SELECT IIF(ST_UNION(layer2.geom) IS NULL,
layer1.geom,
ST_Difference(layer1.geom, ST_UNION(layer2.geom))
)
FROM q layer2
JOIN rtree_q_geom layer2tree ON layer2.rowid = layer2tree.id
WHERE ST_MinX(layer1.geom) <= layer2tree.maxx
AND ST_MaxX(layer1.geom) >= layer2tree.minx
AND ST_MinY(layer1.geom) <= layer2tree.maxy
AND ST_MaxY(layer1.geom) >= layer2tree.miny
AND ST_INTERSECTS(layer1.geom, layer2.geom) = 1
) AS geom
FROM p layer1
``````

Using the spatial index speeds up the processing as expected. Timing on on windows, under wsl2 (under native windows, ~2.5x slower):

As always, there are many roads to Rome.

Par example, there are different ways to write your SQL query to use the spatial index. In my experience joining with it gives the best results. According to the following blog post it also scales the best for large numbers of rows: dealing-with-huge-vector-geopackage.

Another alternative is using geofileops. Disclaimer: I'm the developer. Especially if your question is more general in nature and you would also need to process larger files with more complex polygons than this case, geofileops.erase might be a good fit. Because it uses multiprocessing and some other optimizations under the hood, it can speed up processing significantly. For this case however, with relatively few and simple polygons, there are only minor gains in the test cases with larger numbers of polygons,...

``````import geofileops as gfo

if __name__ == "__main__":
# Because geofileops uses multiprocessing, only call it from a function or check for "__main__"
gfo.erase(
input_path="united.gpkg",
erase_path="united.gpkg",
output_path="output.gpkg",
input_layer="p",
erase_layer="q",
subdivide_coords=-1,
)
``````

The following graph compares the performance of 2 variants of sql statements as well as geofileops.erase for this test case. Note: as geofileops.erase always writes the output file, this is included in its timing:

The following script runs the performance tests. I used the pyogrio I/O engine in `geopandas.read_file` to be able to run the sql statements easily in plain python. Using the `pyogrio` engine also in `to_file` speeds up writing files a lot.

``````from pathlib import Path
import geofileops as gfo
import geopandas as gpd
import matplotlib.pyplot as plt
import shapely.geometry
import numpy as np
import time

def benchmark():
np.random.seed(42)

elapsed_original = []
elapsed_rtree_in = []
elapsed_rtree_join = []
elapsed_gfo = []

ns = [8, 16, 24, 32, 40, 48, 56, 64]
for n in ns:
path = Path(f"united_{n}.gpkg")
if not path.exists():
p = gpd.GeoSeries(
[shapely.geometry.box(j, i, j + 1, i + 1) for i in range(n) for j in range(n)]
)

q = gpd.GeoSeries(
[
shapely.geometry.Point(e)
for e in np.random.uniform(low=0, high=n, size=[n * n * 10, 2])
]
).buffer(distance=0.1)

p.to_file(path, layer="p", engine="pyogrio")
q.to_file(path, layer="q", engine="pyogrio")

t0 = time.time()
sql = """
SELECT ST_Difference(
p.geom,
(SELECT ST_UNION(geom)
FROM q
WHERE ST_Intersects(p.geom,geom)
)
) AS geom
FROM p;
"""
elapsed_original.append(time.time() - t0)

t0 = time.time()
sql = """
SELECT (SELECT IIF( ST_UNION(layer2.geom) IS NULL,
layer1.geom,
ST_Difference(layer1.geom, ST_UNION(layer2.geom))
)
FROM q layer2
JOIN rtree_q_geom layer2tree ON layer2.rowid = layer2tree.id
WHERE ST_MinX(layer1.geom) <= layer2tree.maxx
AND ST_MaxX(layer1.geom) >= layer2tree.minx
AND ST_MinY(layer1.geom) <= layer2tree.maxy
AND ST_MaxY(layer1.geom) >= layer2tree.miny
AND ST_INTERSECTS(layer1.geom, layer2.geom) = 1
) AS geom
FROM p layer1
"""
elapsed_rtree_join.append(time.time() - t0)

t0 = time.time()
sql = """
SELECT (SELECT IIF( ST_UNION(q.geom) IS NULL,
p.geom,
ST_Difference(p.geom, ST_UNION(q.geom))
)
FROM q
WHERE ROWID IN(
SELECT id
FROM rtree_q_geom
WHERE minx <= MbrMaxX(p.geom)
AND maxx >= MbrMinX(p.geom)
AND miny <= MbrMaxY(p.geom)
AND maxy >= MbrMinY(p.geom)
)
AND ST_Intersects(p.geom, geom)
) AS geom
FROM p;
"""
elapsed_rtree_in.append(time.time() - t0)

t0 = time.time()
# Note: use subdivide_coords=-1, otherwise an error is raised that both layers
# are in one input file: https://github.com/geofileops/geofileops/issues/451
result_path = "result_erase.gpkg"
gfo.erase(input_path=path, erase_path=path, output_path=result_path, input_layer="p", erase_layer="q", subdivide_coords=-1, force=True)
elapsed_gfo.append(time.time() - t0)

# Plots
plot_timings(ns, elapsed_original, elapsed_rtree_join, [], [])
plot_timings(ns, [], elapsed_rtree_join, elapsed_rtree_in, elapsed_gfo)

def plot_timings(ns, elapsed_original, elapsed_rtree_join, elapsed_rtree_in, elapsed_gfo):
# Print all passed results
plt.figure(figsize=(8, 4))
if len(elapsed_original) > 0:
print(f"elapsed_original: {elapsed_original}")
plt.scatter([e**2 for e in ns], elapsed_original, c="r")
plt.plot([e**2 for e in ns], elapsed_original, c="r", label="original")
if len(elapsed_rtree_join) > 0:
print(f"elapsed_rtree_join: {elapsed_rtree_join}")
plt.scatter([e**2 for e in ns], elapsed_rtree_join, c='g')
plt.plot([e**2 for e in ns], elapsed_rtree_join, c='g', label="rtree, join")
if len(elapsed_rtree_in) > 0:
print(f"elapsed_rtree_in: {elapsed_rtree_in}")
plt.scatter([e**2 for e in ns], elapsed_rtree_in, c='b')
plt.plot([e**2 for e in ns], elapsed_rtree_in, c='b', label="rtree, in")
if len(elapsed_gfo) > 0:
print(f"elapsed_gfo: {elapsed_gfo}")
plt.scatter([e**2 for e in ns], elapsed_gfo, c='m')
plt.plot([e**2 for e in ns], elapsed_gfo, c='m', label="geofileops.erase")
plt.xlabel("Number of geometries in p")
plt.ylabel("execution time")
plt.legend()
plt.show()

if __name__ == "__main__":
benchmark()
``````

In PostGIS, some predicates like ST_Intersects make use of the spatial index automatically (see the Note at https://postgis.net/docs/ST_Intersects.html).

But the ST_Intersects name in SpatiaLite is just an alias of its Intersects predicate, which doesn't make use of the spatial index. We need to explicitly include in our queries a spatial filter to quickly get rid of geometries whose bounding box rectangles do not intersect.

``````SELECT ST_Difference(
p.geom,
(SELECT ST_UNION(geom)
FROM q
WHERE ROWID IN(
SELECT id
FROM rtree_q_geom
WHERE minx <= MbrMaxX(p.geom)
AND maxx >= MbrMinX(p.geom)
AND miny <= MbrMaxY(p.geom)
AND maxy >= MbrMinY(p.geom)
)
AND ST_Intersects(p.geom, geom)
)
) AS geom FROM p;
``````

As Pieter answer, there are many ways to query the spatial index.

Spatial indexes are trees where each node represent a rectangle that belongs within its parent node. A brief, outdated but useful introduction was described in: http://www.gaia-gis.it/gaia-sins/spatialite-cookbook/html/rtree.html

When a GeoPackage spatial table is written by GDAL, its spatial index is created by default, but can be avoided through a layer creation option: https://gdal.org/drivers/vector/gpkg.html#layer-creation-options.

SpatiaLite implements a new way to query the spatial index, through the SpatialIndex table. It is a virtual table wrapping the rectangles tree, and it is presented in: https://gaia-gis.it/fossil/libspatialite/wiki?name=SpatialIndex.
But neither the GeoPackage standard (nor its GDAL implementation) includes it. Instead, the spatial index is wrapped in SQLite's own rtree virtual tables, which are documented in: https://www.sqlite.org/rtree.html.

Regarding the ROWID attribute, it is documented in https://www.sqlite.org/lang_createtable.html#rowid. In the GeoPackage standard, the Feature ID attribute refers to it.

The complete code I tested is the following:

``````import geopandas as gpd
import shapely.geometry
import numpy as np
import time

np.random.seed(42)

elapsed_with_rtree = []
elapsed_without_rtree = []

ns = [8, 16, 24, 32, 40, 48, 56, 64]

for n in ns:

p = gpd.GeoSeries(
[shapely.geometry.box(j, i, j + 1, i + 1) for i in range(n) for j in range(n)]
)

q = gpd.GeoSeries(
[
shapely.geometry.Point(e)
for e in np.random.uniform(low=0, high=n, size=[n * n * 10, 2])
]
).buffer(distance=0.1)

!rm -f united.gpkg
p.to_file("united.gpkg", layer="p")
q.to_file("united.gpkg", layer="q")

query = """
SELECT ST_Difference(
p.geom,
(SELECT ST_UNION(geom)
FROM q
WHERE ST_Intersects(p.geom,geom)
)
) AS geom FROM p;
"""
query = ' '.join(query.split())

rtree_query = """
SELECT ST_Difference(
p.geom,
(SELECT ST_UNION(geom)
FROM q
WHERE ROWID IN(
SELECT id
FROM rtree_q_geom
WHERE minx <= MbrMaxX(p.geom)
AND maxx >= MbrMinX(p.geom)
AND miny <= MbrMaxY(p.geom)
AND maxy >= MbrMinY(p.geom)
)
AND ST_Intersects(p.geom, geom)
)
) AS geom FROM p;
"""
rtree_query = ' '.join(rtree_query.split())

t0 = time.time()
subprocess.run(['ogrinfo', '-q', '-sql', query, 'united.gpkg'], stdout=subprocess.DEVNULL)
t1 = time.time()
subprocess.run(['ogrinfo', '-q', '-sql', rtree_query, 'united.gpkg'], stdout=subprocess.DEVNULL)
t2 = time.time()

elapsed_without_rtree.append(t1-t0)
elapsed_with_rtree.append(t2-t1)

import matplotlib.pyplot as plt

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

plt.scatter([e**2 for e in ns], elapsed_without_rtree, c='k')
plt.plot([e**2 for e in ns], elapsed_without_rtree, c='r', label="without rtree")

plt.scatter([e**2 for e in ns], elapsed_with_rtree, c='k')
plt.plot([e**2 for e in ns], elapsed_with_rtree, c='b', label="with rtree")

plt.xlabel("Number of geometries in p")
plt.ylabel("ogrinfo execution time")

plt.legend()
``````

• This looks very good. The error I get, when I run your query: `no such table: SpatialIndex`. So far I tried fixing this by creating a gpkg where I explicitly create a spatial index: `ogr2ogr united_w_spatial_index.gpkg united.gpkg -lco SPATIAL_INDEX=YES`. This didn't help, and I realized the original gpkg should also have a spatial index, as confirmed by `ogrinfo -sql "SELECT HasSpatialIndex('p','geom')" united.gpkg` (see output in question). Dec 5, 2023 at 9:25