Left joining CSV to GeoPackage file and outputting GeoJSON in scalable way with ogr2ogr

Question

I have a large GPKG and a CSV file, and I'd like to (left) join the CSV to the GPKG on one of the columns. A visualisation of what I'm trying to achieve:

Reproducible example

We can create an example GPKG file and an example CSV by using this Python code:

import shapely.geometry
import geopandas as gpd
import pandas as pd

p0 = shapely.geometry.Point([0, 0])

rowCount = 4

gdf = gpd.GeoDataFrame.from_dict(
    dict(
        geometry=[p0 for _ in range(rowCount)],
        col0=[f"val{i}" for i in range(rowCount)],
    )
)
df = pd.DataFrame.from_dict(
    dict(
        col0=[f"val{i}" for i in range(rowCount)],
        col1=[f"some_other_value{i}" for i in range(rowCount)],
    )
).sample(frac=1)

gdf.to_file("base.gpkg",layer="base")
df.to_csv("join.csv")

This code will create a GPKG and a CSV file (for contents, see visualisation above).

We can perform the join using the INDIRECT_SQLITE dialect (why not just stick to simple sqlite? see this thread, or the example under the SELECT statement section of these GDAL docs):

ogr2ogr output.geojson base.gpkg \
   -dialect INDIRECT_SQLITE \
    -sql "SELECT \
            base.geom AS Geometry, \
            base.col0 AS col0, \
            csv.col1 AS col1 \
          FROM \
            base \
          LEFT JOIN \
            'join.csv'.join AS csv \
          ON \
            base.col0 = csv.col0" \
    -nln output

The output is this GeoJSON, as expected:

{
"type": "FeatureCollection",
"name": "output",
"features": [
{ "type": "Feature", "properties": { "col0": "val0", "col1": "some_other_value0" }, "geometry": { "type": "Point", "coordinates": [ 0.0, 0.0 ] } },
{ "type": "Feature", "properties": { "col0": "val1", "col1": "some_other_value1" }, "geometry": { "type": "Point", "coordinates": [ 0.0, 0.0 ] } },
{ "type": "Feature", "properties": { "col0": "val2", "col1": "some_other_value2" }, "geometry": { "type": "Point", "coordinates": [ 0.0, 0.0 ] } },
{ "type": "Feature", "properties": { "col0": "val3", "col1": "some_other_value3" }, "geometry": { "type": "Point", "coordinates": [ 0.0, 0.0 ] } }
]
}

Problem

The problem with using INDIRECT_SQLITE the way I do is that runtime seems to scale with the square of the number of rows in the GPKG & CSV files. Runtime visualisation:

I produced the data for this graph in a Jupyter notebook (meaning: lines starting with ! are shell commands rather than Python), this way:

import shapely.geometry
import geopandas as gpd
import pandas as pd
import time

p0 = shapely.geometry.Point([0, 0])
times = {}

for rowCount in [100,400,1000,2000,3000,4000,5000,6000,7000]:

    print(f"rowCount={rowCount}")

    gdf = gpd.GeoDataFrame.from_dict(
        dict(
            geometry=[p0 for _ in range(rowCount)],
            col0=[f"val{i}" for i in range(rowCount)],
        )
    )
    df = pd.DataFrame.from_dict(
        dict(
            col0=[f"val{i}" for i in range(rowCount)],
            col1=[f"some_other_value{i}" for i in range(rowCount)],
        )
    ).sample(frac=1)

    !rm -f base.gpkg
    !rm -f join.csv
    !rm -f output.geojson
    
    gdf.to_file("base.gpkg",layer="base")
    df.to_csv("join.csv")

    t0 = time.time()

    !\
    ogr2ogr output.geojson base.gpkg \
    -dialect INDIRECT_SQLITE \
        -sql "SELECT \
                base.geom AS Geometry, \
                base.col0 AS col0, \
                csv.col1 AS col1 \
            FROM \
                base \
            LEFT JOIN \
                'join.csv'.join AS csv \
            ON \
                base.col0 = csv.col0" \
        -nln output

    times[rowCount]=time.time()-t0

Code for the plot itself:

import matplotlib.pyplot as plt

plt.plot(times.keys(),times.values(),c="blue")
plt.scatter(times.keys(),times.values(),c="red")
plt.xlabel("rowCount")
plt.ylabel("Runtime (sec)")

Question

How can I use ogr2ogr to left join a CSV to a GPKG file and produce a GeoJSON output in a scalable way?

By "scalable" I mean: the expected runtime for a large number of lines to join grows linearly.

Answer 1

One possible solution is to load the CSV file to the GPKG file as a new layer. Once it's in the GPKG, there is no need for INDIRECT_SQLITE! So:

ogr2ogr -update base.gpkg join.csv -nln csv

ogr2ogr base.geojson base.gpkg \
    -sql "SELECT \
           base.geom AS Geometry, \
           base.col0 AS col0, \
           csv.col1 AS col1 \
         FROM base \
         LEFT JOIN csv \
         ON base.col0 = csv.col0" \
    -nln output

Scaling (includes both the CSV -> GPKG & GPKG -> GeoJSON steps):

---

We can also create an index on base.col0 and csv.col0 (source) before the JOIN operation:

ogr2ogr -update base.gpkg join.csv -nln csv

ogrinfo -sql "CREATE INDEX base_col0_idx ON base (col0)" base.gpkg
ogrinfo -sql "CREATE INDEX csv_col0_idx ON csv (col0)" base.gpkg

ogr2ogr base.geojson base.gpkg \
     -sql "SELECT \
            base.geom AS Geometry, \
            base.col0 AS col0, \
            csv.col1 AS col1 \
        FROM base \
        LEFT JOIN csv \
        ON base.col0 = csv.col0" \
     -nln output

Indexing seems to help a bit when the number of geometries becomes large:

and more when the number of geometries is even larger: