#Apache Sedona for Distributed Spatial Joins

There is a specific shape of problem that single-node spatial tools cannot solve, and it is worth naming precisely, because it is much rarer than the number of Spark clusters in production would suggest. It is not “the dataset is large” — DuckDB reads a terabyte of GeoParquet from object storage without complaint. It is not “there are many files” — Athena scans millions of them. The shape that genuinely needs a cluster is a join whose intermediate candidate set does not fit on one machine: 800 million GPS pings against 12 million building footprints, where the pairs to be tested outnumber either input by orders of magnitude.

Apache Sedona is Spark with spatial types, spatial partitioners, spatial indexes, and a join planner that knows what to do with them. This page belongs to the Query Engines & Cloud Analytics section and covers the part that is actually hard: choosing a join strategy, partitioning both sides so candidate pairs meet on the same executor, and controlling the shuffle when geography refuses to distribute evenly.

#Prerequisites

• Apache Spark 3.5+ with Sedona 1.6+ (sedona-spark-shaded plus geotools-wrapper for CRS support), driven from PySpark on Python 3.10+
• Both inputs as GeoParquet in a single CRS, written with covering bbox statistics; mixing CRS across inputs produces a join that runs fine and returns nothing
• Executor memory sized for the candidate set, not the inputs — the intermediate is what fails, and it is usually the number nobody measured
• Access to the Spark UI or event log, because every diagnosis below is a question about the physical plan or the per-task time distribution

#Why a Spatial Join Is Hard to Distribute

An equi-join distributes trivially: hash both sides on the join key, send equal hashes to the same executor, done. A spatial join has no such key. ST_Intersects(a.geom, b.geom) is a geometric relationship, not an equality, and there is no hash function for which intersecting geometries collide.

Absent help, the engine does the only thing correctness allows: it compares every row on the left against every row on the right. For 800 million by 12 million, that is roughly 10¹⁶ predicate evaluations. No cluster size makes that finish.

Every real strategy is a way of avoiding it, and there are exactly two.

Broadcast. Ship the entire small side to every executor, build an in-memory R-tree over it once per executor, and stream the large side past it. Each large-side row queries a local index; nothing shuffles. This is dramatically the faster strategy when it fits — and it stops working abruptly, with an out-of-memory failure, when the small side outgrows executor memory.

Spatial partitioning. Cut space into cells with a shared partitioner, send both sides’ features to the cells they touch, and join locally within each cell. Features spanning a boundary are replicated into every cell they touch, so a deduplication step follows. This scales to any size and costs a shuffle of both inputs.

#The Workflow

#1. Measure both sides before choosing

The number that decides the strategy is the in-memory size of the small side after decoding, which is typically three to six times its compressed GeoParquet size because WKB expands and JVM objects carry overhead. Compare that against executor memory minus the working set, not against total cluster memory. A 3 GB Parquet file is not a 3 GB broadcast.

#2. Broadcast if it fits, with a real margin

Set the broadcast threshold explicitly rather than trusting the planner’s size estimate, which is derived from file statistics and is frequently wrong for geometry columns. Leave headroom: a broadcast that fits at 90% of executor memory will fail the first time the data grows.

#3. Otherwise partition both sides on one shared grid

Both inputs must use the same partitioner instance — the same tree, the same cell boundaries — or features that intersect will land in different partitions and simply never be compared. This is the most common cause of a spatial join that runs cleanly and returns too few rows.

#4. Let the filter step do its job

Sedona’s join runs a filter-and-refine cycle: the index answers the bounding-box question cheaply, and the exact predicate runs only on candidates that survive. The filter’s selectivity depends entirely on how tight the bounding boxes are, which is decided upstream by space-filling-curve sorting and geometry simplification. A join over sprawling multipart geometries has a weak filter stage no matter what the engine does.

#5. Tune for skew, not for parallelism

Read the per-task duration distribution for the join stage. If the median task takes 8 seconds and the maximum takes 400, the problem is a dense cell, and raising spark.sql.shuffle.partitions will not touch it — the dense cell stays one cell. Switch to a density-aware partitioner so cells subdivide where features concentrate.

#Production Implementation

The function below sizes both inputs, picks a strategy from the measurement rather than from the planner’s estimate, and runs the join with an explicit partitioner when broadcasting is not safe.

