Parallel Spatial Joins with Dask-GeoPandas

A spatial join between two out-of-memory datasets is the operation people reach for Dask-GeoPandas to do — and the one that disappoints if you skip the spatial shuffle. This guide runs a points-in-polygons join across partitioned data correctly and in parallel, for analysts joining millions of features that won't fit in a single GeoPandas GeoDataFrame. It sits under Scaling with Dask-GeoPandas in Spatial Analysis & Advanced Query Techniques.

Why This Approach / What Goes Wrong

dask_geopandas.sjoin mirrors the in-memory geopandas.sjoin, and it keeps the exact same join semantics covered in Spatial Joins & Merging — only the distribution changes, not what within or intersects mean. That is precisely where the trap is: parallel correctness depends on how the two inputs are partitioned, something the in-memory API never asks you to think about.

If both inputs are partitioned arbitrarily — typically in row order, so partition 3 holds features scattered across the whole map — the scheduler has no way to know which right partition might match a given left partition. It must test every left partition against every right partition. That all-pairs comparison is O(n²) in partitions and throws away the parallelism you came for; it is routinely slower than just loading the data once and joining on a single machine.

spatial_shuffle is the fix. It repartitions both inputs along a Hilbert space-filling curve so that spatially close features share a partition and each partition gets a compact, non-overlapping bounding box. After the shuffle, a left partition only has to be tested against the handful of right partitions whose boxes actually overlap it; the scheduler prunes everything else. Forgetting the shuffle is the entire difference between a join that scales and one that crawls.

The other classic failure is a CRS mismatch between the two inputs. sjoin evaluates its predicate on raw coordinates, so if the layers disagree on their coordinate reference system every predicate is false — and unlike a broken import, this raises no error. You get an empty or wrong result and no traceback to explain it.

There is a third failure that only appears at scale, and it hits after the join succeeds. A spatial join is not a filter, it is a many-to-many product: each output row is one matched pair. Points inside polygons with within are close to one-to-one, so the output is roughly the size of the input. Polygon-to-polygon with intersects is not — a partition of 400,000 building footprints joined against overlapping zoning districts can emit two or three million rows, and the partition that was sized to fit comfortably in a worker on the way in blows past memory_limit on the way out. The number worth knowing before you launch the full run is the fan-out ratio: output rows divided by left input rows, measured on a single partition. Anything above 1 means your output partitioning is no longer the partitioning you designed, and you should raise npartitions on the left input by that same factor.

Why spatial_shuffle makes a Dask-GeoPandas join scale Two panels compare partition matching. Without spatial_shuffle, features sit in arbitrary row order, so every point partition must be tested against every zone partition — sixteen partition pairs, an all-pairs O(n squared) mesh with no pruning. With spatial_shuffle, both inputs are reordered along a Hilbert curve so partition bounding boxes align by location; each point partition only fans out to the handful of zone partitions that actually overlap it — six pairs here, the rest pruned by the scheduler. Partition alignment decides whether the join scales Without spatial_shuffle With spatial_shuffle (Hilbert) point partitions zone partitions point partitions zone partitions P0 P1 P2 P3 Z0 Z1 Z2 Z3 16 partition pairs tested — O(n²), no pruning P0 P1 P2 P3 Z0 Z1 Z2 Z3 only overlapping pairs — the rest pruned
Why spatial_shuffle makes a Dask-GeoPandas join scale

Prerequisites

conda install -c conda-forge "dask-geopandas=0.4.*" "geopandas=0.14.*" "pyarrow=15.*"

Both inputs should already be stored as partitioned GeoParquet; it is the on-disk format that carries the CRS in metadata and lets the reader skip partitions by bounding box, which keeps the whole pipeline lazy and prunable.

Where a Dask-GeoPandas join stays lazy and where it executes Six stages run left to right. The first five — read_parquet, to_crs, spatial_shuffle, sjoin and groupby.size — sit inside a lazy panel where Dask only records a task graph and no partition is read from disk. The sixth stage, to_parquet, sits outside that panel and is the trigger that executes the whole graph, streaming partitions straight to disk so the driver process never holds the joined result. Badges mark spatial_shuffle as the one full data movement and sjoin as the stage that benefits from partition pruning. The whole join is a graph until something asks for data lazy region — Dask records a task graph; not one partition is read read_parquet lazy handles to_crs per partition spatial_shuffle Hilbert order sjoin matched pairs groupby.size reduction to_parquet writes 128 files spatial_shuffle the one full data movement sjoin only overlapping partitions compared executes the graph streams to disk Only a write or a compute moves data; everything before it is a graph you can still change. Reserve compute for small reductions and send large results straight to Parquet.
Reading, reprojecting, shuffling and joining only build a task graph — the write at the end is what makes the Dask cluster actually touch the data.

Step-by-Step Implementation

1. Read both partitioned inputs and align their CRS.

Reproject before anything else. A Hilbert ordering computed in degrees is not the same ordering you get in metres, and a join across mismatched systems returns silent garbage. Use a proper metric CRS for the region — here ETRS89 / UTM 32N (EPSG:25832) — never Web Mercator (EPSG:3857) for measurement. See Coordinate Systems with PyProj for the CRS conventions.

