Scaling with Dask-GeoPandas

When a GeoDataFrame no longer fits in memory, Dask-GeoPandas partitions it across cores — or across a networked cluster of machines — and runs familiar GeoPandas operations on each chunk in parallel. It keeps the API you already know from GeoPandas DataFrames Explained while adding a spatial partition model so joins and overlays stay local instead of degrading to all-pairs comparisons. This stage of Spatial Analysis & Advanced Query Techniques covers that partition model, the one operation (spatial_shuffle) that makes distributed spatial work tractable, and where the library fits beside the two other out-of-core engines in this section: PostGIS Integration with Python for indexed, transactional access, and DuckDB Spatial Analytics for embedded columnar queries. Reach for Dask-GeoPandas when the work is embarrassingly parallel per feature and the dataset is too large to hold in a single process.

Dask-GeoPandas spatial partitioning A large GeoDataFrame is split into spatially partitioned chunks, each processed in parallel by a worker, and the partial results are concatenated back into one dataset. Large dataset > RAM partition 1 → worker partition 2 → worker partition n → worker Combined result concatenated spatial_shuffle aligns partitions by location so joins stay local
Dask-GeoPandas splits a dataset into spatially coherent partitions, processes them in parallel, and reassembles the output.

Architecture & Data Structures

A dask_geopandas.GeoDataFrame is not a table of geometries — it is a lazy collection of GeoPandas partitions plus a task graph that describes the computation to run over them. Nothing executes when you build it. The graph only runs when you force materialisation with .compute() (collect the whole thing back into a single in-memory GeoDataFrame) or stream it to disk with .to_parquet(). This laziness is the whole point: it lets you describe a pipeline over a dataset far larger than RAM, then have the scheduler walk it partition-by-partition so only a few chunks are resident at once.

Each partition is an ordinary GeoPandas GeoDataFrame backed by a Shapely 2.0 geometry array. Because the type is identical, any operation that is element-wise — buffer, area, centroid, to_crs, attribute filters, column arithmetic — maps across partitions with no coordination and parallelises for free. The library also tracks partition-level bounding boxes (spatial_partitions, a GeoSeries of one box per partition), which is the metadata that lets a spatial join prune whole partitions before touching a single geometry.

import dask_geopandas as dgpd

# Lazily open a partitioned GeoParquet dataset — no geometry is read yet
buildings = dgpd.read_parquet("buildings_partitioned/")

print(buildings.npartitions)                     # 64
print(type(buildings.partitions[0].compute()))   # <class 'geopandas.geodataframe.GeoDataFrame'>
print(buildings.spatial_partitions is not None)  # True once partitions carry bbox metadata

Three pieces of metadata ride along with that plan, and knowing which one is missing explains most confusing behaviour. _meta is an empty GeoDataFrame carrying the column names, the dtypes and the CRS; Dask consults it to answer schema questions without running anything, which is why buildings.dtypes returns instantly on a dataset it has never read. divisions describes how the index is split and is usually a tuple of None values unless something sorted the collection. spatial_partitions is the GeoSeries of one bounding box per partition — the only one of the three that carries geography, None until something computes it, and the metadata every partition-pruning optimisation depends on. A collection with spatial_partitions = None is perfectly usable; it simply never skips a partition, so the work it avoids is zero.

The mental model to hold onto: a Dask-GeoPandas object is a plan, and every partition is a real GeoDataFrame the plan will eventually run. Everything that follows is about keeping that plan cheap — chiefly by making sure spatially related features live in the same partition.

Lifecycle of a lazy Dask-GeoPandas collection Three phases run left to right. In the build phase, read_parquet, column assignment and a filter only append nodes to a task graph and no geometry is touched. The trigger phase is the call to to_parquet or compute that submits the graph. In the execute phase the scheduler walks the graph partition by partition: two partitions have already been flushed to Parquet, two are resident in the worker, and four are still unread on disk. A memory trace underneath stays flat during the build phase and then becomes a sawtooth, rising and falling once per partition rather than climbing to the size of the whole dataset. Nothing runs until the graph is triggered — then it runs one partition at a time 1 · Build the plan read_parquet → assign → filter graph nodes only, no geometry read 2 · Trigger .to_parquet() / .compute() graph handed to the scheduler 3 · Execute partition by partition the scheduler walks the graph in order only a few chunks are ever resident resident window p1 p2 p3 p4 p5 p6 p7 p8 flushed to Parquet, freed in worker memory now still unread on disk memory held in the worker flat — the plan holds no geometry sawtooth — one partition in, one partition out, never the whole dataset
Building the collection only grows a task graph; the memory cost appears at trigger time and stays bounded because the scheduler streams partitions rather than materialising the dataset.

