How to Choose Chunk and Tile Sizes That Actually Help
On this page
Problem statement
Chunk and tile size is the parameter people copy from an example and never revisit, and it decides more about performance than any other choice in a cloud-native pipeline.
Window reads from the same 4096 × 4096 COG with 512-pixel tiles (a Sentinel-2 near-infrared crop, deflate-compressed, served over HTTP ranges; Example 4 measures them) show why:
window requests bytes fetched bytes per pixel
1 px at 0 2 426.0 kB 425,984.0
256 px at 0 2 426.0 kB 6.5
512 px at 0 2 426.0 kB 1.6
512 px at 1792 3 1,472.8 kB 5.6
2048 px at 0 5 6,187.3 kB 1.5
Reading one pixel and reading 262,144 pixels cost exactly the same, because both touch one tile. The tile is the unit of transfer, and every sizing decision follows from that.
Quick answer
Work backwards from the access pattern:
def suggest_tile_size(typical_window_px, dtype_bytes=4, bands=1,
target_bytes=(1e6, 10e6)):
"""A tile that covers the typical read without much waste."""
for candidate in (128, 256, 512, 1024, 2048):
tiles = max(1, -(-typical_window_px // candidate)) ** 2
fetched = tiles * candidate ** 2 * dtype_bytes * bands
useful = typical_window_px ** 2 * dtype_bytes * bands
chunk_bytes = candidate ** 2 * dtype_bytes * bands
ok = target_bytes[0] <= chunk_bytes <= target_bytes[1]
print(f" {candidate:5d} px {tiles:4d} tiles "
f"{chunk_bytes / 1e6:6.2f} MB/tile "
f"{fetched / useful:5.1f}x amplification "
f"{'ok' if ok else ''}")
suggest_tile_size(512) # a 512 × 512 read of float32 128 px 16 tiles 0.07 MB/tile 1.0x amplification
256 px 4 tiles 0.26 MB/tile 1.0x amplification
512 px 1 tiles 1.05 MB/tile 1.0x amplification ok
1024 px 1 tiles 4.19 MB/tile 4.0x amplification ok
2048 px 1 tiles 16.78 MB/tile 16.0x amplification
Two constraints: the chunk should be a few megabytes, and it should not be much larger than the typical read.
The typical read, then size, alignment and shape
1. Establish the typical read
Not the largest, not the smallest — the one that happens most. A tile server reads 256 × 256; an analysis pipeline reads a study area; a time-series extraction reads a pixel.
Those want different layouts, which is why "what chunk size should I use" has no answer without the access pattern.
2. Apply the size constraint
A few megabytes per chunk, uncompressed:
- Under about 1 MB: per-request overhead dominates. Over HTTP that is a round trip per chunk; on object storage it is also a billed request.
- Over about 100 MB: any partial read wastes most of what it fetched, and memory per task grows.
For a float32 single band, 512 × 512 is 1.05 MB and 1024 × 1024 is 4.19 MB. Both are reasonable; 128 × 128 at 0.07 MB is not.
3. Apply the alignment constraint
An unaligned window straddles more tiles than an aligned one. A 512-pixel window at offset 0 touches one 512-tile; the same window at offset 256 touches four.
The measurement shows it: the 512-pixel read at offset 1792 — three and a half tiles in, so it straddles four tiles — fetched 1,472.8 kB, against 426.0 kB for the same size of read wholly inside one tile.
Where you control the reader, align the windows. Where you do not, choose tiles somewhat larger than the typical window so a straddling read touches two rather than four.
4. For n-dimensional data, choose the shape as well as the size
A (time, band, y, x) cube of 120 dates × 4 bands × 512 × 512 pixels has almost the same chunk size — 4.19 MB and 4.06 MB of float32 — for (1, 4, 512, 512) and (120, 4, 46, 46), and they behave completely differently (Example 1 counts them):
| chunk | one date's map | one pixel's history |
|---|---|---|
(1, 4, 512, 512) |
1 read, 4.2 MB | 120 reads, 503.3 MB |
(120, 4, 46, 46) |
144 reads, 585.0 MB | 1 read, 4.1 MB |
If both patterns matter, either store two copies or compromise: weighting maps at half the reads, pixel histories at 0.3 and small-window histories at 0.2, Example 3 picks (30, 4, 128, 128).
5. Measure the amplification, not the time
bytes fetched / bytes used
Time varies with network, cache and load. Amplification is a property of the layout and the query, and it is the number to design against. Above about 10 for your typical read, the layout is wrong.
Measuring read amplification on a Sentinel-2 crop
The data used on this page
A 4096 × 4096 crop of the 10 m near-infrared band of Sentinel-2 scene S2A_30UVD_20250922_0_L2A from Earth Search, written as a plain uncompressed GeoTIFF, nir.tif; Example 4 rewrites it as a COG with 512-pixel tiles. To count what a read costs, the setup also starts a small HTTP server in a separate process that answers byte-range requests for local files and logs the size of every response; cost(fn, path) runs fn against a fresh URL for path and returns the number of requests and the bytes sent. The server has to be its own process: GDAL holds Python's global interpreter lock while it waits for a response, so a server thread inside the same interpreter never gets to answer.
import atexit
import itertools
import os
import subprocess
import sys
import pystac_client
import rasterio
from rasterio.windows import Window
item = next(pystac_client.Client.open("https://earth-search.aws.element84.com/v1")
.search(ids=["S2A_30UVD_20250922_0_L2A"], collections=["sentinel-2-l2a"]).items())
if not os.path.exists("nir.tif"): # a 4096 x 4096 crop of the 10 m near-infrared band
with rasterio.open(item.assets["nir"].href) as src:
win = Window(3000, 3000, 4096, 4096)
profile = {"driver": "GTiff", "dtype": "uint16", "count": 1, "width": 4096, "height": 4096,
"crs": src.crs, "transform": src.window_transform(win)}
with rasterio.open("nir.tif", "w", **profile) as dst:
dst.write(src.read(1, window=win), 1)
SERVER = r"""
import http.server, os, sys
class Counting(http.server.SimpleHTTPRequestHandler):
def do_GET(self):
path = self.translate_path(self.path.split("?")[0])
size = os.path.getsize(path)
first, last = 0, size - 1
if "Range" in self.headers:
a, b = self.headers["Range"].split("=")[1].split("-")
first, last = int(a), min(int(b) if b else size - 1, size - 1)
self.send_response(206)
self.send_header("Content-Range", f"bytes {first}-{last}/{size}")
else:
self.send_response(200)
self.send_header("Content-Length", str(last - first + 1))
self.send_header("Accept-Ranges", "bytes")
self.end_headers()
with open(path, "rb") as f:
f.seek(first)
self.wfile.write(f.read(last - first + 1))
with open("requests.log", "a") as log: # one line per GET: bytes sent
log.write(f"{last - first + 1}\n")
def log_message(self, *args):
pass
server = http.server.ThreadingHTTPServer(("127.0.0.1", 0), Counting)
print(server.server_port, flush=True)
server.serve_forever()
"""
# Its own process: GDAL holds the GIL while it waits for a response, so a server
# thread inside this interpreter would never get to answer.
proc = subprocess.Popen([sys.executable, "-c", SERVER], stdout=subprocess.PIPE, stderr=subprocess.DEVNULL, text=True)
atexit.register(proc.terminate)
BASE = f"http://127.0.0.1:{proc.stdout.readline().strip()}"
_fresh = itertools.count()
def served(path):
"""A URL for a local file that no earlier read has cached."""
return f"{BASE}/{path}?read={next(_fresh)}"
def cost(fn, path):
"""Requests and bytes the server sent while fn(url) ran."""
open("requests.log", "w").close()
with rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR"):
fn(served(path))
sizes = [int(line) for line in open("requests.log")]
return len(sizes), sum(sizes)
print(rasterio.open("nir.tif").shape, f"{os.path.getsize('nir.tif'):,} bytes uncompressed")(4096, 4096) 33,579,368 bytes uncompressed
Example 1 — measuring amplification for a real access pattern
import numpy as np
def amplification(shape, chunks, queries, dtype_bytes=4):
"""Bytes fetched divided by bytes used, per query shape."""
print(f" array {shape}, chunks {chunks}, "
f"{np.prod(chunks) * dtype_bytes / 1e6:.2f} MB per chunk")
results = {}
for name, query in queries.items():
touched = 1
for dim, (size, chunk) in enumerate(zip(shape, chunks)):
span = query.get(dim, size)
touched *= int(np.ceil(span / chunk))
used = np.prod([query.get(d, s) for d, s in enumerate(shape)])
fetched = touched * np.prod(chunks)
results[name] = fetched / max(used, 1)
print(f" {name:26} {touched:6,} chunks "
f"{fetched * dtype_bytes / 1e6:9.1f} MB fetched "
f"{results[name]:8.1f}x")
return results
cube = (120, 4, 512, 512) # 120 dates, 4 bands, 512 x 512 pixels
queries = {"one date's map": {0: 1}, "one pixel's history": {2: 1, 3: 1}}
for chunks in ((1, 4, 512, 512), (120, 4, 46, 46)):
amplification(cube, chunks, queries) array (120, 4, 512, 512), chunks (1, 4, 512, 512), 4.19 MB per chunk
one date's map 1 chunks 4.2 MB fetched 1.0x
one pixel's history 120 chunks 503.3 MB fetched 262144.0x
array (120, 4, 512, 512), chunks (120, 4, 46, 46), 4.06 MB per chunk
one date's map 144 chunks 585.0 MB fetched 139.5x
one pixel's history 1 chunks 4.1 MB fetched 2116.0x
Running this for the two or three queries you actually make, over three or four candidate chunk shapes, produces a table that decides the question in a minute.
Example 2 — aligning windows to the tile grid
import numpy as np
from rasterio.windows import Window
def align_window(window, block_w, block_h):
"""Expand a window outward to whole tiles."""
col0 = int(window.col_off // block_w) * block_w
row0 = int(window.row_off // block_h) * block_h
col1 = int(np.ceil((window.col_off + window.width) / block_w)) * block_w
row1 = int(np.ceil((window.row_off + window.height) / block_h)) * block_h
return Window(col0, row0, col1 - col0, row1 - row0)
def read_aligned(ds, window, band=1):
"""Read whole tiles, then slice to the requested window."""
block_h, block_w = ds.block_shapes[band - 1]
outer = align_window(window, block_w, block_h)
data = ds.read(band, window=outer)
r0 = int(window.row_off - outer.row_off)
c0 = int(window.col_off - outer.col_off)
print(f" requested {int(window.width)}x{int(window.height)}, "
f"read {int(outer.width)}x{int(outer.height)} "
f"({outer.width * outer.height / (window.width * window.height):.1f}x)")
return data[r0:r0 + int(window.height), c0:c0 + int(window.width)]
with rasterio.open("nir.tif") as src:
data, profile = src.read(1), src.profile
profile.update(driver="COG", blocksize=512, compress="deflate", predictor=2)
with rasterio.open("cog512.tif", "w", **profile) as dst:
dst.write(data, 1)
window = Window(1792, 1792, 300, 300) # an unaligned 300 x 300 request
with rasterio.open("cog512.tif") as ds:
part = read_aligned(ds, window)
assert np.array_equal(part, ds.read(1, window=window))
def plain(url):
with rasterio.open(url) as ds:
ds.read(1, window=window)
def aligned(url):
with rasterio.open(url) as ds:
read_aligned(ds, window)
for name, fn in (("as requested", plain), ("aligned", aligned)):
n, sent = cost(fn, "cog512.tif")
print(f" {name:13} {n} requests {sent / 1000:9,.1f} kB") requested 300x300, read 1024x1024 (11.7x)
as requested 3 requests 1,472.8 kB
requested 300x300, read 1024x1024 (11.7x)
aligned 3 requests 1,472.8 kB
The two reads cost exactly the same — 3 requests and 1,472.8 kB — because GDAL fetches whole tiles whichever window you ask for: the 300 × 300 request straddles four 512-pixel tiles, so 1024 × 1024 pixels cross the wire either way. Aligning in the reader does not save bytes; what it does is make the cost visible, since the printed ratio (11.7x here) is the amplification you are paying. The saving comes from aligning the requests — choosing windows that start on a tile boundary.
Example 3 — picking a chunking for a cube from its queries
import numpy as np
def choose_cube_chunks(shape, dtype_bytes=4, target_mb=(2, 20),
weights=None):
"""Score candidate chunk shapes against weighted query patterns."""
n_time, n_band, height, width = shape
weights = weights or {"map": 0.5, "pixel_series": 0.3, "window_series": 0.2}
candidates = []
for t in (1, 5, 30, n_time):
for s in (64, 128, 256, min(height, width)):
chunk = (t, n_band, min(s, height), min(s, width))
mb = np.prod(chunk) * dtype_bytes / 1e6
if not target_mb[0] <= mb <= target_mb[1]:
continue
cost = (
weights["map"] * np.ceil(n_time / t) ** 0 *
np.ceil(height / chunk[2]) * np.ceil(width / chunk[3]) +
weights["pixel_series"] * np.ceil(n_time / t) +
weights["window_series"] * np.ceil(n_time / t) *
np.ceil(64 / chunk[2]) * np.ceil(64 / chunk[3])
)
candidates.append((cost, chunk, mb))
candidates.sort()
for cost, chunk, mb in candidates[:6]:
print(f" {str(chunk):28} {mb:6.2f} MB weighted cost {cost:8.1f}")
return candidates[0][1] if candidates else None
print(choose_cube_chunks((120, 4, 512, 512))) (30, 4, 128, 128) 7.86 MB weighted cost 10.0
(5, 4, 256, 256) 5.24 MB weighted cost 14.0
(120, 4, 64, 64) 7.86 MB weighted cost 32.5
(1, 4, 512, 512) 4.19 MB weighted cost 60.5
(30, 4, 128, 128)
Weighting the query patterns by how often they occur is what turns an argument into a calculation. If maps are half your reads and pixel series a third, the optimum is not the shape that is best at either.
Example 4 — what one window costs, measured
print(f"cog512.tif: {os.path.getsize('cog512.tif'):,} bytes, 512-pixel tiles")
def read_window(size, offset):
def fn(url):
with rasterio.open(url) as ds:
ds.read(1, window=Window(offset, offset, size, size))
return fn
print(f"{'window':>16} {'requests':>9} {'bytes fetched':>14} {'bytes per pixel':>16}")
for size, offset in ((1, 0), (256, 0), (512, 0), (512, 1792), (2048, 0)):
n, sent = cost(read_window(size, offset), "cog512.tif")
print(f"{size:5} px at {offset:4} {n:9} {sent / 1000:11,.1f} kB {sent / size ** 2:16,.1f}")cog512.tif: 32,989,693 bytes, 512-pixel tiles
window requests bytes fetched bytes per pixel
1 px at 0 2 426.0 kB 425,984.0
256 px at 0 2 426.0 kB 6.5
512 px at 0 2 426.0 kB 1.6
512 px at 1792 3 1,472.8 kB 5.6
2048 px at 0 5 6,187.3 kB 1.5
Why both very small and very large chunks fail
Why the tile is the unit of transfer
Chunks and tiles are compressed independently, so a reader cannot decompress part of one. The minimum transfer for any element is its whole containing chunk.
That is why the measurement shows 1 pixel, 256 × 256 and 512 × 512 pixels costing the same 426.0 kB. It is also why per-pixel access patterns are pathological: a thousand scattered pixels can cost a thousand chunks.
Every sizing rule is a consequence: make chunks small enough that a targeted read is targeted, and large enough that the per-chunk overhead is amortised.
Why very small chunks fail
Each chunk costs a request, a compression header and a metadata entry. Over HTTP the request is a round trip — 100 ms against a remote store.
Reading a 4096 × 4096 raster in 128-pixel tiles is 1,024 requests. At 100 ms each, serially, that is 102 seconds for 67 MB of data. The same raster in 512-pixel tiles is 64 requests.
Object stores also bill per request, so tiny chunks cost money as well as latency.
Why very large chunks fail
A 4096 × 4096 single-chunk raster requires reading all 67 MB to obtain any pixel. That is the pathological case that tiling was invented to avoid.
Less obviously, large chunks hurt parallelism. A task processes a chunk, so chunk size sets the granularity of the work. Eight workers and four chunks means half the workers are idle.
Why alignment is free performance
An unaligned window straddles chunk boundaries in both dimensions, so it can touch four chunks where an aligned one touches one — a factor of four in bytes for the same output.
The measured example: a 512-pixel read at offset 1792 fetched 1,472.8 kB, against 426.0 kB for a read inside one tile.
Where you generate the windows — a tiling job, a batch extraction — aligning them is a few lines and often the single largest available improvement.
Whole-chunk reads, parallel granularity and multiples of 16
- A 1-pixel read costs a whole chunk. Read blocks, not points.
- A few megabytes per chunk is the broad optimum.
- Align windows to the chunk grid where you control them.
- Chunk shape matters as much as size for n-dimensional data.
- Amplification, not time, is the metric. Time varies with the network.
- Chunk size sets parallel granularity — too few chunks starves workers.
- Tile dimensions must be multiples of 16 in GeoTIFF.
- Compressed chunk sizes vary; design against the uncompressed size.
Internal links
- Cloud-optimised GeoTIFF explained — where tile size is set
- Chunked arrays and Zarr explained — chunk shape for n-dimensional data
- How to write a cloud-optimised GeoTIFF in Python — applying the choice
- How to read a COG from a URL without downloading the whole file — aligning reads
- Reading a COG from a URL is slow or downloads everything — when the layout is wrong
- Lazy loading explained: why xarray reads nothing until you ask — chunk size and task count
- How to scale a GeoPandas job with dask-geopandas — partition sizing for tables
- How to process a very large GeoPackage in chunks with Python — chunking without a chunked format
FAQ
What chunk size should I use?
A few megabytes uncompressed — 512 × 512 for a single-band float32 raster is 1.05 MB and a good default.
Why does reading one pixel cost so much?
Because chunks are compressed independently, so the minimum transfer is one whole chunk. One pixel and a 512 × 512 window both fetched 426.0 kB from the 512-tile COG in Example 4.
Does chunk shape matter as well as size?
For n-dimensional data, enormously. Time-chunked and space-chunked cubes of identical chunk size differ by two orders of magnitude on the same query.
What is read amplification?
Bytes fetched divided by bytes used. It is a property of the layout and the query, unlike time, and it is the right number to design against.
Should I align my read windows to the tile grid?
Yes where you can. An unaligned window can touch four tiles where an aligned one touches one.
Are smaller chunks safer?
No. Below about a megabyte, per-request overhead and task count dominate. A 4096-pixel raster in 128-pixel tiles is 1,024 requests.
How does chunk size affect parallelism?
Chunks are the unit of work. Too few chunks and workers sit idle; too many and the scheduler spends its time on bookkeeping.