Add Overture benchmarks to Python benchmarks folder

[paleolimbot]

As referenced in #680, we have some tests that use spatial joins involving OvertureMaps datasets (address, building, landuse). I also have a few that I put together for the release post for 0.3.0 at apache/sedona#2675 (and there are some in our overture docs guide that mostly exercise pruning).

We should keep track of these in the benchmarks directory (even if informally) since coming up with these queries is non-trivial!

[paleolimbot]

Here's one that's related to the blog post I wrote on counting buildings by zipcode a while back. I used census tracts here because that's the first shapefile I found:

# curl -L https://www2.census.gov/geo/tiger/GENZ2024/shp/cb_2024_us_tract_5m.zip -o tracts.shp.zip
# curl -L https://github.com/geoarrow/geoarrow-data/releases/download/v0.2.0/microsoft-buildings_point.parquet -o buildings.parquet
import sedona.db

sd = sedona.db.connect()
sd.options.interactive = True
sd.options.memory_limit = "unlimited"

sd.read_pyogrio("tracts.shp.zip").to_view("tracts")
sd.read_parquet("buildings.parquet").to_view("buildings")

# 6.2s
sd.sql("""
SELECT
    t."GEOIDFQ",
    COUNT(b.geometry) AS building_count
FROM tracts AS t
LEFT JOIN buildings AS b
    ON ST_Contains(t.wkb_geometry, ST_Transform(b.geometry, 'epsg:4269'))
GROUP BY t."GEOIDFQ"
ORDER BY building_count DESC
""").to_memtable().to_view("result", overwrite=True)

I'm not sure what DuckDB is doing differently here but DuckDB (1.5.0 pre) takes 90s to do this:

import duckdb

duckdb.load_extension("spatial")
duckdb.sql("""
CREATE OR REPLACE VIEW tracts AS
SELECT *
FROM ST_Read('tracts.shp.zip')
""")

duckdb.sql("""
SELECT
    t."GEOIDFQ",
    COUNT(b.geometry) AS building_count
FROM tracts AS t
LEFT JOIN "buildings.parquet" AS b
    ON ST_Contains(t.geom, ST_Transform(b.geometry, 'epsg:4326', 'epsg:4269'))
GROUP BY t."GEOIDFQ"
ORDER BY building_count DESC
""").to_table("results")

[paleolimbot]

This one is from the release post: deduplicating buildings from a large ish Parquet file. I think DuckDB is faster at this one because of its hash join but I'm not sure.

# pip install "apache-sedona[db]" --force-reinstall --pre --extra-index-url=https://pypi.fury.io/sedona-nightlies/
# curl -L https://github.com/geoarrow/geoarrow-data/releases/download/v0.2.0/microsoft-buildings_point.parquet -o buildings.parquet
import sedona.db

sd = sedona.db.connect()
sd.options.memory_limit = 'unlimited'

url = "buildings.parquet"
sd.read_parquet(url).to_view("buildings")
sd.sql(
    """
    SELECT ROW_NUMBER() OVER () AS building_id, geometry
    FROM buildings
    """
).to_parquet("buildings_idx.parquet")

sd.sql(
    """
    SELECT
        l.building_id,
        r.building_id AS nearest_building_id,
        l.geometry AS geometry,
        ST_Distance(l.geometry, r.geometry) AS dist
    FROM "buildings_idx.parquet" AS l
    JOIN "buildings_idx.parquet" AS r
        ON l.building_id <> r.building_id
        AND ST_DWithin(l.geometry, r.geometry, 1e-10)
    """
).to_memtable().to_view("duplicates", overwrite=True)

# 30.8s
sd.view("duplicates").count()
#> 1017476

# pip install --force-reinstall --pre duckdb
import duckdb

duckdb.load_extension("spatial")

# 19s
df = duckdb.sql("""
    SELECT
        l.building_id,
        r.building_id AS nearest_building_id,
        l.geometry AS geometry,
        ST_Distance(l.geometry, r.geometry) AS dist
    FROM "buildings_idx.parquet" AS l
    JOIN "buildings_idx.parquet" AS r
        ON l.building_id <> r.building_id
        AND ST_DWithin(l.geometry, r.geometry, 1e-10)
    """)
len(df)
#> 1017476

[paleolimbot]

Here's a slightly smaller one where spilling doesn't really enter into it: find the steepest spots given a mesh of elevation values (i.e., find points between 10 and 100 m away and calculate all the slopes). This one correctly identifies the cliff-like parts of Nova Scotia (and a few quarries / cell towers too).

# pip install "apache-sedona[db]"
import sedona.db

sd = sedona.db.connect()
sd.options.interactive = True

# curl -L https://github.com/geoarrow/geoarrow-data/releases/download/v0.2.0/ns-water_elevation.parquet -o elevation.parquet
sd.read_parquet("elevation.parquet").to_view("elevation")

sd.sql("""
SELECT
    ST_Distance(l.geometry, r.geometry) AS dist,
    abs(ST_Z(l.geometry) - ST_Z(r.geometry)) AS dz,
    l.geometry AS geometry
FROM elevation as l
JOIN elevation as r ON ST_DWithin(l.geometry, r.geometry, 100)
WHERE ST_Distance(l.geometry, r.geometry) > 10
""").to_view("knn_slopes", overwrite=True)

# 3.0s
result = sd.sql("""
SELECT dz / dist AS slope, dist, geometry FROM knn_slopes ORDER BY dz / dist DESC LIMIT 10000
""").to_memtable()
result

DuckDB is slightly faster at this one:

import duckdb

duckdb.load_extension("spatial")

duckdb.sql("""
SELECT
    ST_Distance(l.geometry, r.geometry) AS dist,
    abs(ST_Z(l.geometry) - ST_Z(r.geometry)) AS dz,
    l.geometry AS geometry
FROM "elevation.parquet" as l
JOIN "elevation.parquet" as r ON ST_DWithin(l.geometry, r.geometry, 100)
WHERE ST_Distance(l.geometry, r.geometry) > 10
""").to_view("knn_slopes")

# 2.4s
result = duckdb.sql("""
SELECT dz / dist AS slope, dist, geometry FROM knn_slopes ORDER BY dz / dist DESC LIMIT 10000
""").to_arrow_table()