Environment Configuration & Dependency Resolution

conda install -c conda-forge "dask-geopandas=0.4.*" "geopandas>=1.0" "dask>=2024.1" "pyarrow=15.*"
# Scaling beyond one machine adds the distributed scheduler and dashboard:
conda install -c conda-forge "distributed>=2024.1" "bokeh>=3.1"

pyarrow is not optional. Partitioned GeoParquet is the on-disk format that makes Dask-GeoPandas practical — it stores each partition as a separate row group or file, carries the CRS in metadata, and lets the reader skip partitions by their bounding box — and that I/O path runs entirely through Arrow. Pin pyarrow to the same major version across every process that reads or writes the dataset; a writer/reader skew is the most common source of opaque Parquet errors in a distributed cluster.

For a single multi-core workstation the default threaded scheduler is usually the right choice, because the heavy lifting happens inside GEOS, which releases the GIL. You only need dask.distributed when one box runs out of cores or RAM, or when you want the live task dashboard. Every worker imports the full GEOS/PROJ stack, so the same binary-dependency caveats from installing and configuring GeoPandas apply to each one — keep the environment identical across workers or a mismatched PROJ build will surface as silent reprojection differences between partitions.

# Single machine: threads are the default and usually best for GEOS-bound work.
# Multi-machine or when you want the dashboard:
from dask.distributed import Client

client = Client(n_workers=8, threads_per_worker=1, memory_limit="4GB")
print(client.dashboard_link)   # http://127.0.0.1:8787/status

Three version boundaries change the answers on this page, and all three fail in ways that do not name the real culprit.

The dask query planner. From dask 2024.3 the expression-based query planner (dask-expr) became the default backend for dask.dataframe. Dask-GeoPandas builds released before that change cannot see through the new expression layer: you get an import-time failure, or — worse — a collection that silently loses its geometry dtype at the first map_partitions. The switch is read exactly once, when dask.dataframe is first imported, so setting it after any Dask import does nothing at all and looks like the flag is broken. Move dask and dask-geopandas as a pair; the config escape hatch is for the window when you cannot.

GeoPandas 1.0. The sjoin argument op= was removed in favour of predicate=, the bundled example datasets were dropped, and pyogrio became the default I/O engine. Dask-GeoPandas forwards these calls straight to the per-partition GeoPandas method, so a pipeline written against an older release raises inside a worker, and the traceback points at an internal lambda rather than at your line.

Cluster-wide uniformity. Every process — scheduler, each worker, and the client — must run the same versions of dask, distributed, geopandas, shapely and pyarrow. distributed warns about a mismatch at connect time and then behaves unpredictably, most often as pickling errors on geometry arrays. Build workers from one environment file or one container image, and assert it at startup rather than reading logs later.

# This must sit ABOVE the first `import dask.dataframe` / `import dask_geopandas`
import dask

dask.config.set({"dataframe.query-planning": False})  # only for older dask-geopandas

import dask_geopandas as dgpd
from dask.distributed import Client

if __name__ == "__main__":                 # required: Windows and macOS use spawn
    client = Client(n_workers=8, threads_per_worker=1, memory_limit="4GB")
    client.get_versions(check=True)        # raises if a worker's stack differs
    parcels = dgpd.read_parquet("parcels_partitioned/")
    print(parcels.npartitions)

Vectorized Operations & Core Workflow

The canonical workflow is: open a partitioned dataset lazily, chain per-partition transforms, then trigger execution by writing the result back out. Element-wise operations need no shuffle and stay fully parallel. The one recipe that does require a shuffle — the distributed spatial join — has its own deep dive in Parallel Spatial Joins with Dask-GeoPandas.

import dask_geopandas as dgpd

# Parcels stored in a metric CRS (ETRS89 / UTM 32N) so area is in real m²
parcels = dgpd.read_parquet("parcels_partitioned/")   # CRS: EPSG:25832

