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.readparquet("cities.parquet"), xcol="lon", y_col="lat") result = sf.lazy().filter(pl.col("population") > 100000).rangequery(-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()
    .rangequery(minx=-10.0, miny=35.0, maxx=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().knnjoin(querydf, xcol="qx", ycol="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()     .withindistancejoin(landmarks, xcol="lon", ycol="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.frompolygons(pl.DataFrame({"id": range(10000), "geom": polygons}), geometry_col="geom")
overlaps = sf.intersectspairs(keycol="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 indexmode="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.

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} \\ (\log2 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 \log2 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