Partitioning Strategies for Dask-GeoPandas

The same join, the same machine, the same predicate: one run finishes in four minutes, the other is still going an hour later — and the only difference is which rows ended up sitting next to which. This guide is for anyone whose Dask-GeoPandas job turned out slower than the single-machine version it replaced, and it covers the layout decisions that fix it. It sits under Scaling with Dask-GeoPandas in Spatial Analysis & Advanced Query Techniques, and it builds the aligned dataset that Parallel Spatial Joins with Dask-GeoPandas then joins.

Why This Approach / What Goes Wrong

A partition layout answers two independent questions. How big is a partition sets the cost of one task — worker memory, spill pressure, scheduler overhead. What is inside a partition sets the cost of every operation that has to compare features across partitions: joins, overlays, cross-partition dissolves, nearest-neighbour searches. read_parquet answers neither question on your behalf. It inherits whatever split the writer produced, which is almost always ingest order — collection date, alphabetical county, the order an API paged its results. Row order is a perfectly good layout for a full scan and a disastrous one for anything spatial.

The failure is measurable before you run a single join, and the metric is worth internalising. Every Dask-GeoPandas collection can carry spatial_partitions: a GeoSeries holding one bounding box per partition, computed from that partition's total_bounds. Divide the sum of those box areas by the area of the dataset's overall extent and you get a coverage ratio. A layout where each partition covers a distinct tile scores close to 1.0. A row-order layout scores close to npartitions, because every partition's box is effectively the whole map. That ratio is roughly the multiplier on how many partition pairs the executor has to consider before it can start testing geometries — at 96 partitions and a coverage of 71, no pair can be pruned and the "distributed" join degenerates into work no single machine would ever have done.

Partition bounding boxes before and after a Hilbert shuffle Two map panels show the same twenty-four features. On the left, row-order partitions assign colours at random across the map, so all four partition bounding boxes are drawn as near-full-extent dashed rectangles stacked on top of each other; the summed box area is 3.9 times the extent. On the right, a Hilbert shuffle groups neighbouring features, so each partition's dashed box is a compact tile and the summed box area is only 1.2 times the extent. The ratio approximates how many partition pairs a join must consider. Score a layout before you run anything: box area ÷ extent area Row-order partitions After spatial_shuffle (Hilbert) Σ box area ÷ extent = 3.9× Σ box area ÷ extent = 1.2× Nothing can be pruned when every box spans the extent — a real 96-partition dataset scores near 71×
The coverage ratio is a one-line health check on a partition layout: near 1 means compact tiles, near npartitions means every box covers the map.

The fix — reordering features along a Hilbert space-filling curve so neighbours share a partition — is a full data movement: every row potentially crosses the network to a different worker. That cost is real, and it is why the shuffle is worth treating as a build step rather than something you do at the top of every script. Compute it once, write the reordered collection back out as partitioned GeoParquet, and every downstream job starts from an already-aligned dataset.

Amortising the shuffle across three jobs Two timelines. The upper timeline shows three jobs run against an unsorted dataset: each repeats read, a long shuffle, and a join, so the shuffle is paid three times and the run ends late. The lower timeline pays one layout pass — read, shuffle, write — and then runs the same three jobs with only a short read and a join each, finishing considerably earlier. The gap between the two finish lines is marked as time saved on every run after the first. Where the shuffle cost goes when three jobs share one dataset No persisted layout — every job re-shuffles read shuffle join read shuffle join read shuffle join job 1job 2job 3 One layout pass, then three cheap jobs read shuffle write read join read join read join built once, reusedjob 1job 2job 3 saved per run The shuffle moves every row across the network — pay it once, then let later jobs start from the aligned copy.
Treat the Hilbert ordering as a build artefact, not a step in every script: one shuffle, written to disk, serves every later job.

Prerequisites

conda install -c conda-forge "dask-geopandas=0.4.*" "geopandas>=1.0" \
  "dask>=2024.1" "pyarrow=15.*" "shapely=2.0.*"

Keep pyarrow pinned to one major version across every process that touches the dataset. A writer/reader skew is the most common cause of a layout that reads back without its partition boxes.

Step-by-Step Implementation

The worked example is a national address-point layer — roughly 42 million points delivered as 96 GeoParquet files in county-name order.

1. Measure the layout you were given before changing it.

import dask_geopandas as dgpd

address_points = dgpd.read_parquet("address_points_raw/")

rows = address_points.map_partitions(len).compute()
print("partitions:", address_points.npartitions)   # partitions: 96
print("rows total:", int(rows.sum()))              # rows total: 42180000
print("min / median / max rows:",
      int(rows.min()), int(rows.median()), int(rows.max()))
# min / median / max rows: 11204 214905 9880311