import dask_geopandas as dgpd

# ~12M GPS pings and ~50k administrative zones, both as partitioned GeoParquet
pings = dgpd.read_parquet("gps_pings/")
zones = dgpd.read_parquet("admin_zones/")

# Joins require a shared CRS — reproject both to the same metric system.
# to_crs runs per-partition and in parallel; it is still lazy here.
pings = pings.to_crs(epsg=25832)
zones = zones.to_crs(epsg=25832)

2. Spatially shuffle so partitions align by location.

Shuffle both sides with a comparable npartitions so their bounding boxes line up. The shuffle is a full data movement, so do it once, here, before the join — not repeatedly.

# Hilbert-curve repartitioning: spatially close features land together.
pings = pings.spatial_shuffle(npartitions=128)
zones = zones.spatial_shuffle(npartitions=128)

The two npartitions values do not have to match, and forcing them to is a common mistake — 50,000 zone polygons split 128 ways gives partitions of 390 features each, which is far below the point where a task is worth scheduling. What has to line up is the geography: both sides ordered on the same curve so their boxes tile the same extent. Size each side from its own byte budget, using the measurement recipe in Partitioning Strategies for Dask-GeoPandas, and let the shuffle handle the alignment.

3. Run the spatial join lazily.

# Assign each ping the zone that contains it.
# predicate="within" tests point-within-polygon; the point layer is the left side.
joined = dgpd.sjoin(pings, zones, how="inner", predicate="within")

Two arguments carry all the semantics here, and both behave the way they do in-memory. predicate is evaluated left relative to right, so within on a point-left/polygon-right call keeps points inside polygons; swap the operands without swapping the predicate and every test fails. how="inner" drops pings that fall in no zone; how="left" keeps them with null zone columns, which is what you want when the unmatched count is itself the answer — how many pings landed outside the coverage area. Be aware of what a left join does to dtypes: the null-filled right-hand columns force an integer zone_id up to float64, so a later groupby("zone_id") silently groups on floats and any downstream write records the wrong type. Cast it back explicitly with .astype("Int64") (the nullable integer) immediately after the join. A right join is not supported by the distributed implementation — reverse the operands and use a left join instead, which computes the same thing.

4. Aggregate and materialize. Stream the result to disk rather than pulling it into the driver process. Reserve .compute() for small reductions; use .to_parquet() whenever the result is itself large.

counts = (
    joined.groupby("zone_id")
    .size()
    .rename("ping_count")
    .reset_index()
)
counts.to_parquet("ping_counts/", write_index=False)

5. (Optional) Use a distributed client for a real cluster.

A single workstation's threaded scheduler handles a surprising amount of GEOS-bound work, because the heavy lifting releases the GIL. Add dask.distributed only when one box runs out of cores or RAM, or when you want the live dashboard.

# from dask.distributed import Client
# client = Client(n_workers=8, threads_per_worker=2, memory_limit="4GB")
# print(client.dashboard_link)   # http://127.0.0.1:8787/status
# ...the same code above now runs across workers...

6. When one side is small, broadcast it instead of shuffling.

The shuffle is the most expensive stage in the pipeline above, and it is avoidable whenever the right-hand layer fits in a worker's memory with room to spare. 50,000 administrative zones is a few hundred megabytes of geometry — send a copy to every worker once, and each left partition joins against the whole zone layer locally. No data movement, no partition alignment to get wrong, and the left side keeps whatever layout it arrived with.

import geopandas as gpd
import dask_geopandas as dgpd
from dask.distributed import Client

client = Client(n_workers=8, threads_per_worker=2, memory_limit="8GB")

pings = dgpd.read_parquet("gps_pings/").to_crs(epsg=25832)   # no shuffle needed
zones = gpd.read_parquet("admin_zones.parquet").to_crs(epsg=25832)

# Ship the small layer to every worker once, not once per task
zones_remote = client.scatter(zones, broadcast=True)


def join_against_zones(part: gpd.GeoDataFrame, zone_layer) -> gpd.GeoDataFrame:
    return gpd.sjoin(part, zone_layer, how="inner", predicate="within")


sample = pings.partitions[0].head(50)
joined = pings.map_partitions(
    join_against_zones, zones_remote,
    meta=join_against_zones(sample, zones).head(0),
)
joined.to_parquet("pings_with_zone/", write_index=False)

The rule of thumb: broadcast when the right-hand layer is under roughly a tenth of one worker's memory_limit, because every worker holds a full copy for the duration of the job on top of its share of the left side. Above that, the shuffle is cheaper than the replication. Note broadcast=True on the scatter — without it, Dask sends the layer to a few workers and the rest pull it over the network on demand, which reintroduces the transfer you were avoiding. Skipping client.scatter entirely and closing over zones in the function is the version that looks simplest and performs worst: the whole layer is then serialised into every single task.

Verification

Confirm the join actually matched features and that the totals reconcile against the input. Checking a reduction is cheap — it never materializes the full joined dataset.

