Fast subtraction of geometries in GPKG file

Question

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?

Answer 1

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;
        """
        _ = gpd.read_file(path, sql=sql, engine="pyogrio")
        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
        """
        _ = gpd.read_file(path, sql=sql, engine="pyogrio")
        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;
        """
        _ = gpd.read_file(path, sql=sql, engine="pyogrio")
        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()

Answer 2

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()