# Requires: apache-sedona[spark]>=1.6, pyspark>=3.5  (Python 3.10+)
from __future__ import annotations

from dataclasses import dataclass

from pyspark.sql import DataFrame, SparkSession
from pyspark.sql import functions as F
from sedona.spark import SedonaContext

# Decoded geometry is 3-6x its compressed Parquet size; broadcast only if the
# decoded estimate leaves room for the executor's own working set.
DECODE_FACTOR = 5.0
BROADCAST_SAFETY = 0.35     # use at most this fraction of executor memory


@dataclass(frozen=True)
class JoinPlan:
    strategy: str
    small_side_mb: float
    reason: str


def build_session(app: str, *, executor_memory_gb: int = 16) -> SparkSession:
    """A Sedona-enabled session with the spatial serialisers registered."""
    builder = (
        SparkSession.builder.appName(app)
        .config("spark.serializer", "org.apache.spark.serializer.KryoSerializer")
        .config("spark.kryo.registrator", "org.apache.sedona.core.serde.SedonaKryoRegistrator")
        .config("spark.executor.memory", f"{executor_memory_gb}g")
        .config("spark.sql.adaptive.enabled", "true")
        # Adaptive skew handling splits oversized shuffle partitions; it helps
        # but does not replace a density-aware spatial partitioner.
        .config("spark.sql.adaptive.skewJoin.enabled", "true")
    )
    return SedonaContext.create(builder.getOrCreate())


def plan_join(
    spark: SparkSession, small_path: str, *, executor_memory_gb: int = 16
) -> JoinPlan:
    """Choose broadcast or partitioned from measured size, not from estimates."""
    stats = spark.read.format("geoparquet").load(small_path)
    compressed_mb = (
        spark.sparkContext._jvm.org.apache.hadoop.fs.FileSystem
        .get(spark.sparkContext._jsc.hadoopConfiguration())
        .getContentSummary(
            spark.sparkContext._jvm.org.apache.hadoop.fs.Path(small_path)
        ).getLength() / (1 << 20)
    )
    decoded_mb = compressed_mb * DECODE_FACTOR
    budget_mb = executor_memory_gb * 1024 * BROADCAST_SAFETY

    if decoded_mb <= budget_mb:
        return JoinPlan(
            "broadcast", decoded_mb,
            f"decoded ~{decoded_mb:.0f} MB fits the {budget_mb:.0f} MB broadcast budget",
        )
    return JoinPlan(
        "partitioned", decoded_mb,
        f"decoded ~{decoded_mb:.0f} MB exceeds the {budget_mb:.0f} MB budget",
    )


def spatial_join(
    spark: SparkSession,
    large_path: str,
    small_path: str,
    *,
    partitions: int = 512,
    executor_memory_gb: int = 16,
) -> tuple[DataFrame, JoinPlan]:
    """Join two GeoParquet datasets on ST_Intersects with an explicit strategy."""
    plan = plan_join(spark, small_path, executor_memory_gb=executor_memory_gb)

    large = spark.read.format("geoparquet").load(large_path).alias("l")
    small = spark.read.format("geoparquet").load(small_path).alias("s")

    if plan.strategy == "broadcast":
        small = F.broadcast(small)
    else:
        # KDB-tree partitioning subdivides by density, so a cell over a city
        # holds a comparable feature count to a cell over farmland.
        spark.conf.set("sedona.join.gridtype", "kdbtree")
        spark.conf.set("sedona.join.numpartition", str(partitions))
        spark.conf.set("spark.sql.shuffle.partitions", str(partitions))

    joined = large.join(
        small,
        F.expr("ST_Intersects(l.geometry, s.geometry)"),
        how="inner",
    )
    return joined, plan

Before trusting a run, confirm the planner actually engaged a spatial strategy — this check catches the nested-loop fallback that turns an eight-minute job into an eight-hour one:

# Requires: pyspark>=3.5 — assert the physical plan is not quadratic
def assert_spatial_strategy(df) -> str:
    """Fail loudly if the planner fell back to a nested-loop join."""
    plan = df._jdf.queryExecution().executedPlan().toString()
    if "BroadcastIndexJoin" in plan:
        return "broadcast-index"
    if "RangeJoin" in plan or "DistanceJoin" in plan:
        return "spatial-partitioned"
    raise RuntimeError(
        "no spatial join strategy in the plan — this will run as a nested loop.\n"
        + plan[:600]
    )

