How to Scale a GeoPandas Job with dask-geopandas

Problem statement

dask-geopandas partitions a GeoDataFrame and runs operations across partitions, in parallel or out of core. It is the right tool for data that does not fit in memory, and it is slower than plain GeoPandas for data that does.

Measured on 954,489 Overture building footprints in Rotterdam and Den Haag (Example 1):

                      total area    buffer 50 m, then area
plain geopandas          0.093 s          13.781 s
dask,  1 partition       0.173 s          14.747 s
dask,  4 partitions      0.115 s           4.283 s
dask,  8 partitions      0.152 s           2.515 s
dask, 16 partitions      0.130 s           1.951 s
dask, 64 partitions      0.319 s           1.878 s

Summing areas, plain GeoPandas wins at every partition count, because the whole job is a tenth of a second and dask's coordination is a comparable cost. Buffering first, dask at 16 partitions is 7.1 times faster, because each polygon now carries real work for a core to do.

Knowing when that inverts is the whole skill.

Quick answer

Use dask-geopandas when at least one is true:

the data does not fit in memory
the per-row work is heavy (buffering, overlay, complex predicates)
the data is already partitioned on disk
the operation is embarrassingly parallel and takes minutes

Use plain GeoPandas otherwise:

import dask_geopandas as dgpd

ddf = dgpd.read_parquet("buildings/")      # one partition per file, read in parallel
result = ddf.geometry.buffer(50).area.sum().compute()
print(f"{ddf.npartitions} partitions; buffered area {result / 1e6:,.1f} km²")

8 partitions; buffered area 9,166.4 km²

The files in buildings/ were written in space-filling-curve order, so each partition already covers a compact area and no spatial_shuffle() is needed; on arbitrary row ranges, consider it before an expensive spatial join (Example 2 measures when it pays).

Prove the scale problem, size partitions, shuffle spatially

1. Establish whether you have a scale problem

import shapely

columns_mb = gdf.memory_usage(deep=True).sum() / 1e6
coords_mb = shapely.get_num_coordinates(gdf.geometry.values).sum() * 16 / 1e6
print(f"{len(gdf):,} rows; {columns_mb:,.0f} MB reported by memory_usage, "
      f"plus about {coords_mb:,.0f} MB of coordinates it cannot see")

954,489 rows; 83 MB reported by memory_usage, plus about 114 MB of coordinates it cannot see

memory_usage(deep=True) counts a geometry column as an array of pointers; the coordinates live inside GEOS objects it does not look into, so count them separately. At this size, with plenty of memory, the question is not whether the data fits but whether the work is heavy enough to share out — Example 1 measures both cases.

2. Size the partitions by memory, not by core count

A partition should be big enough that the work dwarfs the task overhead, and small enough that several fit in memory at once.

Around 100 MB to 1 GB per partition is a reasonable range. npartitions = cores is a common default and usually wrong — it optimises for parallelism and ignores memory.

3. Spatially shuffle before any spatial operation

ddf = ddf.spatial_shuffle()
print(ddf.npartitions, "partitions, each covering a compact area")

8 partitions, each covering a compact area

By default, partitions are arbitrary row ranges, so a spatial join must compare every left partition with every right one. After a spatial shuffle, rows are grouped by location, and only overlapping partitions need comparing.

It reduces the number of partition pairs a join has to consider — measured in Example 2, from 256 to 86 for a self-join at 16 partitions. Whether that saves time depends on the size of the join; at a million buildings, it did not.

4. Keep the computation lazy until the end

result = (ddf[ddf["subtype"] == "residential"]
            .assign(area=lambda d: d.geometry.area)
            .groupby("class")["area"].sum()
            .compute())
print((result.sort_values(ascending=False).head(4) / 1e6).round(2))       # km²

class
house          27.51
apartments     15.06
residential     0.73
garage          0.33
Name: area, dtype: float64

Every .compute() executes the graph from the start. Build the whole expression, compute once.

5. Prefer partitioned Parquet as the input

