--- title: "Parallel raster processing in Python with dask and rioxarray: a whole Sentinel-2 tile without running out of memory" description: "Process large satellite rasters in Python without running out of memory: dask chunks, task graphs, diagnostics and real timings on Sentinel-2." url: https://www.technicalgis.com/tutorials/parallel-raster-processing-dask-python/ author: TechnicalGIS Team publisher: TechnicalGIS (https://www.technicalgis.com/about/) published: 2026-11-03 updated: 2026-11-03 language: en-GB --- # Parallel raster processing in Python with dask and rioxarray: a whole Sentinel-2 tile without running out of memory > Process large satellite rasters in Python without running out of memory: dask chunks, task graphs, diagnostics and real timings on Sentinel-2. - Level: Intermediate - Tested with: Python 3.12 - Source code: [Download scripts](https://www.technicalgis.com/code/parallel-raster-processing-dask-python/) ## Key takeaways 1. Opening sixteen 10980 × 10980 Sentinel-2 bands (3.86 GB of uint16, 7.72 GB as float32) with rioxarray and dask chunks took 0.1 seconds and read no pixels: you build the whole calculation first, as a task graph, and nothing runs until you call compute() or write the result. 2. Match dask chunks to the file. The Sentinel-2 COGs store 512 px internal tiles, so 1024, 2048 or 4096 px chunks read whole tiles. At 2048 px the median NDVI of 8 scenes was a graph of 6,136 tasks; at 256 px it was 314,346 tasks, and the scheduler overhead made it several times slower. 3. Eager NumPy needed 12.4 GB to take the median NDVI of only 4 scenes, and grew by about 1.15 GB per scene. Dask did all 8 scenes in 19 seconds with a peak of 6.4 GB, and gave the same mean NDVI to four decimal places. 4. Chunk size and threads both matter, and so does rioxarray's read lock. The fastest run was 1024 px chunks on 16 threads (12.8 s); 256 px chunks took 55 s on 16 threads, 4096 px chunks needed 21 GB, and rioxarray's default lock, which lets only one thread read at a time, made the 2048 px run five times slower (108 s against 21 s). 5. Writing the 482 MB median straight to a tiled, compressed GeoTIFF with rio.to\_raster(windowed=True) took 21 seconds and raised memory by at most 2.5 GB, against 6.6 GB for computing the same result into memory. ## Before you start - Comfortable with NumPy and xarray; the Sentinel-2 Lake Mead tutorial introduces STAC searches and Cloud Optimised GeoTIFFs - Python 3.10+ with dask, distributed, xarray, rioxarray, rasterio, matplotlib, networkx, psutil, pystac-client and planetary-computer - About 5 GB of free disk space and 8 GB of RAM or more; the script downloads 4.6 GB of Sentinel-2 bands once, with no account or key **Short answer:** to process a satellite raster that is too big for memory in Python, open it with `rioxarray.open_rasterio(path, chunks={"band": 1, "x": 2048, "y": 2048})`. You get a lazy dask array instead of a NumPy array. Write the calculation with ordinary xarray operations, which builds a task graph, then either call `.compute()` on a result that fits in memory, or write it straight to disk with `.rio.to_raster(path, tiled=True, windowed=True, lock=threading.Lock())`. Choose chunks that are whole multiples of the file’s internal tiles and roughly 10–100 MB each, and keep the axis you reduce over (here, time) in one chunk. A single Sentinel-2 band is 10980 × 10980 pixels: 120 million values, 241 MB as uint16 and 482 MB as float32. Stack a few dates and bands, convert to float for an index, and a plain NumPy script runs out of memory long before it runs out of data. This tutorial takes **eight full Sentinel-2 scenes of tile 30UXC** (from Oxford and Aylesbury to central London, March to November 2025), computes NDVI for each date and reduces them to a **median NDVI map**: 1.93 billion input pixels, processed on a laptop-sized memory budget. You will learn how to: 1. find and download whole Sentinel-2 tiles with **STAC**, and check how they are stored; 2. open them **lazily** with rioxarray and choose a **chunk size**; 3. build the median NDVI as a **task graph** and see what dask will run; 4. watch a run with **dask diagnostics**: the task stream, CPU and memory; 5. compare **peak memory** with the same calculation in eager NumPy; 6. **time** five chunk sizes and five thread counts, the read lock and a `LocalCluster`; 7. **write the result to disk** chunk by chunk as a tiled, compressed GeoTIFF. All outputs below are real results from [`dask_raster.py`](https://www.technicalgis.com/code/parallel-raster-processing-dask-python/) and its helper `raster_pipeline.py`, run on 3 November 2026 with Python 3.12, dask 2026.8.0, xarray 2026.9.0 and rioxarray 0.23.0 on a Windows 11 workstation with 16 cores (32 threads) and 68 GB of RAM. The whole script took 25 minutes, most of it in the timing grid of Step 6. ## Step 1: download eight whole Sentinel-2 scenes Dask works on remote files too, but benchmarks over the internet mostly measure the network. So the script downloads bands B04 (red) and B08 (near infrared) for eight low-cloud dates once, with a STAC search on the Microsoft Planetary Computer. Where two processings of the same date exist, it keeps the latest: `dask_raster.py` ```python catalog = pystac_client.Client.open(STAC_URL, modifier=pc.sign_inplace) def find_item(date): items = catalog.search(collections=["sentinel-2-l2a"], datetime=date, query={"s2:mgrs_tile": {"eq": TILE}, "s2:nodata_pixel_percentage": {"lt": 1}}).item_collection() return sorted(items, key=lambda i: i.properties["s2:generation_time"])[-1] jobs = [] for d in DATES: it = find_item(d) jobs += [(it.assets[b].href, band_path(d, b)) for b in ["B04", "B08"]] with ThreadPoolExecutor(4) as ex: fetched = sum(ex.map(lambda j: download(*j), jobs)) ``` `Output — Step 1` ```text 2025-03-06 S2B_MSIL2A_20250306T110829_R137_T30UXC_20250306T133711 tile cloud 3.7% no data 0.0002% 2025-04-02 S2A_MSIL2A_20250402T110651_R137_T30UXC_20250402T201716 tile cloud 1.0% no data 0.0001% 2025-04-05 S2B_MSIL2A_20250405T110619_R137_T30UXC_20250405T132414 tile cloud 0.3% no data 0.0000% 2025-04-30 S2C_MSIL2A_20250430T110641_R137_T30UXC_20250430T151423 tile cloud 6.5% no data 0.0000% 2025-05-02 S2A_MSIL2A_20250502T110701_R137_T30UXC_20250502T181413 tile cloud 2.1% no data 0.0000% 2025-06-19 S2C_MSIL2A_20250619T110641_R137_T30UXC_20250619T151116 tile cloud 0.2% no data 0.0000% 2025-07-11 S2A_MSIL2A_20250711T110651_R137_T30UXC_20250711T153113 tile cloud 5.2% no data 0.0000% 2025-11-21 S2B_MSIL2A_20251121T111259_R137_T30UXC_20251121T120801 tile cloud 4.9% no data 0.0005% 16 files on disk, 4.58 GB compressed each file: 10980 x 10980 px, uint16, deflate, internal tiles (512, 512), overviews [2, 4, 8, 16] uncompressed: 8 dates x 2 bands x 10980 x 10980 uint16 = 3.86 GB (7.72 GB as float32) ``` All eight dates come from relative orbit 137, which covers the whole tile, and all three Sentinel-2 satellites (2A, 2B and 2C) contribute. The line that matters most for dask is `internal tiles (512, 512)`: the file is stored in 512 × 512 blocks, and any read touches whole blocks. We chose clear dates on purpose, so that the tutorial is about dask rather than cloud. For a real composite you would mask cloud first; see [Sentinel-2 cloud and shadow masking in Python](https://www.technicalgis.com/tutorials/sentinel-2-cloud-masking-python/). ## Step 2: open lazily and choose a chunk size `rioxarray.open_rasterio` with a `chunks=` argument returns a dask-backed `DataArray`. It reads the file headers, not the pixels. We open each band of each date with the same chunks and concatenate them along a new time dimension: `raster_pipeline.py` ```python def open_stack(dates, chunk, lock=False): chunks = {"band": 1, "y": chunk, "x": chunk} red, nir = [], [] for d in dates: red.append(rioxarray.open_rasterio(band_path(d, "B04"), chunks=chunks, lock=lock).squeeze("band", drop=True)) nir.append(rioxarray.open_rasterio(band_path(d, "B08"), chunks=chunks, lock=lock).squeeze("band", drop=True)) time_index = xr.Variable("time", np.array(dates, dtype="datetime64[ns]")) return xr.concat(red, dim=time_index), xr.concat(nir, dim=time_index) ``` The `lock` argument matters more than it looks. rioxarray’s default puts one shared lock around all reads; we pass `lock=False` so that every thread can read and decompress at the same time. Step 6 measures the difference. `Output — Step 2` ```text opened 16 files in 0.10 s without reading any pixels Size: 2GB dask.array chunk size options for a 10980 x 10980 band: 256 px: 1849 chunks per band, 0.3 MB per float32 chunk, 2 MB per 8-date block for the median 512 px: 484 chunks per band, 1.0 MB per float32 chunk, 8 MB per 8-date block for the median 1024 px: 121 chunks per band, 4.2 MB per float32 chunk, 34 MB per 8-date block for the median 2048 px: 36 chunks per band, 16.8 MB per float32 chunk, 134 MB per 8-date block for the median 4096 px: 9 chunks per band, 67.1 MB per float32 chunk, 537 MB per 8-date block for the median ``` Three rules decide the chunk size: - **Align with the file’s blocks.** All five options are multiples of 512, so no chunk splits an internal tile. A 1000 px chunk would make neighbouring chunks decompress the same tiles twice. - **Aim for roughly 10–100 MB per chunk of the working data type.** Smaller chunks mean more tasks, and each task costs the scheduler a few hundred microseconds whatever its size. Larger chunks mean fewer opportunities to run in parallel, and more memory per thread. - **Think about the reduction.** A median needs every date of a pixel in one place, so the whole time axis goes into one chunk (`time: -1` in Step 3). With 2048 px chunks, one 8-date float32 block is 134 MB; with 4096 px it is 537 MB per thread, which adds up quickly with 16 threads. ![Left: an NDVI map of the 110 by 110 km Sentinel-2 tile 30UXC west of London, in red, yellow, green and blue, divided by black lines into a 6 by 6 grid of 2048 pixel chunks; the right-hand column and bottom row of chunks are narrow strips labelled 2048 by 740 and 740 by 2048. One central chunk is outlined in dark red. Right: that chunk at full 10 m resolution, showing individual fields, woods and villages, with dashed lines dividing it into 4 by 4 squares of 512 pixels.](https://www.technicalgis.com/_astro/01-chunk-grid.7zm15VUU_Ux05W.webp) _The 2048 px chunk grid over tile 30UXC (NDVI on 5 April 2025). Edge chunks are only 740 px wide, because 10980 is not a multiple of 2048. Each full chunk holds exactly 16 of the file's 512 px internal tiles, so dask never reads part of a tile._ ## Step 3: build the median NDVI as a task graph The calculation is ordinary xarray. Sentinel-2 L2A data processed with baseline 04.00 or later carry an offset of 1000 on every digital number, and 0 means no data: `raster_pipeline.py` ```python OFFSET = 1000 # processing baseline 04.00+ adds 1000 to every DN def median_ndvi(red, nir): red = red.where(red > 0).astype("float32") - OFFSET # DN 0 is no data nir = nir.where(nir > 0).astype("float32") - OFFSET ndvi = (nir - red) / (nir + red) ndvi = ndvi.where(np.isfinite(ndvi)) ndvi = ndvi.chunk({"time": -1}) # a median needs every date in one chunk out = ndvi.median("time", skipna=True).astype("float32") out.attrs = {} # drop the uint16 _FillValue/scale attributes return out ``` Calling `median_ndvi` returns instantly. The result is still a dask array, and its graph tells us how much work dask has planned: `Output — Step 3` ```text Size: 482MB dask.array 2048 px chunks: 6136 tasks in 81 layers; output 482 MB in 36 chunks 256 px chunks: 314346 tasks, 1849 output chunks 512 px chunks: 82296 tasks, 484 output chunks 1024 px chunks: 20586 tasks, 121 output chunks 2048 px chunks: 6136 tasks, 36 output chunks 4096 px chunks: 1546 tasks, 9 output chunks two dates in 2 x 2 chunks: 180 tasks as written, 52 after dask fuses chains of element-wise steps graph as written: 180 tasks, 228 dependencies, 14 tasks on the longest path task types: concatenate 16, elementwise 104, median 4, open file 4, read chunk 48, rechunk 4 ``` The number of tasks grows with the square of the chunk count: halving the chunk size from 512 to 256 px quadruples the graph to over 300,000 tasks. Before running anything, dask optimises the graph, fusing chains of element-wise steps (mask, subtract, divide) into single tasks; for the full 2048 px graph that took 6,136 tasks down to 1,513. The full graph is too big to draw, but the same pipeline on two dates and 2 × 2 chunks is small enough. You can draw your own with `med.data.visualize()` (it needs Graphviz), or with networkx as the script does: `dask_raster.py` ```python small = median_ndvi(*open_stack(DATES[:2], 5490)) dsk = dict(small.data.__dask_graph__()) G = nx.DiGraph() for k in dsk: G.add_node(k) G.add_edges_from((dep, k) for dep in get_dependencies(dsk, k) if dep in dsk) ``` ![A network diagram read from top to bottom. Four purple open-file nodes at the top lead to three rows of blue read-chunk nodes, a row of orange concatenate nodes, several rows of green element-wise nodes, four red rechunk nodes and finally four dark red median nodes at the bottom, joined by grey lines.](https://www.technicalgis.com/_astro/02-task-graph.DU9D_6S6_218TkP.webp) _The task graph for two dates in 2 × 2 chunks: 180 tasks before optimisation. Each output chunk of the median (bottom) depends only on the same chunk of each input, so the four columns of work can run in parallel and finish independently._ The shape of the graph is the key to why dask saves memory: each output chunk depends only on its own column of inputs. Dask can read one column, reduce it and release the inputs before it reads the next, so only a few columns need to be in memory at once. ## Step 4: watch the run with dask diagnostics Nothing has been read yet. `compute()` runs the graph on the threaded scheduler. Two context managers record what happens: dask’s own `Profiler` logs the start and end of every task on every thread, and a small psutil sampler from `raster_pipeline.py` records the process memory and CPU every 0.1 s: `dask_raster.py` ```python with dask.config.set(scheduler="threads", num_workers=8): with Profiler() as prof, RSSSampler(interval=0.1) as rprof, ProgressBar(dt=10, out=sys.stdout): result = med.compute() ``` `Output — Step 4` ```text [ ] | 0% Completed | 294.90 us [####################### ] | 59% Completed | 10.01 s [########################################] | 100% Completed | 20.01 s computed in 20.1 s: 1513 tasks; process memory 0.35 GB at start, 6.64 GB at peak; mean CPU 694% (800% = 8 busy threads) median NDVI: mean 0.584; 100.00% of pixels valid; 6,874 km² above 0.6 and 192 km² below 0.1 ``` ![Top: a task stream with eight horizontal rows, one per thread, filled for about 19 seconds with short blue bars for chunk reads, thin green bars for element-wise steps and longer dark red bars for the median, spread evenly across all eight threads. Bottom: process memory in blue rising from 0.35 GB to about 6 GB over 18 seconds and dropping to 4.4 GB at the end; CPU in orange jumps to between 650% and 900% within half a second and stays there until about 17.5 seconds, then falls to zero.](https://www.technicalgis.com/_astro/03-dask-diagnostics.Bg3saWtv_2epdl4.webp) _Dask diagnostics for the median NDVI of eight scenes on eight threads. The task stream (top) shows every thread busy reading chunks (blue) and taking medians (dark red) from start to finish. CPU (bottom, orange) stays near 800% until the last chunks finish; memory (blue) climbs as finished output chunks accumulate._ The diagnostics answer the two questions that matter for any dask job: - **Are all the threads busy?** Yes. The task stream has almost no gaps, and the mean CPU was 694% of the 800% available. Gaps in the stream would mean that tasks are waiting on each other, on a lock or on the disk. - **Where does the memory go?** The working set is a few chunks per thread. The memory line rises steadily because the result itself (482 MB of float32) is assembled in memory, and because the operating system does not hand memory back immediately; the peak of 6.6 GB is well under the 7.7 GB the float32 inputs alone would need. The computed median is a reasonable picture of southern England in 2025: mean NDVI 0.58, with 6,874 km² of the 12,056 km² tile above 0.6 (crops and grass in leaf) and 192 km² below 0.1 (water, roofs and roads). ## Step 5: peak memory, eager NumPy against dask To see what dask saves, `raster_pipeline.py` also does the calculation the eager way: read every band into a NumPy array, stack them, and call `np.nanmedian`. Each run starts in a fresh Python process and records the peak memory above that process’s starting point: `raster_pipeline.py` ```python def eager_median_ndvi(dates): def read(d, b): with rasterio.open(band_path(d, b)) as src: return src.read(1) red = np.stack([read(d, "B04") for d in dates]).astype("float32") nir = np.stack([read(d, "B08") for d in dates]).astype("float32") red[red == 0] = np.nan nir[nir == 0] = np.nan red -= OFFSET nir -= OFFSET ndvi = (nir - red) / (nir + red) del red, nir ndvi[~np.isfinite(ndvi)] = np.nan return np.nanmedian(ndvi, axis=0).astype("float32") ``` `Output — Step 5` ```text eager NumPy, 1 date(s): 18.9 s, peak 8.95 GB (mean NDVI 0.4903) eager NumPy, 2 date(s): 35.3 s, peak 9.92 GB (mean NDVI 0.5125) eager NumPy, 4 date(s): 62.9 s, peak 12.37 GB (mean NDVI 0.5475) dask, 1 date(s): 3.9 s, peak 2.12 GB (mean NDVI 0.4903) dask, 2 date(s): 6.1 s, peak 2.50 GB (mean NDVI 0.5125) dask, 4 date(s): 10.4 s, peak 4.12 GB (mean NDVI 0.5475) dask, 8 date(s): 19.0 s, peak 6.36 GB (mean NDVI 0.5836) eager peak grows by 1.15 GB per date: about 17 GB for 8 dates and 88 GB for the 70 or so scenes a year of this tile ``` ![Left: memory over time for three runs. The eager NumPy run with 4 dates, in red, climbs slowly in a saw-tooth to 3 GB over 43 seconds, then jumps in steps to 12.4 GB at 62 seconds. Two dask runs, with 4 and 8 dates, in green and blue, rise to about 4 and 6 GB and finish within 11 and 19 seconds. Right: peak memory against the number of dates. The eager NumPy line rises from 9 GB for 1 date to 12.4 GB for 4 dates, with a dotted trend reaching 17 GB at 8 dates; the dask line rises from 2.1 to 6.4 GB between 1 and 8 dates.](https://www.technicalgis.com/_astro/04-memory-eager-vs-dask.Zrd7oQcK_Z1vE006.webp) _Memory of eager NumPy and dask for the same median NDVI. Eager NumPy needed 12.4 GB for four scenes and was not run for eight; the dotted line extends its trend. Dask processed eight scenes in less memory than eager NumPy needed for one._ Three things stand out: - **The results agree.** The mean NDVI matches to four decimal places for every stack size, so dask changes how the work is scheduled, not what is computed. - **Eager NumPy is greedy even for one date.** Reading a single date (two bands) is about 0.5 GB of uint16, yet the peak was 9 GB. The float32 copies, the NDVI array and its temporaries account for a few gigabytes; the rest is GDAL’s block cache (see Step 6) and `np.nanmedian` itself, which copies the whole array before it sorts. Dask calls the same `nanmedian`, but on one 2048 px block at a time. - **The eager trend does not stop.** At 1.15 GB per date, a year of this tile (about 70 usable dates) would need about 88 GB. Dask’s peak grows too, because the result and a few blocks per thread are in memory at once, but the per-thread working set depends on the chunk size, not on the size of the tile. Dask was also five to six times faster here, mostly because eager NumPy reads and decompresses the files on one thread while dask reads 8 chunks at once. ## Step 6: time chunk sizes and thread counts Which chunk size and how many threads? The only reliable answer is to measure. The script runs the 8-date median for five chunk sizes and up to five thread counts, each in a fresh process, then repeats the default with rioxarray’s read lock switched on and once on a distributed `LocalCluster`: `dask_raster.py` ```python for c in [256, 512, 1024, 2048, 4096]: for w in [1, 2, 4, 8, 16]: if w == 1 and c not in (512, 2048): continue # single-threaded runs are slow; two are enough to show scaling r = run_bench("dask", 8, chunk=c, workers=w) lock_run = run_bench("dask", 8, chunk=2048, workers=8, lock=True) dist_run = run_bench("distributed", 8, chunk=2048, workers=2) ``` Inside each process the benchmark is just the pipeline from Steps 2 and 3 with a different scheduler setting: `raster_pipeline.py` ```python red, nir = open_stack(dates, chunk, lock=lock) with dask.config.set(scheduler="threads", num_workers=workers): result = median_ndvi(red, nir).compute().values ``` `Output — Step 6` ```text 256 px chunks, 2 threads: 125.0 s, peak 4.94 GB 256 px chunks, 4 threads: 88.8 s, peak 4.92 GB 256 px chunks, 8 threads: 60.3 s, peak 5.04 GB 256 px chunks, 16 threads: 55.4 s, peak 5.23 GB 512 px chunks, 1 threads: 136.5 s, peak 4.58 GB 512 px chunks, 2 threads: 71.6 s, peak 4.57 GB 512 px chunks, 4 threads: 37.2 s, peak 4.69 GB 512 px chunks, 8 threads: 20.4 s, peak 4.80 GB 512 px chunks, 16 threads: 18.5 s, peak 4.88 GB 1024 px chunks, 2 threads: 60.0 s, peak 4.45 GB 1024 px chunks, 4 threads: 34.0 s, peak 4.47 GB 1024 px chunks, 8 threads: 18.9 s, peak 4.70 GB 1024 px chunks, 16 threads: 12.8 s, peak 5.39 GB 2048 px chunks, 1 threads: 122.3 s, peak 5.01 GB 2048 px chunks, 2 threads: 67.1 s, peak 5.07 GB 2048 px chunks, 4 threads: 35.9 s, peak 5.33 GB 2048 px chunks, 8 threads: 21.4 s, peak 6.54 GB 2048 px chunks, 16 threads: 14.2 s, peak 9.00 GB 4096 px chunks, 2 threads: 67.9 s, peak 8.22 GB 4096 px chunks, 4 threads: 40.0 s, peak 10.44 GB 4096 px chunks, 8 threads: 23.1 s, peak 13.94 GB 4096 px chunks, 16 threads: 18.0 s, peak 20.92 GB 2048 px chunks, 8 threads, lock=True: 108.2 s, peak 5.03 GB LocalCluster, 4 processes x 2 threads, 2048 px chunks: 21.2 s (memory is in the worker processes) fastest: 1024 px chunks with 16 threads, 12.8 s; 2048 px: 1 thread 122 s -> 8 threads 21 s (5.7x) wrote ndvi_median_2025_30UXC.tif: 370 MB on disk in 21.1 s; memory rose by at most 2.49 GB for a 482 MB result check: 10980 x 10980 float32, EPSG:32630, deflate, 484 blocks of (512, 512) ``` ![Left: wall-clock time against threads from 1 to 16 on a logarithmic axis, one line per chunk size. All lines fall as threads are added; the 256 px line is highest, from 125 s on 2 threads to 55 s on 16; the 512, 1024, 2048 and 4096 px lines fall together from 120 to 140 s on 1 thread to 13 to 18 s on 16, with 1024 px lowest. A black cross for 2048 px with lock=True sits at 108 s on 8 threads. Right: peak memory against threads. The 256, 512 and 1024 px lines stay flat between 4.5 and 5.4 GB; the 2048 px line rises from 5 to 9 GB and the 4096 px line from 8.2 to 20.9 GB as threads increase.](https://www.technicalgis.com/_astro/05-timing-chunks-threads.BFDb3wZL_1c5PaW.webp) _Time and peak memory for the median NDVI of eight full scenes, by chunk size and thread count. Very small chunks are slow at any thread count; large chunks multiply memory by the number of threads. The cross is the 2048 px run with rioxarray's read lock switched on._ What the grid shows: - **Threads help almost linearly up to 8.** With 2048 px chunks, 1 thread took 122 s and 8 threads 21 s, a speed-up of 5.7?; 16 threads took 14 s (8.6?). The gain flattens because the workstation has 16 physical cores, the reads compete for the disk cache, and other jobs were running. - **Tiny chunks are slow.** 256 px chunks make a graph of 314,346 tasks, and the scheduler, not the CPU, becomes the bottleneck: 55 s on 16 threads, four times slower than 1024 px. If the task stream from Step 4 shows lots of thin bars and gaps, your chunks are too small. - **Big chunks are hungry.** Peak memory is roughly a few 8-date blocks per thread. At 4096 px each block is 537 MB of float32 and its NDVI temporaries, so 16 threads needed 20.9 GB, more than eager NumPy needed for four dates. At 1024 px, memory stayed between 4.5 and 5.4 GB whatever the thread count. - **The sweet spot is broad.** 512, 1024 and 2048 px chunks were within a few seconds of each other at each thread count. 1024 px on 16 threads was fastest (12.8 s); 2048 px on 8 threads (21 s, 6.5 GB) is a sensible default on an ordinary laptop. - **Check the read lock.** By default (`lock=None` or `True`), rioxarray puts one shared lock around every read, so only one thread can read and decompress a block at a time, from any file. That is safe for any GDAL driver, but here it made the run five times slower (108 s against 21 s) and left the other seven threads waiting. Our `open_stack` passes `lock=False`, which opens a separate file handle for each read and is safe for local GeoTIFFs and COGs. - **A LocalCluster was no faster here.** Four worker processes with two threads each took 21.2 s, the same as eight threads in one process. NumPy and GDAL release the GIL for the heavy work, so threads are enough. The distributed scheduler earns its start-up cost when tasks hold the GIL (pure-Python code), when you need its dashboard or memory limits, or when the work spreads over several machines. The floor of about 4.5 GB, even on one thread, is mostly not dask. It includes the 482 MB result, Python and the libraries, and GDAL’s block cache, which by default may grow to 5% of physical memory (3.4 GB on this machine). Set the `GDAL_CACHEMAX` environment variable (for example to `512` MB) if you need a smaller footprint. ## Step 7: write the result to disk chunk by chunk For a result that is too big for memory, or one you simply want as a file, skip `compute()` and let rioxarray write it. With `windowed=True`, each chunk is written to its own window of the GeoTIFF as soon as it is computed. GDAL cannot write to one file from several threads at once, so pass a lock; the computation still runs in parallel, and only the writes take turns: `dask_raster.py` ```python med = median_ndvi(*open_stack(DATES, 2048)).rio.write_nodata(np.nan, encoded=False) with dask.config.set(scheduler="threads", num_workers=8): med.rio.to_raster(dest, tiled=True, blockxsize=512, blockysize=512, compress="deflate", predictor=3, windowed=True, lock=threading.Lock()) ``` `Output — Step 7` ```text wrote ndvi_median_2025_30UXC.tif: 370 MB on disk in 21.1 s; memory rose by at most 2.49 GB for a 482 MB result check: 10980 x 10980 float32, EPSG:32630, deflate, 484 blocks of (512, 512) ``` The write took 21 seconds, about as long as computing the result in Step 4, and the process memory rose by at most 2.5 GB, against 6.6 GB when the result was gathered into memory. Each finished chunk goes to disk and is released, so this pattern works for results far larger than RAM. The file is a 370 MB deflate-compressed GeoTIFF with `predictor=3` (the floating-point predictor), in 484 internal blocks of 512 px, the same block size as the inputs, so the next dask job can read it efficiently. Two details are easy to miss. `rio.write_nodata(np.nan, encoded=False)` marks NaN as no data in the output, and `out.attrs = {}` in `median_ndvi` drops the input’s uint16 `_FillValue`, which would otherwise be written as a nodata value of 0 into a float file. To make the result a full Cloud Optimised GeoTIFF with overviews, write it with `driver="COG"` instead, which builds the overviews at the end, or add them afterwards with `gdaladdo`. ![Left: the median NDVI of the whole 110 by 110 km tile for March to November 2025, mostly green, blue and yellow farmland and woodland, with the red and orange urban area of north-west and central London and the Thames in the bottom right, and dark red reservoirs near Heathrow; a small black rectangle marks a zoom area in the centre. Right: the 12 by 9 km zoom at full resolution around Aylesbury, showing individual fields in shades of blue, green and orange, the town in red, and a railway running diagonally.](https://www.technicalgis.com/_astro/06-median-ndvi-result.BoGC36Fh_2c4mYF.webp) _The median NDVI of eight scenes from March to November 2025, read back from the GeoTIFF written in Step 7. The whole 110 km tile (left) is read at one tenth of full resolution; the zoom around Aylesbury (right) is one 1200 ? 900 px window at full 10 m resolution, which rasterio reads from just the blocks it needs._ ## Pitfalls and limitations - **Chunks that split the file’s tiles.** A chunk size that is not a multiple of the internal block size makes neighbouring chunks decompress the same blocks. Check `src.block_shapes` with rasterio, or `gdalinfo`, before choosing. - **Too many tasks.** Every task costs the scheduler time. If `len(arr.__dask_graph__())` runs into hundreds of thousands, use bigger chunks. - **Reductions across chunks.** A median, percentile or sort along an axis needs that axis in one chunk; dask raises an error if it is not. Rechunk once, as in Step 3, rather than repeatedly. - **`compute()` loads the result.** Dask keeps the inputs out of memory, not the output. If the result itself is too big, write it with `to_raster`, `to_zarr` or `to_netcdf` instead of computing it. - **Remote files.** Reading COGs over HTTP from the Planetary Computer works with the same code, but each chunk then costs a network round trip, and the best chunk size and thread count change. The signed URLs also expire after about an hour, so long jobs need fresh ones. - **Memory numbers on Windows.** We measured the resident set size of the process with psutil. The operating system and the memory allocator do not hand memory back immediately, so the memory curves show a high-water mark, not what dask is holding at each moment. - **Benchmarks are specific.** These timings are for one workstation, one calculation and files on a local SSD. A different reduction, a network file system or a machine with less memory will give different best settings, so time your own pipeline on a sample before you run it on a large area. ## Frequently asked questions ### How do I process a large raster in Python without running out of memory? Open it with `rioxarray.open_rasterio(path, chunks=...)` or `xarray.open_dataset(..., chunks=...)` so that it becomes a lazy dask array, write your calculation as normal, and either compute a result that fits in memory or write the result to disk with `rio.to_raster(..., windowed=True, lock=...)`. Dask then streams the raster through memory in chunks. In this tutorial, eight full Sentinel-2 scenes went through in 6.4 GB, while eager NumPy needed 12.4 GB for four. ### What chunk size should I use with dask and rioxarray? Use a multiple of the file’s internal block size (512 px for Sentinel-2 COGs, often 256 or 512 for other GeoTIFFs), and aim for chunks of roughly 10–100 MB in the data type you compute in. For Sentinel-2 that means 1024 or 2048 px. In our tests 512, 1024 and 2048 px chunks were all within a few seconds of each other, 256 px chunks were four times slower, and 4096 px chunks needed up to 21 GB on 16 threads. Keep any axis you reduce over (such as time for a median) in a single chunk. ### Is dask faster than NumPy for raster processing? For data that fits in memory, NumPy on one core can be just as fast for simple element-wise maths, because NumPy has no scheduler overhead. Dask wins when the work is large enough to keep several threads busy, when reading and decompressing the files is a big part of the cost, or when the data does not fit in memory at all. Here dask was 5 to 6 times faster than eager NumPy, mostly because it read eight chunks in parallel. ### Should I use the threaded scheduler or a dask distributed LocalCluster? For one machine and raster maths built on NumPy and GDAL, start with the default threaded scheduler: those libraries release the GIL, so threads run in parallel, share memory and have no start-up cost. In our test a `LocalCluster` of 4 processes ? 2 threads took 21.2 s, the same as 8 threads. Use `distributed` when your tasks run pure-Python code that holds the GIL, when you want its live dashboard and per-worker memory limits, or when you scale out to a cluster. On Windows and macOS, create the cluster inside `if __name__ == "__main__":`, because worker processes are spawned by re-importing your script. ### Can I use dask with rasterio directly, without xarray? Yes, with `dask.array.from_delayed` or `da.map_blocks` and your own windowed reads, but you then manage coordinates, CRS and nodata yourself. rioxarray does that for you and keeps the georeferencing through every step, which is why it is the usual choice. For writing, `rioxarray` handles the windowed writes and the lock. ### How do I see what dask is doing? Wrap `compute()` in `dask.diagnostics.ProgressBar()` for a progress bar and `Profiler()` for a task stream, as in Step 4, or start a `distributed.Client` and open its dashboard (it needs Bokeh) to watch tasks, memory and workers live in the browser. `arr.data.visualize()` draws the task graph if Graphviz is installed. ## Next steps - Summarise the median NDVI by area: [zonal statistics in Python](https://www.technicalgis.com/tutorials/zonal-statistics-python/) works with chunked rasters too. - Build a per-date series instead of one median: [NDVI time series in Python with Sentinel-2 and odc-stac](https://www.technicalgis.com/tutorials/ndvi-time-series-odc-stac/) uses odc-stac, which returns the same kind of dask-backed `DataArray`. - Run a neural network over a whole tile: [deep learning inference on large rasters](https://www.technicalgis.com/tutorials/deep-learning-inference-large-rasters/) tiles a large image for a model, the same idea as chunks. - Apply the same chunking to elevation: [raster terrain analysis with rasterio](https://www.technicalgis.com/tutorials/raster-terrain-analysis-rasterio/) computes slope and aspect, which need overlapping chunks (`map_overlap` in dask). - Start with STAC and COGs: [Sentinel-2 imagery of Lake Mead with STAC](https://www.technicalgis.com/tutorials/sentinel-2-stac-lake-mead/). - Turn composites into an animation: [satellite timelapse in Python](https://www.technicalgis.com/tutorials/satellite-timelapse-python/) builds one median composite per year. ## Data and licences - **Sentinel-2 L2A**: contains modified Copernicus Sentinel data 2025, processed by ESA. Bands B04 and B08 of tile 30UXC for eight dates, downloaded as Cloud Optimised GeoTIFFs from the [Microsoft Planetary Computer](https://planetarycomputer.microsoft.com/dataset/sentinel-2-l2a) on 3 November 2026. Copernicus Sentinel data are free and open under the Copernicus data licence. - Code: [`dask_raster.py` and `raster_pipeline.py`](https://www.technicalgis.com/code/parallel-raster-processing-dask-python/), tested with Python 3.12.0, dask 2026.8.0, distributed 2026.8.0, xarray 2026.9.0, rioxarray 0.23.0, rasterio 1.5.2, NumPy 2.5.4, pandas 3.1.0, matplotlib 3.11.2, networkx 3.7, psutil 7.2.2, pystac-client 0.9.1 and planetary-computer 1.0.0. ## Key terms - **Chunk**: One block of a dask array, such as 2048 × 2048 pixels of one band on one date. Each chunk is an ordinary NumPy array, and dask runs one task per chunk per step. - **Lazy evaluation**: Recording a calculation without running it. Dask and xarray build a task graph as you write the code, and only read data and compute when you ask for a result. - **Task graph**: The network of small tasks that make up a dask calculation, with an edge wherever one task needs the output of another. Dask's scheduler walks the graph and runs independent tasks in parallel. - **Scheduler**: The part of dask that decides which task runs where and when. The threaded scheduler uses threads in one process; the distributed scheduler (LocalCluster or a real cluster) uses separate worker processes. - **Internal tiles (blocks)**: The fixed-size blocks in which a tiled GeoTIFF or COG stores its pixels. Reading part of a block means reading and decompressing the whole block, so dask chunks should be whole multiples of it. - **Out-of-core**: Processing data that is larger than memory by streaming it through in pieces, rather than loading it all at once. ## Data for this tutorial - [Microsoft Planetary Computer](https://www.technicalgis.com/data-sources/#src-microsoft-planetary-computer) · Sentinel-2 L2A as Cloud Optimised GeoTIFFs via STAC - [Copernicus Data Space Ecosystem](https://www.technicalgis.com/data-sources/#src-copernicus-data-space-ecosystem) · The official Sentinel-2 archive Licences and formats for each source are listed in the [GIS data sources catalogue](https://www.technicalgis.com/data-sources/). --- Source: https://www.technicalgis.com/tutorials/parallel-raster-processing-dask-python/