A 900× spread between the smallest and largest partition is the first problem, and it is entirely independent of the spatial one: whichever worker draws the 9.9-million-row partition becomes the straggler every stage waits on. Both problems have the same remedy, which is convenient — a spatial shuffle splits the sorted feature sequence into equal-count chunks, so it fixes the size skew and the locality in one pass. That is also why it is worth measuring first: if the sizes are already even and the coverage ratio is already near 1, the dataset was written by someone who did this, and re-shuffling buys nothing.

2. Size partitions from measured bytes, not from a rule of thumb.

Bytes per row is not a property of the library, it is a property of your geometry. A point layer costs tens of bytes per feature; a coastline layer with 40,000-vertex polygons costs kilobytes. Measure one partition and extrapolate. WKB length is the closest cheap proxy for coordinate memory — live Shapely 2.0 geometries cost roughly twice their WKB size once the GEOS objects and the object array are counted.

import math

sample = address_points.partitions[0].compute()

wkb_bytes = int(sample.geometry.to_wkb().apply(len).sum())
attr_bytes = int(sample.drop(columns="geometry").memory_usage(deep=True).sum())
bytes_per_row = 2 * (wkb_bytes + attr_bytes) / len(sample)

TARGET_BYTES = 192 * 1024**2          # ~192 MiB of live GeoDataFrame per partition
total_rows = int(address_points.shape[0].compute())
target_partitions = max(1, math.ceil(total_rows * bytes_per_row / TARGET_BYTES))

print(f"{bytes_per_row:.0f} bytes/row -> {target_partitions} partitions")
# 419 bytes/row -> 88 partitions

The target band is set by two costs pulling against each other. Each task carries scheduler overhead on the order of a millisecond, so partitions small enough to finish in under a tenth of a second waste more time being scheduled than being computed. Meanwhile a worker holds several partitions at once — input, intermediate, output — so the byte budget should be a fraction of memory_limit, not all of it.

Partition sizing matrix for 42 million address points A four-row table for a 42.18 million row dataset at 419 bytes per row. Sixty thousand rows per partition is 25 megabytes across 703 partitions and is scheduler-bound. Two hundred and forty thousand rows is 96 megabytes across 176 partitions and is comfortable. Four hundred and eighty thousand rows is 192 megabytes across 88 partitions and is the target band. One and a half million rows is 600 megabytes across 28 partitions, causing spilling and stragglers. A footer gives the sizing formula. Two costs pull in opposite directions 42.18 M rows measured at 419 bytes per row live rows per partition live size partitions what actually happens 60 k25 MB703 240 k96 MB176 480 k192 MB88 1.5 M600 MB28 scheduling costs more than the work comfortable · plenty of parallelism target band workers spill · one straggler stalls all npartitions = ceil(total_rows × bytes_per_row ÷ target_bytes) Measure bytes_per_row on one partition — it varies by orders of magnitude between points and polygons.
Row count is the knob you turn, but bytes and task count are the costs you are actually balancing.

3. Score the layout you have.

calculate_spatial_partitions() runs one pass over the data and attaches a GeoSeries of per-partition boxes in place. With those boxes you can compute the coverage ratio and count how many partition pairs actually overlap — the number a join will have to work through.

from shapely.geometry import box

address_points.calculate_spatial_partitions()      # in place: sets .spatial_partitions
parts = address_points.spatial_partitions          # one bounding box per partition

extent = box(*parts.total_bounds)
coverage = parts.area.sum() / extent.area