dgpd.read_parquet on a partitioned directory reads each file as a partition, in parallel, with no shuffle needed. Reading a single GeoPackage and repartitioning does the expensive part serially before the parallelism starts.

A scaled spatial join on Overture buildings for Rotterdam to Den Haag

The data used on this page

The buildings are Overture Maps' building footprints for Rotterdam, Delft and Den Haag — every footprint whose bounding box touches 4.2–4.6° E, 51.85–52.15° N, about 950,000 of them — read straight from Overture's public S3 bucket with DuckDB, projected to the Dutch RD New grid (EPSG:28992) so areas are in square metres, and saved twice: as one buildings.parquet file and as a buildings/ directory of 8 files written in Hilbert-curve order. The 65 towns and cities (Overture locality division areas) in the same box are the right-hand side of the join. Overture keeps only recent releases, so the setup looks up the newest one and the counts can move a little between releases.

import os
import re

import dask_geopandas as dgpd
import duckdb
import geopandas as gpd
import requests

BBOX = (4.2, 51.85, 4.6, 52.15)          # Rotterdam, Delft and Den Haag


def latest_overture_release():
    listing = requests.get("https://overturemaps-us-west-2.s3.amazonaws.com/",
                           params={"list-type": "2", "prefix": "release/", "delimiter": "/"},
                           timeout=60).text
    return sorted(re.findall(r"<Prefix>release/([^<]+)/</Prefix>", listing))[-1]


def overture(theme, kind, where, path):
    """Copy the features inside BBOX from Overture's public bucket to a local GeoParquet file."""
    if os.path.exists(path):
        return
    con = duckdb.connect()
    con.execute("install httpfs; load httpfs; install spatial; load spatial;")
    con.execute("set s3_region = 'us-west-2'; set enable_progress_bar = false")
    x0, y0, x1, y1 = BBOX
    release = latest_overture_release()
    con.execute(f"""
        copy (select * exclude (bbox) from read_parquet(
                's3://overturemaps-us-west-2/release/{release}/theme={theme}/type={kind}/*.parquet',
                hive_partitioning = 0)
              where bbox.xmin < {x1} and bbox.xmax > {x0} and bbox.ymin < {y1} and bbox.ymax > {y0}
                and {where})
        to '{path}' (format parquet)""")
    with open(path + ".release", "w") as note:
        note.write(release)


overture("buildings", "building", "true",
         "overture_buildings.parquet")
overture("divisions", "division_area", "subtype = 'locality' and country = 'NL'", "overture_localities.parquet")

if not os.path.exists("buildings.parquet"):
    raw = gpd.read_parquet("overture_buildings.parquet", columns=["id", "subtype", "class", "height", "geometry"])
    raw.to_crs(28992).to_parquet("buildings.parquet")          # RD New, metres
if not os.path.exists("localities.parquet"):
    places = gpd.read_parquet("overture_localities.parquet")
    places["name"] = places["names"].map(lambda n: n["primary"])
    places[["id", "name", "geometry"]].to_crs(28992).to_parquet("localities.parquet")
if not os.path.exists("buildings"):
    ordered = gpd.read_parquet("buildings.parquet")
    ordered = ordered.iloc[ordered.geometry.hilbert_distance().argsort()]
    dgpd.from_geopandas(ordered.reset_index(drop=True), npartitions=8).to_parquet("buildings")

gdf = gpd.read_parquet("buildings.parquet")
print(f"Overture release {open('overture_buildings.parquet.release').read()}: "
      f"{len(gdf):,} buildings, {len(gpd.read_parquet('localities.parquet'))} localities; "
      f"buildings/ holds {len(os.listdir('buildings'))} files")

Overture release 2026-09-23.1: 954,489 buildings, 65 localities; buildings/ holds 8 files

Example 1 — deciding whether dask is worth it

import time
import geopandas as gpd
import dask_geopandas as dgpd


