pycanopy 0.3.3

Declarative spatial query layer for Polars

A spatial query layer for Polars. Rust core, Python API.

---

[!NOTE] Highly competitive on Apache SpatialBench (single node spatial query benchmark): fastest on 11/24 testcases, within 5% of winning on 14/24 testcases

_{Apache SpatialBench SF1 · lower is better · bars past the cap truncated with their value · TIMEOUT / ERROR annotated}

---

Installation

pip install pycanopy

Pre-built wheels for Linux, macOS, and Windows. No Rust toolchain required.

import polars as pl
from pycanopy import SpatialFrame

sf = SpatialFrame(pl.read_parquet("cities.parquet"), x_col="lon", y_col="lat")
result = sf.lazy().filter(pl.col("population") > 100_000).range_query(-10.0, 35.0, 40.0, 70.0).collect()

---

Why PyCanopy

The driving motivator behind creating this library was to provide the optimizations of relational DBs (query planning, indexing, etc) in a fast, Polars-like interface meant for in-memory spatial work.

Edit [June 19 2026]: Apache SedonaDB released a cool Python DataFrame API. There are similarities between their API and this tool but some key differences are that this library uses (1) a Polars-native query engine and (2) a cost model that decides whether and how to index.

	PyCanopy	GeoPandas	DuckDB	SedonaDB	Spatial Polars
Polars-native API	✓	✗	✗	✗	✓
Spatial query planner (reorder, pushdown, etc)	✓	✗	✓	✓	✗
Index vs scan decided by cost model	✓	✗	✗	✗	✗
Dynamic index selection	✓	✗	✗	✗	✗

---

Example Operations

Inspecting the query plan

lf = (
    sf.lazy()
    .range_query(min_x=-10.0, min_y=35.0, max_x=40.0, max_y=70.0)
    .filter(pl.col("population") > 100_000)
)
print(lf.explain())
# RANGE_QUERY [(-10, 35) → (40, 70)]
# FROM
#   FILTER [(col("population")) > (dyn int: 100000)]
#   FROM
#     DF [N=100,000; path: EXPR]

The optimizer moved the scalar filter below the range query. It runs first on all rows, then the spatial index is probed on the smaller survivor set.

kNN join

query_df = pl.DataFrame({"qx": [2.35, 13.4], "qy": [48.85, 52.5]})

result = sf.lazy().knn_join(query_df, x_col="qx", y_col="qy", k=3).collect()

For each row in query_df, returns the 3 nearest rows in the SpatialFrame. Large probes are streamed in morsels automatically.

Proximity join with aggregation

import pycanopy as pc

stats = (
    sf.lazy()
    .within_distance_join(landmarks, x_col="lon", y_col="lat", distance=0.5)
    .group_by(["landmark"])
    .agg(count=pc.agg.count(), avg_fare=pc.agg.mean("fare"))
)

The full pair frame is never materialised. Each probe morsel folds into per-group partials and combines at the end.

Polygon intersects self-join

from shapely.geometry import box

polygons = [box(i, 0, i + 1.5, 1.0) for i in range(10_000)]
sf = SpatialFrame.from_polygons(pl.DataFrame({"id": range(10_000), "geom": polygons}), geometry_col="geom")

overlaps = sf.intersects_pairs(key_col="id")

Returns all intersecting polygon pairs with overlap area and IoU. key_col replaces positional indices with values from that column.

[!NOTE] For the full operation catalogue, index modes, streaming joins, and API reference see the docs site.

---

Benchmarks

Apache SpatialBench

Run on a single m7i.2xlarge (8 vCPU, 32 GB), the same hardware used by Apache SpatialBench. PyCanopy is measured live with index_mode="auto". Results were produced using the benchmark harness in bench/spatial_bench.

PyCanopy wins a total of 11/24 testcases and lands within 5% of winning 14/24 testcases (there is some variance among benchmark runs).

SF1 (~6M trips)

_{Apache SpatialBench SF1 · lower is better · linear axis, bars past the cap truncated with their value · TIMEOUT / ERROR annotated}

SF10 (~60M trips)

_{Apache SpatialBench SF10 · lower is better · linear axis, bars past the cap truncated with their value · TIMEOUT / ERROR annotated}

All times in seconds. Bold = fastest on that query. SedonaDB, DuckDB, and GeoPandas baselines from published SpatialBench results.

SF1 Query PyCanopy SedonaDB DuckDB GeoPandas q1 1.39 0.66 0.96 12.78 q2 3.74 8.07 9.95 20.74 q3 1.23 0.80 1.17 13.59 q4 7.44 8.41 9.83 25.24 q5 1.71 5.10 1.80 47.08 q6 5.51 8.59 9.36 24.43 q7 2.15 1.66 1.82 137.00 q8 1.04 1.10 1.08 16.08 q9 0.23 0.23 50.15 0.28 q10 8.65 18.79 207.84 46.13 q11 9.90 32.98 TIMEOUT 51.01 q12 14.86 14.55 ERROR TIMEOUT