#Reference Matrix

Measured on 812 million GPS pings joined against 11.6 million building footprints, both GeoParquet in EPSG:3857, on a 12-executor cluster with 16 GB executors.

The last row is the one worth internalising. Adding a bounding-box predicate on the large side before the join — restricting to the area of interest — reduced the work below what any partitioning strategy achieved, because it eliminated rows rather than distributing them. Filtering before joining beats every join optimisation, and it depends on the layout work described in partitioning GeoParquet for cost control.

#Failure Modes and Gotchas

The silent nested loop. The query runs. It produces correct results. It takes eleven hours. Nothing in the logs says “I gave up on your index”, because from the engine’s perspective nothing went wrong. Assert on the physical plan in the job itself, as the snippet above does, so the fallback is a failure rather than a surprise on the invoice.

Different partitioners on the two sides. Two grids with different cell boundaries put intersecting features in different partitions, so they are never compared. The join completes quickly and returns a plausible-looking but incomplete result — the worst possible failure, because nothing is obviously broken. Share one partitioner instance across both inputs.

Forgetting to deduplicate after a partitioned join. A polygon straddling four cells is replicated into all four, so a point near the boundary can match it more than once. Sedona handles this for its own join operators, but a hand-rolled partition-then-join pipeline must deduplicate on the feature key pair explicitly.

Mixed CRS across inputs. Spatial predicates treat coordinates as plain numbers, so joining a layer in EPSG:4326 with one in EPSG:3857 produces zero matches — not an error. Normalise both sides upstream, the way CRS metadata preservation requires, and assert the CRS in the job.

Reaching for a cluster when a single node would do. Sedona’s operational cost is real: a cluster to size, a Spark version to track, executors to tune, and a much longer feedback loop than a single-node query. If the workload is a scan with a filter, or a join where one side is genuinely small, DuckDB on one large VM will finish sooner and cost less. Reserve the cluster for the joins that need it.

#Frequently Asked Questions

When is Sedona the right tool instead of DuckDB or Athena?

When the join itself, not the scan, is what exceeds one machine. DuckDB handles remarkably large spatial workloads on a single large VM, and Athena is excellent for scan-and-filter queries over partitioned GeoParquet. Sedona earns its operational cost when you are joining two large geometry collections — hundreds of millions of points against millions of polygons — where the intermediate candidate set will not fit in one node’s memory no matter how the scan is optimised.

Why is a spatial join slow even though both inputs are small?

Almost always because the planner fell back to a nested-loop join. A spatial predicate is not an equality, so the engine cannot hash on it; unless a spatial index or a broadcast strategy is engaged, it compares every row against every row. Check the physical plan for a broadcast or index join node. If you see a Cartesian or nested-loop node, the fix is to make the small side broadcastable or to partition both sides — not to add executors, which only makes the same quadratic work parallel.

What causes skew in a spatial join and how do you fix it?

Geographic data is never uniform: a grid cell over a capital city can hold a thousand times the features of a rural cell, so one task processes a thousand times the work and the whole stage waits for it. Raising parallelism does not help, because the partitioning is what is uneven. The fix is a density-aware partitioner — a quadtree or KDB-tree that subdivides dense regions further — so every partition holds a comparable number of features regardless of area.

Does Sedona read GeoParquet natively?

Yes. Sedona reads GeoParquet through Spark’s Parquet reader and decodes the geometry column without a conversion step, which means partition pruning and column projection work the way they do for any Parquet table. Files written with covering bbox statistics and sorted on a space-filling curve give the reader the row-group skipping it needs, so the layout work that pays off for DuckDB and Athena pays off here too.

#Related

• Broadcast vs Partitioned Spatial Joins in Sedona — choosing between the two strategies from measured sizes
• Tuning Spark Shuffle for Large Spatial Joins — partition counts, skew handling, and spill control
• DuckDB Spatial Extension for GeoParquet — the single-node alternative that handles more than most teams expect
• Space-Filling Curves for Spatial Partitioning — the upstream layout that makes the filter stage selective
• Spatial Partitioning with Quadtree Indexes — the density-aware structure used to defeat skew

← Back to Query Engines & Cloud Analytics for GeoParquet