def should_i_use_dask(gdf, operation, npartitions=8, sample=50_000):
    """Time the operation both ways on a sample and extrapolate."""
    subset = gdf.sample(min(sample, len(gdf)), random_state=0)

    start = time.perf_counter()
    operation(subset)
    plain = time.perf_counter() - start

    ddf = dgpd.from_geopandas(subset, npartitions=npartitions)
    start = time.perf_counter()
    result = operation(ddf)
    if hasattr(result, "compute"):
        result.compute()
    parallel = time.perf_counter() - start

    print(f"  on {len(subset):,} rows: plain {plain:.3f}s, "
          f"dask {parallel:.3f}s ({parallel / plain:.1f}x)")
    if parallel < plain:
        print("  -> use dask: the work per row is heavy enough to parallelise")
    else:
        print("  -> use plain geopandas: dask overhead exceeds the work")
    return parallel < plain

Two operations on the same buildings: summing their areas, which is one cheap GEOS call per polygon, and buffering each by 50 m first, which is a real geometric construction per polygon.

area = lambda frame: frame.geometry.area.sum()
buffered = lambda frame: frame.geometry.buffer(50).area.sum()
for label, operation in (("area", area), ("buffer 50 m", buffered)):
    print(label)
    should_i_use_dask(gdf, operation)

area
  on 50,000 rows: plain 0.022s, dask 0.054s (2.5x)
  -> use plain geopandas: dask overhead exceeds the work
buffer 50 m
  on 50,000 rows: plain 0.768s, dask 0.158s (0.2x)
  -> use dask: the work per row is heavy enough to parallelise

The sample settles the question in seconds. Here is the same comparison on all the buildings, at several partition counts:

for label, operation in (("area", area), ("buffer 50 m", buffered)):
    start = time.perf_counter()
    operation(gdf)
    print(f"{label:12} plain geopandas   {time.perf_counter() - start:6.3f} s")
    for n in (1, 4, 8, 16, 64):
        ddf = dgpd.from_geopandas(gdf, npartitions=n)
        start = time.perf_counter()
        operation(ddf).compute()
        print(f"{label:12} dask {n:3d} partitions {time.perf_counter() - start:6.3f} s")

area         plain geopandas    0.093 s
area         dask   1 partitions  0.173 s
area         dask   4 partitions  0.115 s
area         dask   8 partitions  0.152 s
area         dask  16 partitions  0.130 s
area         dask  64 partitions  0.319 s
buffer 50 m  plain geopandas   13.781 s
buffer 50 m  dask   1 partitions 14.747 s
buffer 50 m  dask   4 partitions  4.283 s
buffer 50 m  dask   8 partitions  2.515 s
buffer 50 m  dask  16 partitions  1.951 s
buffer 50 m  dask  64 partitions  1.878 s

Example 2 — a spatial join that scales

import dask_geopandas as dgpd
import geopandas as gpd


def scaled_spatial_join(left_path, right_path, predicate="intersects"):
    """Partitioned join; the small right side is read by GeoPandas."""
    left = dgpd.read_parquet(left_path)
    # Reading both sides with dgpd.read_parquet fails under pandas 3: each frame stores its
    # spatial partitions in .attrs, and the join compares the two with ==.
    right = dgpd.from_geopandas(gpd.read_parquet(right_path), npartitions=1)

    left.calculate_spatial_partitions()
    right.calculate_spatial_partitions()
    joined = dgpd.sjoin(left, right, predicate=predicate)
    print(f"  left {left.npartitions} partitions, right {right.npartitions}; "
          f"{joined.npartitions} partition pairs to join")
    return joined

joined = scaled_spatial_join("buildings", "localities.parquet")
counts = joined.groupby("name").size().compute().sort_values(ascending=False)
print(counts.head(3).to_string())

  left 8 partitions, right 1; 8 partition pairs to join
name
Rotterdam     181982
Den Haag      168965
Zoetermeer     59702

The partition-pair count is where the shuffle is supposed to help: a join compares every left partition with every right one whose bounding box overlaps it. Measured on a self-join of the buildings at 16 partitions — rows in random order, then spatially shuffled:

shuffled_rows = gdf.sample(frac=1, random_state=0).reset_index(drop=True)
start = time.perf_counter()
spatial = dgpd.from_geopandas(shuffled_rows, npartitions=16).spatial_shuffle()
spatial.calculate_spatial_partitions()
print(f"spatial_shuffle itself       {time.perf_counter() - start:6.2f} s")
arbitrary = dgpd.from_geopandas(shuffled_rows, npartitions=16)
arbitrary.calculate_spatial_partitions()
for label, frame in (("arbitrary row ranges", arbitrary), ("spatially shuffled", spatial)):
    pairs = dgpd.sjoin(frame, frame, predicate="intersects")
    start = time.perf_counter()
    n = len(pairs)
    print(f"{label:22} {pairs.npartitions:4} pairs, {n:,} matches in {time.perf_counter() - start:5.2f} s")

spatial_shuffle itself        11.14 s
arbitrary row ranges    256 pairs, 2,172,929 matches in  3.04 s
spatially shuffled       86 pairs, 2,172,929 matches in  5.80 s

reset_index matters: from_geopandas sorts by the index by default, so shuffling the rows without resetting it hands dask the original, spatially ordered file back.

The shuffle cut the pairs to a third and made the join slower — 5.80 s against 3.04 s, after 11.14 s spent shuffling — because at a million polygons each pair is quick and the shuffle is a full data movement. It pays when the pairs are expensive — larger data, heavier predicates, partitions that do not fit in memory together — which is another instance of the same rule.

Example 3 — writing partitioned output

import dask_geopandas as dgpd


def process_and_write(input_path, output_path, buffer_m=50):
    """Read partitions, transform, write partitions — nothing in memory."""
    ddf = dgpd.read_parquet(input_path)
    ddf = ddf.set_geometry(ddf.geometry.buffer(buffer_m))   # not ddf["geometry"] = ...
    ddf["area_m2"] = ddf.geometry.area
    ddf.to_parquet(output_path, write_index=False)
    print(f"  wrote {ddf.npartitions} partitions to {output_path}")
    return output_path

import os

start = time.perf_counter()
process_and_write("buildings", "buffered")
print(f"  {time.perf_counter() - start:.1f} s; {len(os.listdir('buffered'))} files, "
      f"{len(dgpd.read_parquet('buffered')):,} rows read back")

  wrote 8 partitions to buffered
  8.6 s; 8 files, 954,489 rows read back

set_geometry is not a style choice. Assigning the buffered series with ddf["geometry"] = ... builds without complaint and then fails when the partitions are written, with AttributeError: 'Series' object has no attribute 'area': inside the partitions the assigned column has become a plain pandas Series.

Partitioned in, partitioned out, and nothing ever holds the whole dataset. That is the case dask-geopandas is genuinely for: a dataset larger than memory, transformed row-wise.

to_crs on a dask GeoDataFrame reprojects per partition, which is correct — reprojection is row-wise. Operations that are not row-wise, such as a dissolve across the whole dataset, need a shuffle and are much more expensive.

Task overhead, partition counts and the shuffle

Why dask is slower on small data

Every dask operation builds a task graph, serialises the tasks, schedules them, and collects the results. That costs milliseconds per task regardless of how much work the task does.

Computing the area of 954,489 polygons took 0.093 s of actual work in the measured run. Split into 64 partitions, the scheduling overhead alone exceeds it — hence 0.319 s. Buffering the same polygons took 13.8 s, and there the overhead disappears into the work: 1.95 s at 16 partitions.

The crossover is where per-partition work exceeds per-task overhead. For cheap row-wise operations that needs very large partitions; for expensive ones — an overlay, a complex predicate — it arrives much sooner.

Why more partitions is not better

The measurement shows 64 partitions no faster than 16 for the buffer, and three times slower than 4 for the area sum, which is the overhead curve turning up.

Each partition adds tasks, and each task adds scheduling. Beyond the point where every core is busy, more partitions only add coordination.

Size partitions by memory: large enough that the work is substantial, small enough that several fit at once. Then check that the count is at least the number of cores, and stop.