import dask_geopandas as dgpd

result = dgpd.read_parquet("ping_counts/").compute()
print("Zones with pings:", len(result))                        # Zones with pings: 47213
print("Total joined pings:", int(result["ping_count"].sum()))  # Total joined pings: 11984502

# With predicate="within" each ping matches at most one zone, so the joined
# total must not exceed the input ping count.
total_pings = dgpd.read_parquet("gps_pings/").shape[0].compute()
assert result["ping_count"].sum() <= total_pings

A joined total of zero almost always means the two inputs were in different CRSs — go back to step 1 and assert crs.is_projected on a single partition with pings.partitions[0].compute().crs.

Run the cheaper check before the full job, not after it. One partition, joined on its own, tells you the fan-out ratio and the boundary-duplication rate for the price of a few seconds — and both numbers scale linearly, so a single partition predicts the whole run.

import dask_geopandas as dgpd

pings = dgpd.read_parquet("gps_pings/").to_crs(epsg=25832)
zones = dgpd.read_parquet("admin_zones/").to_crs(epsg=25832)

# Join ONE aligned partition pair and read the shape of the result
left = pings.partitions[0].compute()
probe = dgpd.sjoin(
    dgpd.from_geopandas(left, npartitions=1), zones, predicate="within"
).compute()

fan_out = len(probe) / len(left)
dupes = int(probe.index.duplicated().sum())
print(f"fan-out {fan_out:.2f}x, {dupes} pings matched more than one zone")
# fan-out 1.00x, 3 pings matched more than one zone

A fan-out of 1.00 confirms the join really is one-to-one and that output partitions will be the size you planned. The three duplicates are pings that landed exactly on a shared zone boundary — a rounding artefact of coordinates stored at limited precision, not a data error, and the reason the assertion above uses <= rather than ==. If the fan-out comes back above 1, multiply your left-side npartitions by it before the real run so the output partitions land in the target band rather than the input ones.

Reading a failed reconciliation check back to its cause A decision tree starts from the observation that the verification numbers look wrong and splits into three symptoms. A total of zero rows traces to the two inputs sitting in different coordinate systems, fixed by reprojecting both before the shuffle. A total larger than the input row count traces to points on shared zone borders matching two polygons, fixed by dropping duplicates on the point id or applying a tie-break rule. Correct totals with a slow run trace to one or both sides never being spatially shuffled, fixed by shuffling both with a comparable partition count. Three ways the numbers come out wrong — and where each one starts Verification numbers look wrong symptom cause fix joined total is 0 total exceeds input rows totals right, run is slow the two layers sit in different coordinate systems points on a shared border match two zone polygons one or both sides never got the spatial shuffle reproject both to one metric CRS before the shuffle de-duplicate on the ping id or apply a tie-break rule shuffle both sides with a comparable partition count
Each of the three ways this join goes wrong has a single distinguishing symptom, so the reconciliation numbers alone point at the fix.

Edge Cases & Debugging

Frequently Asked Questions

Do I really need spatial_shuffle on both inputs? Yes, when both are large. Shuffling only one side leaves the other in arbitrary row order, so partition bounding boxes never line up and the executor falls back to comparing many partition pairs — you pay the shuffle and keep the all-pairs cost. The exception is the broadcast case in step 6, where the small side is not partitioned at all.

How do I decide between shuffling and broadcasting? Compare the size of the smaller layer against one worker's memory_limit. Under roughly a tenth of it, broadcast: every worker holds a copy, and you skip a full data movement. Above it, the replication costs more than the shuffle it avoids, and the shuffle also leaves you with an aligned dataset that later jobs reuse.

Why within and not intersects? For a strict points-in-polygons assignment, within matches a point only when it lies inside a polygon, avoiding the boundary double-counting that intersects produces. Choose the predicate that matches your question, then handle boundary ties explicitly.

Should I use Dask-GeoPandas or PostGIS for this join? If the data is shared, updated, and queried by many clients, load it into PostGIS and let its GiST index serve the join. For a one-shot batch over files you already hold, Dask-GeoPandas avoids standing up a database. For in-process SQL analytics over local GeoParquet, DuckDB Spatial Analytics is the lighter tool.

Can I join a partitioned layer against one that is not partitioned at all? Yes — that is exactly the broadcast path. Read the small side with plain GeoPandas, scatter it, and join each partition against it inside map_partitions. Do not convert it to a Dask collection first just to make the types match: a one-partition collection reintroduces the alignment problem for nothing.

How do I get the number of points per polygon without materialising the join? Chain the aggregation onto the lazy join and only compute the reduction, as in step 4. groupby(...).size() produces one row per zone — small enough to pull into the driver — and the executor never has to hold the full matched-pair table anywhere. Materialise the pair table only when you actually need per-pair attributes.

Does the joined output keep the left side's partition layout? The row-to-partition mapping survives, but the sizes do not: each output partition holds however many matched pairs its input partition produced. When the fan-out is uneven — dense zones producing more matches — the output is skewed even though the input was balanced. Repartition before writing if downstream jobs depend on even files.

Related