hits = parts.sindex.query(parts, predicate="intersects")
overlapping_pairs = int((hits[0] != hits[1]).sum() // 2)

print(f"coverage {coverage:.1f}x, overlapping pairs {overlapping_pairs}")
# coverage 71.4x, overlapping pairs 4560

4,560 is exactly 96 × 95 / 2 — every partition overlaps every other, so nothing can be pruned. Keep this snippet around as a regression check: run it on any dataset before you plan a large join, and treat anything above roughly 2× coverage as a layout that has not been built yet. It is also the honest way to compare shuffle settings, because it reports the property a join actually depends on rather than a wall-clock number that varies with cache state and worker count.

4. Reproject to a metric CRS, then shuffle onto the Hilbert curve.

Order the reprojection first. hilbert_distance ranks each feature by where its centre falls on a space-filling curve laid over the dataset's total bounds, so the ordering depends on the coordinate system it is computed in. In EPSG:4326 a degree of longitude shrinks toward the poles, which stretches the curve's cells into anisotropic strips; worse, the authority definition of EPSG:4326 is latitude-first, so any hand-rolled pyproj.Transformer needs always_xy=True to keep the coordinate order GeoPandas expects. Use an equal-area metric CRS for a continental extent — here NAD83 / Conus Albers (EPSG:5070) — or a local UTM zone for a regional one, and never Web Mercator (EPSG:3857), whose area distortion would bias the curve toward high latitudes. The projection mechanics live in Coordinate Systems with PyProj.

address_points = address_points.to_crs(epsg=5070)   # NAD83 / Conus Albers, metres

# level=16 gives a 65536 x 65536 grid over the extent — fine enough that
# distinct addresses rarely share a curve position.
aligned = address_points.spatial_shuffle(
    by="hilbert",
    level=16,
    npartitions=88,        # from step 2
)

print(aligned.npartitions)                     # 88
print(aligned.spatial_partitions is not None)  # True — the shuffle recomputes them

spatial_shuffle sets the collection's index to the curve position and rebuilds the partition boxes on the way out (calculate_partitions=True is the default). Passing by="morton" swaps in a Z-order curve, which is cheaper to compute but leaves longer jumps between consecutive cells; by="geohash" expects geographic coordinates and is the wrong choice on a projected dataset.

5. Persist the layout so the next job inherits it.

Writing the shuffled collection back out is what turns a one-off shuffle into a reusable asset. Dask-GeoPandas writes one file per partition and records each file's bounding box in the GeoParquet metadata, so read_parquet can rebuild spatial_partitions without touching a single geometry.

aligned.to_parquet("address_points_hilbert/", write_index=True)

# A later, unrelated job — no shuffle, no scan:
reopened = dgpd.read_parquet(
    "address_points_hilbert/",
    gather_spatial_partitions=True,     # default; reads boxes from metadata
)
print(reopened.npartitions)                       # 88
print(reopened.spatial_partitions.iloc[0].bounds)
# (-2278400.0, 1148700.0, -2016100.0, 1401950.0)

Keeping write_index=True preserves the Hilbert distance as the index, so a reader can restore divisions as well as boxes. This is also the point at which the format choice pays off — a directory of GeoParquet files carries CRS and per-file extent in its metadata, which is one of the practical gaps covered in GeoParquet vs Shapefile for Storage.

How a partition layout survives to disk and back A left-to-right flow. A shuffled collection of eighty-eight Hilbert-ordered partitions is written by to_parquet as one file per partition into a directory; each file carries a bounding box in its GeoParquet metadata. Reading the directory back with gather_spatial_partitions set to true restores the partition boxes from that metadata without reading any geometry, so the next job starts already aligned. The layout travels in the GeoParquet metadata shuffled collection 88 partitions to_parquet() one file per partition address_points_hilbert/ part.0.parquet bbox part.1.parquet bbox part.87.parquet bbox read_parquet() gather_spatial_partitions = True Partition boxes are restored from metadata — the next job starts aligned without reading one coordinate Files written by a tool that omits per-file boxes force a full pass through calculate_spatial_partitions() instead.
Writing the shuffled collection is what makes the layout durable: the boxes ride along in the file metadata and come back for free.

Verification

Three invariants matter: no rows were lost, the partitions landed in the target size band, and the layout actually got tighter.

from shapely.geometry import box

reopened = dgpd.read_parquet("address_points_hilbert/")

# 1. The shuffle moves rows, it never drops them
assert int(reopened.shape[0].compute()) == 42180000

# 2. Partition sizes are even — no straggler
rows = reopened.map_partitions(len).compute()
print("min / max rows:", int(rows.min()), int(rows.max()))
# min / max rows: 479318 479319
assert rows.max() / rows.min() < 1.05

# 3. The layout is compact and the CRS is still metric
parts = reopened.spatial_partitions
coverage = parts.area.sum() / box(*parts.total_bounds).area
print(f"coverage {coverage:.2f}x  (was 71.4x)")   # coverage 1.18x  (was 71.4x)
assert coverage < 1.5, "shuffle did not produce compact partitions"
assert reopened.crs.is_projected

A coverage that stays high after a shuffle almost always means the shuffle ran on geographic coordinates, or that a handful of very large geometries are inflating their partitions' boxes.

Edge Cases & Debugging

Frequently Asked Questions

How do I pick npartitions when I have no idea of the data size? Work backwards from bytes, not from core count. Measure bytes_per_row on a single materialised partition as in step 2, choose a target of roughly 100–250 MiB of live data per partition, and divide. Then sanity-check that the result is at least two or three times your total thread count, so the scheduler has slack to balance stragglers.

Does spatial_shuffle sort the rows inside each partition too? It sorts globally by the space-filling curve position and then splits that sorted sequence into partitions, so rows within a partition are in curve order as well. That is a side benefit for compression — neighbouring rows have similar coordinates — but the ordering the executor cares about is the partition boundary, not the row order inside it.

Do I need to re-shuffle after adding a column or filtering rows? No. Element-wise work never moves a row between partitions, so the layout survives to_crs, attribute arithmetic, and boolean filters. Only a size problem calls for action after a filter, and repartition fixes that without a full data movement.

Is level=16 always the right resolution? It is the default and it is right for continental extents in metres, where it resolves the bounding box into cells of a few tens of metres. Lower it only when the extent is small enough that the curve resolution outruns the coordinate precision, and remember the level is applied to the dataset's total bounds — a single outlier point on the far side of the world coarsens every cell.

Related