Why the spatial shuffle matters

Dask partitions a GeoDataFrame by row order. Rows adjacent in the file may be anywhere on the map.

A spatial join then cannot rule out any partition pair: any left partition might contain a row that intersects any right partition. So the join degenerates to comparing all pairs.

spatial_shuffle reorders rows by a space-filling curve and repartitions, so each partition covers a compact region with a small bounding box. The join can then skip most pairs on bounding boxes alone — 170 of 256 in the self-join measured in Example 2. Skipping them saves time only when each pair is expensive; at a million buildings it was not, and the shuffle cost more than it saved.

Why partitioned Parquet is the right input

Reading a single large file is serial: one process decodes it, then repartitions in memory. The parallelism starts after the slowest part.

A partitioned directory gives each worker a file. Reading is parallel from the start, and if the partitions were written spatially sorted, the shuffle is unnecessary too.

That makes the storage layout part of the compute design, which is the recurring theme of every cloud-native format.

Cheap operations, one compute and the dashboard

• Plain GeoPandas is faster for cheap operations — an area sum over 954,489 polygons — and dask is faster for heavy ones: 7.1× on a 50 m buffer of the same polygons.
• Size partitions by memory, roughly 100 MB to 1 GB.
• More partitions is not better. 64 was no faster than 16 for the buffer, and slower for the area sum.
• spatial_shuffle only when the pairs are expensive. On these polygons it cut the partition pairs to a third and still made the join slower (Example 2).
• One .compute() at the end, not per step.
• Partitioned Parquet input makes reading parallel from the start.
• Row-wise operations parallelise well. Dissolve, sort and global joins need shuffles.
• The dask scheduler dashboard shows whether time goes to tasks or to coordination.
• High memory after a threaded run is not always a leak. After Example 1 this process held 9.5 GB with nothing left to hold; ctypes.CDLL("libc.so.6").malloc_trim(0) brought it back to 1.4 GB. glibc keeps memory freed by dask's worker threads in per-thread arenas, so on Linux a notebook that runs several of these examples can reach 14 GB before the last one.

Internal links

• How to write partitioned GeoParquet and query it with filters — the input layout
• dask-geopandas uses more memory or time than plain GeoPandas — when it goes wrong
• Threads, processes and the GIL: parallelism in Python GIS explained — the alternatives
• Why GeoPandas is slow: the four real bottlenecks — check these first
• How to speed up batch GIS jobs with parallel processing in Python — file-level parallelism
• Lazy loading explained: why xarray reads nothing until you ask — the same laziness for arrays
• How to process a very large GeoPackage in chunks with Python — manual chunking
• Spatial indexes explained: R-trees and why spatial joins are fast — what the shuffle enables

FAQ

Is dask-geopandas faster than GeoPandas?

It depends on the work, not only the size. On 954,489 polygons, plain GeoPandas summed the areas faster at every partition count; buffering them first, dask at 16 partitions was 7.1× faster. Dask pays off when the data exceeds memory or the per-row work is heavy.

How many partitions should I use?

Enough that each holds 100 MB to 1 GB. Sizing by core count is a common default and usually wrong.

Why is my spatial join so slow?

Count the partition pairs first (pairs.npartitions). On arbitrary row ranges every partition meets every other — 256 pairs at 16 partitions here. A spatial shuffle cuts that, but it is a full data movement: on these polygons the shuffled join was slower, and it pays only when each pair is expensive.

When should I use dask-geopandas?

When the data does not fit in memory, when the per-row work is expensive, or when the data is already partitioned on disk.

Does to_crs work on a dask GeoDataFrame?

Yes — reprojection is row-wise, so it parallelises cleanly per partition.

Why did adding partitions make it slower?

Beyond the point where every core is busy, extra partitions only add scheduling overhead. In the measurement here, 64 partitions summed areas in 0.319 s against 0.115 s for 4.

What input format should I use?

Partitioned GeoParquet. Reading a single large file is serial, so the parallelism starts only after the slowest step.