7

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");

enter image description here

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:

enter image description here

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

enter image description here

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?

8
  • 4
    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
  • 4
    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)".
    – user30184
    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?
    – zabop
    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?
    – zabop
    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

2 Answers 2

7
+50

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

enter image description here

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:

enter image description here

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

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

Figure of time processing with and without rtree query.

6
  • 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).
    – zabop
    Dec 5, 2023 at 9:25
  • 1
    Asked about this issue, in a more generalized way here: gis.stackexchange.com/questions/471410/….
    – zabop
    Dec 5, 2023 at 10:15
  • 1
    Hi. Yes, sorry, my answer is wrong. I answered as if GeoPackage spatial index can be queried the same way as SpatiaLite, but the name of the virtual table is not the same and the search_frame is not implemented. The query changes a bit, let me do some test and will edit my answer. Dec 5, 2023 at 10:48
  • Thanks! I think this story worth a thread on its own. Here it is: gis.stackexchange.com/questions/471430/….
    – zabop
    Dec 5, 2023 at 15:36
  • 1
    For information, based on some tests I did in the past, using the spatial index using a join seemed faster, but not sure if it is the case in all situations. Based on the following (rather old) blog post it also scales better if the number of rows returned are large: dealing-with-huge-vector-geopackage
    – Pieter
    Dec 7, 2023 at 20:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.