SF10 Query PyCanopy SedonaDB DuckDB GeoPandas q1 8.52 3.04 4.58 ERROR q2 9.39 8.89 8.26 ERROR q3 6.88 4.09 5.17 TIMEOUT q4 17.34 7.52 8.51 ERROR q5 14.60 50.81 14.40 ERROR q6 11.07 9.11 10.67 ERROR q7 22.73 14.44 14.03 ERROR q8 7.30 7.24 7.57 TIMEOUT q9 0.34 0.38 942.98 0.49 q10 27.26 42.02 ERROR ERROR q11 37.21 97.52 ERROR ERROR q12 175.31 145.66 ERROR TIMEOUT

---

How It Works

The engine has dedicated components for query planning / execution and ultimately returns a Polars DataFrame.

Query flow

flowchart LR
    A[User chain] --> B[SpatialOptimizer] --> C[SpatialExecutor] --> F[pl.DataFrame]

Logical planning

• Predicate pushdown: scalar filters run first, reducing rows before any spatial work.
• Fusion: consecutive range/contains predicates get interleaved into a single Rust call.
• Join side: indexes on the side that makes the join most efficient.
• Projection pushdown: a terminal .select() narrows both join sides before the gather.
• IO path: low-selectivity queries return results as a direct slice, bypassing the Polars expression pipeline.
• EXPR path: runs the spatial engine as a Polars map_batches expression over the query set.

Cost model

index_mode determines how we use the cost model:

Mode	Behaviour
auto (default)	build index when cost model allows it
eager	always build the selected index type, skip the cost check
none	always scan

When index_mode="auto", the planner picks the minimum-cost option ($Q$ queries, $N$ items):

$$ \text{winner} = \arg\min \begin{cases} \text{Cost}{\text{probe}}(\text{built index}) & \text{build already paid} \ \text{Cost}{\text{build}} + \text{Cost}{\text{probe}}(\text{best new index}) \ \text{Cost}{\text{probe}}(\text{brute force}) \end{cases} $$

Selectivity (fraction of the dataset expected to match):

$$ \text{sel} = \begin{cases} \text{hist}(\text{bbox}) / N & \text{range (32×32 density histogram)} \ k / N & \text{kNN} \ 1 / N & \text{contains} \end{cases} $$

Probe cost ($Q$ warm queries against a built index):

$$ \text{Cost}{\text{probe}} = Q \times \begin{cases} N \cdot c{\text{scan}} & \text{brute force} \ (\log_2 N + \text{sel} \cdot N) \cdot c_{\text{tree}} & \text{KD-tree or R-tree} \ \text{sel} \cdot N \cdot c_{\text{grid}} & \text{grid} \end{cases} $$

Build cost (paid once):

$$ \text{Cost}{\text{build}} = \begin{cases} 0 & \text{brute force} \ N \cdot c{\text{build}} & \text{grid} \ N \log_2 N \cdot c_{\text{build}} & \text{KD-tree or R-tree} \end{cases} $$

The empirical constants ($c_{\text{scan}}$, $c_{\text{tree}}$, $c_{\text{grid}}$, $c_{\text{build}}$) are calibrated from benchmark runs in bench/ops.

Index selection

select_index is a rule-based pre-filter that picks a candidate index type:

flowchart TD
    A[Query arrives] --> B{N < 500\nor sel > 50%?}
    B -- yes --> BF[Brute force]
    B -- no --> C{kNN with\nk/N > 10%?}
    C -- yes --> BF
    C -- no --> D{Polygon\ndataset?}
    D -- yes --> RT[R-tree]
    D -- no --> E{Range query\nand uniform?}
    E -- yes --> GR[Grid]
    E -- no --> KD[KD-tree]

All index types share the same coordinate arrays with no duplication.

Why Rust

The hot paths need packed immutable index structures, zero-copy array slices at the Python boundary, and loop-level parallelism. C++ would require a separate FFI layer and would lose the native Polars plugin integration that PyO3/Maturin provides for free.

---

Acknowledgements

Some works that inspired this project:

• Polars: a columnar DataFrame engine that PyCanopy builds on
• geo-index: provides packed, immutable, zero-copy KD-tree and R-tree structures used
• Spatial Polars: an earlier effort to bring spatial functionality to Polars
• Apache Sedona: state-of-the-art spatial SQL engine + benchmark for evals

---

License

MIT