# All three lines are per-partition, fully parallel, and still lazy —
# no geometry has been touched yet.
parcels["area_ha"] = parcels.geometry.area / 1e4
parcels["compact"] = parcels.geometry.length ** 2 / parcels.geometry.area
large = parcels[parcels["area_ha"] > 1.0]

# .to_parquet() is what finally runs the graph and streams results to disk,
# one partition at a time, without ever holding the whole dataset in memory.
large.to_parquet("large_parcels/", write_index=False)

Two habits keep this fast. First, prefer .to_parquet() over .compute() whenever the result is itself large — .compute() pulls everything into the driver process and defeats the purpose. Reserve .compute() for reductions (a count, a total area, a single small answer). Second, do your filtering as early in the chain as possible: a predicate like area_ha > 1.0 shrinks every partition before the expensive operations downstream see it.

Not every GeoPandas method has a Dask-GeoPandas wrapper. overlay, make_valid, explode, sample_points, anything from a third-party library — for all of these the escape hatch is map_partitions, which applies an ordinary function to each partition and stitches the results back into a lazy collection. The one thing that reliably goes wrong is meta. Dask must know the output schema before it runs anything, so it either infers it by calling your function on an empty frame — which for geometry work either raises or lies — or it trusts the meta you hand it. Give it a plain pandas.DataFrame and you get back a Dask DataFrame whose geometry column is object dtype: .geometry raises AttributeError: No geometry data set yet, and .to_parquet writes a column of WKB blobs with no GeoParquet metadata and no CRS. Nothing warns you until a downstream reader gets a table it cannot interpret as spatial. Build meta by running the function on a real sample and truncating it, so the dtypes and the CRS come from the same code path that will run in production.

import geopandas as gpd
import dask_geopandas as dgpd

# A small reference layer held on the driver; parcels are far too large for RAM
floodplain_boundary = gpd.read_file("floodplain_boundary.gpkg").to_crs(epsg=25832)
parcels = dgpd.read_parquet("parcels_partitioned/")     # EPSG:25832


def clip_to_floodplain(part: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Repair, then intersect one partition against the floodplain."""
    part = part.copy()
    part["geometry"] = part.geometry.make_valid()       # Shapely 2.0 / GEOS 3.8+
    return part.overlay(floodplain_boundary, how="intersection", keep_geom_type=True)


# Run the real function on a real sample, then truncate: correct dtypes AND CRS
sample = parcels.partitions[0].head(50)
meta = clip_to_floodplain(sample).head(0)

flooded = parcels.map_partitions(clip_to_floodplain, meta=meta)
flooded.to_parquet("parcels_in_floodplain/", write_index=False)

Two constraints govern what may go inside that function. It is serialised and shipped to every worker, so anything it closes over travels with it — a few-megabyte floodplain_boundary is fine, a national road network is not, and the fix for the second case is client.scatter() so the object moves once rather than once per task. More importantly, the function sees only its own rows. Any operation whose answer depends on features in another partition — a dissolve that spans a boundary, a nearest-neighbour lookup, a running total, a rank — returns a per-partition answer that looks entirely plausible and is wrong. Those operations need the shuffle described next.

Geometry / Data Processing Details

The decisive operation in the whole library is spatial_shuffle. By default, partitions are whatever the reader produced — often row-order chunks with no spatial meaning, so partition 3 might hold parcels from every corner of the map. That is fatal for a spatial join or a cross-partition dissolve: with no spatial alignment, every partition of the left side must be compared against every partition of the right, which is O(n²) in partitions and throws away all the parallelism you came for.

spatial_shuffle fixes this by repartitioning on a space-filling curve (Hilbert by default, Morton optional). It sorts every feature by its position along the curve so that spatially close features land in the same — or an adjacent — partition, then rebuilds the partition bounding boxes. After the shuffle, a join only has to compare partitions whose boxes actually overlap, and the scheduler prunes the rest.

import dask_geopandas as dgpd

sensors = dgpd.read_parquet("sensors_partitioned/")   # EPSG:25832

# Repartition by spatial proximity (Hilbert curve) BEFORE any spatial join.
# Shuffle both sides of a join with a comparable npartitions so their
# bounding boxes line up.
sensors = sensors.spatial_shuffle(npartitions=64)

# The shuffle is itself lazy; inspecting the new partition boxes forces it.
print(sensors.spatial_partitions.head())
Partitioning before and after spatial_shuffle Two map grids of sixteen cells coloured by partition. On the left, row-order partitions scatter each partition's cells across the whole map, so every partition's bounding box spans the full extent and overlaps every other, forcing all-to-all join comparisons. On the right, spatial_shuffle sorts features along a Hilbert curve so each partition becomes a compact non-overlapping tile and a join only compares neighbouring partitions. Row-order partitions a chunk holds features from anywhere After spatial_shuffle features sorted along a Hilbert curve each partition's bbox ≈ the full extent four compact, non-overlapping tiles Every pair overlaps Only neighbours overlap 6 partition pairs to check 4 partition pairs to check
Row-order chunks scatter each partition across the map, so every bounding box overlaps and a join degrades to all-to-all. spatial_shuffle sorts features along a Hilbert curve into compact tiles, so a join only compares partitions whose boxes actually touch.

The join and dissolve semantics are exactly the in-memory ones covered in Spatial Joins & Merging; Dask-GeoPandas only changes how the work is distributed, not what the predicates mean. The cost you are trading is a full data movement — the shuffle physically reorders rows across partitions — so do it once, before a run of spatial operations, not repeatedly. Choosing how to lay those partitions out, and measuring whether the layout you were handed is any good, is its own decision, worked through in Partitioning Strategies for Dask-GeoPandas.

A dissolve is a shuffle in disguise. dissolve(by="district_id") cannot union anything until every row sharing a key sits on one worker, so Dask-GeoPandas implements it as a groupby: a partial aggregation inside each partition, a shuffle on the key, then the final union. The parameter that decides whether it survives is split_out, the number of partitions the result gets. The default of 1 funnels every group into a single output partition — fine for a few hundred districts, fatal for half a million postcodes, because that one partition has to hold every dissolved geometry at once and the worker holding it dies while the rest of the Dask cluster idles.

import dask_geopandas as dgpd

parcels = dgpd.read_parquet("parcels_hilbert/")    # already spatially shuffled

# split_out sizes the RESULT, not the input. Raise it when the key has
# many distinct values; leave it at 1 only for a genuinely small group count.
districts = parcels.dissolve(by="district_id", split_out=16)

# Area of the merged shape — computed after the union, deliberately.
districts["area_ha"] = districts.geometry.area / 1e4
districts.to_parquet("districts/", write_index=True)

The ordering in those last two lines is a real decision, not style. area after the dissolve measures the merged footprint; a groupby("district_id")["area_ha"].sum() before it measures the total of the parts. They agree only when the parts neither overlap nor share slivers, and the cheap version is the right one whenever you do not need the merged geometry at all.

There is no distributed sjoin_nearest. Nearest-neighbour matching is the operation the partition model genuinely cannot express: the closest feature to a point sitting near a partition edge may live in the next partition, and a task only ever sees its own rows. Run sjoin_nearest per partition through map_partitions and every row gets an answer — quietly wrong for the ones near a boundary, with an error that grows with the gap between the true match and the edge. The workable pattern is a halo. Bound the search distance, buffer each partition's own bounding box by that distance, pull the candidate features falling inside the halo, and join the partition against partition-plus-halo. The single-machine mechanics of the operation itself are in Nearest Neighbor & KD-Tree Search.

import geopandas as gpd
import dask_geopandas as dgpd
from shapely.geometry import box

MAX_SEARCH_M = 500          # never match a hydrant further away than this

sensors = dgpd.read_parquet("sensors_hilbert/")        # EPSG:25832, out of core
hydrants = gpd.read_parquet("hydrants.parquet")        # small enough to hold in RAM
hydrant_sindex = hydrants.sindex


def nearest_hydrant(part: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    if part.empty:
        return part.assign(index_right=-1, hydrant_id="")
    # Halo = this partition's extent grown by the maximum search distance
    halo = box(*part.total_bounds).buffer(MAX_SEARCH_M)
    candidates = hydrants.iloc[hydrant_sindex.query(halo)]
    return gpd.sjoin_nearest(
        part, candidates, how="left", max_distance=MAX_SEARCH_M
    )


sample = sensors.partitions[0].head(50)
matched = sensors.map_partitions(
    nearest_hydrant, meta=nearest_hydrant(sample).head(0)
)
matched.to_parquet("sensors_with_hydrant/", write_index=False)

The halo is only correct because the radius is bounded — every candidate that could win is inside it by construction. An unbounded nearest join ("the closest hydrant, however far away") has no partition-local answer at all: push that one into PostGIS where a GiST index serves it as an indexed <-> query, or reduce the candidate layer until it fits on one machine.

CRS Alignment & Projection Pipeline

Coordinate reference systems are the quietest failure mode at scale, because a mismatch produces wrong numbers rather than an exception, and at partition granularity you may never eyeball the bad rows. Two rules keep results correct.

Reproject before you partition or shuffle. A Hilbert ordering computed in degrees is not the same ordering you get in metres, and area or distance computed in a geographic CRS is meaningless — so establish the canonical projected CRS first. Use a proper metric CRS for the region (here ETRS89 / UTM 32N, EPSG:25832), never Web Mercator (EPSG:3857) for measurement. Set the CRS conventions with Coordinate Systems with PyProj.

Keep both sides of a join in the same CRS, and assert it before computing. to_crs runs per-partition and in parallel, but if two datasets carry different CRSs the join will still run and silently return garbage matches.

import dask_geopandas as dgpd

roads = dgpd.read_parquet("roads_partitioned/")   # EPSG:4326 (lon/lat)

# Reproject the whole lazy collection to a metric CRS in parallel.
roads_utm = roads.to_crs(epsg=25832)

# Verify on a single partition without materialising the entire dataset.
sample = roads_utm.partitions[0].compute()
assert sample.crs.to_epsg() == 25832, "reprojection did not take"
assert sample.crs.is_projected, "metric CRS required before area/length"

Checking one partition with .partitions[0].compute() is the cheap way to validate a CRS invariant on a huge dataset — you confirm the projection took without paying to materialise all 64 partitions.

One trap belongs to the lazy API specifically. estimate_utm_crs() picks a zone from whatever extent it is shown, so calling it inside map_partitions gives each partition its own zone. After a spatial shuffle that is the worst possible arrangement: partitions are compact and geographically separate, which is exactly the condition under which every one of them picks a different answer. The collection then holds partitions in mutually incompatible grids, to_parquet records the CRS it saw first, and coordinates disagree by hundreds of kilometres with no exception anywhere. Resolve the zone once, on the driver, from a sample — then hand that fixed CRS to the parallel call.

import dask_geopandas as dgpd

sensors = dgpd.read_parquet("sensors_wgs84/")     # EPSG:4326, lon/lat

# Resolve ONE zone on the driver from a cheap sample...
zone = sensors.partitions[0].head(1000).estimate_utm_crs()
print(zone.to_epsg())                             # 32632

# ...then apply that single fixed CRS across every partition.
sensors_m = sensors.to_crs(zone)
assert sensors_m.crs.to_epsg() == zone.to_epsg()

The sampling shortcut is only safe when the first partition is representative of the extent, which after a Hilbert shuffle it is emphatically not — a shuffled partition covers one corner of the map. Take the estimate from the unshuffled collection, or compute total_bounds across the whole collection first. And if the dataset genuinely straddles several zones, no UTM zone is right at any granularity: use one equal-area CRS over the whole extent and accept its distortion uniformly. The zone-selection logic itself is covered in choosing a UTM zone automatically in Python.

Production Export & Integration

Landing the result in a database. There is no distributed to_postgis, and you would not want one: a hundred workers opening a hundred connections is how a Dask cluster exhausts max_connections and takes the database down with it. Write partition by partition instead, with a deliberately bounded worker count, and build the SQLAlchemy engine inside the task — engines hold sockets and do not survive being pickled and shipped to a worker.

import geopandas as gpd
import pandas as pd
import dask_geopandas as dgpd
from sqlalchemy import create_engine

DSN = "postgresql+psycopg://gis@db.internal:5432/parcels"


def write_partition(part: gpd.GeoDataFrame) -> pd.DataFrame:
    engine = create_engine(DSN)          # built in the task, never serialised
    part.to_postgis("flood_parcels", engine, if_exists="append",
                    index=False, chunksize=10_000)
    engine.dispose()
    return pd.DataFrame({"rows": [len(part)]})


flooded = dgpd.read_parquet("parcels_in_floodplain/")
written = flooded.map_partitions(
    write_partition, meta=pd.DataFrame({"rows": pd.Series(dtype="int64")})
).compute()
print("rows written:", int(written["rows"].sum()))   # rows written: 4180233

Create the table once up front and add its GiST index after the load: building the index in one pass at the end is far cheaper than maintaining it through several million inserts, and it removes the write contention that otherwise turns parallel workers into a queue. The connection and dtype mechanics are in Connecting GeoPandas to PostGIS with SQLAlchemy, and the index itself in Spatial Indexing in PostGIS with GiST.

Where the ceiling is. Dask-GeoPandas parallelises the per-feature part of a workload and nothing else, and three separate costs put a ceiling on it. Task overhead comes first: the scheduler spends on the order of a millisecond per task, so a stage whose partitions each finish in 50 ms burns more time being coordinated than computed — the remedy is fewer, larger partitions, never more workers. The shuffle comes second: it is the one stage that moves every row, it scales with network and disk rather than with cores, and past the point where the interconnect saturates, adding workers makes it slower. The driver comes third: every .compute() funnels its result through a single process, so a pipeline that ends in a large .compute() keeps a single-machine memory ceiling however big the Dask cluster is.

The band where the library earns its overhead starts at roughly ten million features, or a working set several times one machine's RAM. Below that, a single GeoPandas process running Shapely 2.0's vectorized operations usually wins outright — the array-level speedups described in Shapely 1.x vs Shapely 2 vectorization removed most of the per-geometry Python overhead that made distribution attractive in the first place. Measure the single-machine version before you distribute anything; it is a common outcome that the honest benchmark ends the project.

Windows / Platform Edge Cases & Debugging

Most Dask-GeoPandas problems are performance cliffs or silent wrong answers rather than crashes, so the debugging discipline is to assert invariants and read the task dashboard rather than trust the output.

Frequently Asked Questions

When should I use Dask-GeoPandas instead of plain GeoPandas? Only when the dataset genuinely exceeds comfortable memory, or a per-feature operation is CPU-bound enough to want every core. Below a few million features, in-memory GeoPandas is simpler and often faster once you account for scheduler overhead — Dask earns its keep on tens of millions of features and up.

Do I always need spatial_shuffle? No. Element-wise work (buffer, area, to_crs, attribute filters) never needs it. You need a shuffle only before operations that compare features across partitions — spatial joins, cross-partition dissolves, nearest-neighbour joins. Shuffle once, then run all of them.

Dask-GeoPandas or PostGIS for large spatial joins? If the data is shared, updated, and queried by many clients, put it in PostGIS and let its GiST index serve indexed joins. If it is a one-shot batch computation over files you already have, Dask-GeoPandas avoids standing up a database. The PostGIS integration guide covers the server-side path.

How is this different from DuckDB Spatial? DuckDB is an in-process columnar engine that excels at SQL-style filtering, aggregation, and GeoParquet scans on a single machine. Dask-GeoPandas keeps the Python/GeoPandas object model and scales the same per-partition operations across many cores or machines. Use DuckDB for analytical queries, Dask for parallel geometric processing that stays in Python.

What partition size should I target? Roughly 100–300 MB per partition in memory. Too many small partitions swamp the scheduler with task overhead; too few large ones exhaust worker RAM. Set npartitions to land in that band and confirm on the dashboard that workers are not spilling.

How do I debug a failure that only happens on one partition? Bring it back to a single process. bad = collection.partitions[17].compute() gives you an ordinary GeoDataFrame you can pass straight to the failing function under a normal debugger, and dask.config.set(scheduler="synchronous") runs the whole graph in the calling thread so the traceback points at your line instead of arriving deserialised from a worker. Reproduce first, distribute second.

Does .persist() remove the need to write the shuffled dataset to disk? Only for the lifetime of that cluster. .persist() keeps the computed partitions in distributed memory, which is the right move when several operations in one script share an aligned dataset — but it evaporates when the Dask cluster shuts down, and it competes with the working set of the very jobs you are running. For anything reused across scripts or across days, write the aligned collection to partitioned GeoParquet and let the next job read it back.

Can I use Dask-GeoPandas for raster work? No — it partitions tabular vector data, and a raster has no rows to distribute. Chunked, larger-than-memory raster processing belongs to Dask's array side, which is what Xarray & rioxarray raster cubes wrap for geospatial use. The two compose in one script: cube in xarray, features in Dask-GeoPandas, joined at the zonal-statistics step.

Related