Reading a COG from a URL Is Slow or Downloads Everything
Problem statement
A remote read that should transfer a megabyte transfers the whole file, or takes seconds when it should take milliseconds. There are five causes, and they are distinguishable.
Measured against a byte-counting server on a 4096 × 4096 raster:
whole scene at 1/16 resolution
file with overviews 1 request 0.115 MB
file without overviews 3 requests 24.672 MB <- the whole file
four scattered windows
default GDAL settings 10 requests 4.610 MB
GDAL_DISABLE_READDIR_ON_OPEN 6 requests 4.155 MB
The first is a property of the file. The second is a property of your configuration. Diagnosing which you have is the whole job.
Quick answer
import rasterio
with rasterio.open(url) as ds:
print(f"blocks {ds.block_shapes[0]}")
print(f"overviews {ds.overviews(1) or 'NONE'}")
print(f"compress {ds.profile.get('compress')}")blocks (1, 4096)
overviews NONE
compress deflate
A block height of 1 means the file is striped, not tiled. An empty overview list means every downsampled read costs full resolution. Either makes the file behave like an ordinary GeoTIFF, and no client setting fixes it.
Step-by-step solution
1. Cause one: the file has no overviews
The largest single effect. Without overviews, GDAL produces a downsampled image by reading every full-resolution pixel and averaging.
whole scene at 1/16
with overviews 0.115 MB
without 24.672 MB
There is no client-side fix. Either rewrite the file with overviews, or accept full-resolution reads.
Where you cannot rewrite the original, an external .ovr sidecar alongside it works — but only if the server hosts it and CPL_VSIL_CURL_ALLOWED_EXTENSIONS permits GDAL to look.
2. Cause two: the file is striped
one 512 x 512 window
tiled 1.557 MB
striped 3.244 MB
A striped file stores whole rows, so a narrow window reads full-width strips. The waste grows with the raster's width: on a 10,980-pixel-wide Sentinel-2 band, a 512-pixel window reads 21 times more pixels than it needs.
Again, no client fix.
3. Cause three: sidecar probing
default settings 10 requests
GDAL_DISABLE_READDIR_ON_OPEN 6 requests
Four requests were GDAL looking for .aux.xml, .ovr, .msk and world files that do not exist. Each 404 is a full round trip.
rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR",
CPL_VSIL_CURL_ALLOWED_EXTENSIONS=".tif")This is the highest-value client-side setting, and against a high-latency store it often dominates everything else.
4. Cause four: reading full resolution and downsampling in NumPy
data = ds.read(1)[::16, ::16] # reads everything
data = ds.read(1, out_shape=(h // 16, w // 16),
resampling=Resampling.average) # reads an overviewThe first line is the single most common way to accidentally download a whole file. It looks like a downsampled read and is not.
5. Cause five: an access pattern that is not tile-shaped
1 px 470.8 kB 470,793 bytes/px
256 px 470.8 kB 7.2 bytes/px
512 px 1,585.3 kB 6.0 bytes/px
One pixel costs a whole tile. Sampling a raster at a thousand scattered points can cost a thousand tiles; reading the block containing them costs a handful.
Code examples
Example 1 — a diagnostic that names the cause
import time
import rasterio
from rasterio.enums import Resampling
def diagnose_remote(url, settings=None):
"""Identify which of the five causes applies."""
env = settings or dict(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR",
CPL_VSIL_CURL_ALLOWED_EXTENSIONS=".tif")
problems = []
start = time.time()
with rasterio.Env(**env):
with rasterio.open(url) as ds:
open_ms = (time.time() - start) * 1000
block_h, block_w = ds.block_shapes[0]
overviews = ds.overviews(1)
print(f" {ds.width} x {ds.height} {ds.dtypes[0]}, "
f"{ds.count} band(s)")
print(f" blocks {(block_h, block_w)}, overviews "
f"{overviews or 'NONE'}, compress "
f"{ds.profile.get('compress')}")
print(f" open took {open_ms:.0f} ms")
if block_h == 1 or block_h == ds.height:
problems.append(
f"striped ({block_h} x {block_w}): a narrow window reads "
f"full-width rows — {ds.width / 512:.0f}x more pixels than "
"a 512-wide window needs")
if not overviews:
problems.append(
"no overviews: any downsampled read costs the full "
"resolution. Rewrite the file or add an .ovr sidecar")
elif max(overviews) * max(block_h, block_w) < max(ds.width, ds.height):
problems.append(
f"coarsest overview is {max(ds.width, ds.height) / max(overviews):.0f} px "
"— a full-extent view still needs several requests")
if open_ms > 500:
problems.append(
"slow open: check GDAL_DISABLE_READDIR_ON_OPEN and the "
"latency to the host")
for p in problems:
print(f" ! {p}")
return problemsExample 2 — proving that the settings are the problem
import time
import rasterio
from rasterio.windows import Window
CONFIGS = {
"defaults": {},
"no readdir": {"GDAL_DISABLE_READDIR_ON_OPEN": "EMPTY_DIR"},
"no readdir + extensions": {
"GDAL_DISABLE_READDIR_ON_OPEN": "EMPTY_DIR",
"CPL_VSIL_CURL_ALLOWED_EXTENSIONS": ".tif",
},
"everything": {
"GDAL_DISABLE_READDIR_ON_OPEN": "EMPTY_DIR",
"CPL_VSIL_CURL_ALLOWED_EXTENSIONS": ".tif",
"VSI_CACHE": "TRUE",
"GDAL_HTTP_MERGE_CONSECUTIVE_RANGES": "YES",
},
}
def compare_settings(url, windows):
"""Time the same reads under several GDAL configurations."""
for name, env in CONFIGS.items():
start = time.time()
with rasterio.Env(**env):
with rasterio.open(url) as ds:
for window in windows:
ds.read(1, window=window)
print(f" {name:26} {(time.time() - start) * 1000:7.0f} ms")Run this in a fresh process per configuration where you can — GDAL caches across Env blocks, so the second configuration measured in one process benefits from the first's cache.
Example 3 — restructuring an access pattern
import numpy as np
import rasterio
from rasterio.windows import Window
def sample_many(url, xs, ys, band=1, max_block_px=50_000_000):
"""One block read instead of a read per point."""
with rasterio.open(url) as ds:
rows, cols = rasterio.transform.rowcol(ds.transform, xs, ys)
rows, cols = np.asarray(rows), np.asarray(cols)
row0, row1 = max(rows.min(), 0), min(rows.max() + 1, ds.height)
col0, col1 = max(cols.min(), 0), min(cols.max() + 1, ds.width)
span = (row1 - row0) * (col1 - col0)
print(f" {len(xs)} points span {col1 - col0} x {row1 - row0} px "
f"({span / 1e6:.1f} Mpx)")
if span > max_block_px:
print(" ! too spread out for one block — cluster the points and "
"read one block per cluster")
return None
block = ds.read(band, window=Window(col0, row0, col1 - col0, row1 - row0))
return block[rows - row0, cols - col0]ds.sample() issues one read per point. For points inside one tile the cache absorbs most of that, and for points spread over a raster it does not — a block read is unambiguous.
Explanation
Why a missing overview is unfixable from the client
Downsampling requires reading the values being averaged. Without overviews those values exist only at full resolution, so producing a 256 × 256 thumbnail from a 4096 × 4096 image genuinely requires reading 16.7 million pixels.
No client setting changes that arithmetic. The saving from overviews is that the averaging was done once, at write time, and stored — 0.115 MB instead of 24.672 MB.
The only remote workaround is an external .ovr file, which contains the same pre-computed pyramid alongside the original. It needs to be hosted next to the file and permitted by CPL_VSIL_CURL_ALLOWED_EXTENSIONS.
Why striping wastes so much
A strip spans the full width of the image. Reading a 512-pixel-wide window from a 4096-wide striped file transfers 4096 pixels per row and uses 512 — eight times the data.
On a Sentinel-2 band at 10,980 pixels wide the ratio is 21. Measured on the 4096-wide test raster, the striped file transferred 3.244 MB against the tiled file's 1.557 MB.
Striping is the default in many writers, which is why "it is a GeoTIFF on S3" is not the same as "it is a COG".
Why the sidecar probes matter so much
GDAL looks for auxiliary files when opening: statistics, external overviews, masks, world files. Locally these are stat calls costing microseconds.
Over HTTP each is a request. Four missing sidecars is four round trips — 400 ms against a store 100 ms away, before a single byte of data.
The measured effect was ten requests down to six. In wall-clock terms on a high-latency link, that is often the largest single improvement available, and it is one line of configuration.
Why scattered access is the worst pattern
Tiles are compressed independently, so the minimum transfer for any pixel is its whole tile — 470 kB on the file measured here.
Reading a thousand scattered pixels therefore costs up to a thousand tiles, or 470 MB, to obtain a thousand values. Reading a block containing them costs a few tiles.
The general principle: make the access pattern match the storage layout. Where points are genuinely spread across a large raster, cluster them and read one block per cluster.
Edge cases or notes
- Check
ds.overviews(1)first. An empty list is the most common cause. block_shapes[0][0] == 1means striped, and no client setting helps.GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIRremoved 4 of 10 requests here.ds.read(1)[::16, ::16]reads everything. Useout_shape.- A 1-pixel read costs a whole tile.
- Open once and reuse the handle; a fresh open costs two requests and clears the cache.
- GDAL caches across
Envblocks, so benchmark configurations in separate processes. - Requester-pays buckets need
AWS_REQUEST_PAYER=requester, or every read fails.
Internal links
- How to read a COG from a URL without downloading the whole file — the settings in full
- Cloud-optimised GeoTIFF explained — why tiles and overviews matter
- How to write a cloud-optimised GeoTIFF in Python — fixing the file
- GDAL cannot open an S3 or HTTPS path — when the open fails outright
- How to read spatial data from S3 and other object storage — credentials and requester-pays
- How to choose chunk and tile sizes that actually help — matching access to layout
- How to profile a slow Python GIS script and find the real bottleneck — general profiling
- Cloud-native geospatial explained — the wider architecture
FAQ
Why does reading a small window download the whole file?
Either the file is striped, so a window reads full-width rows, or you asked for a downsampled read from a file with no overviews.
How do I check whether a remote file is really a COG?
Open it and print ds.block_shapes[0] and ds.overviews(1). A block height of 1 means striped; an empty overview list means no pyramid.
What is the single most useful GDAL setting?
GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR. It removed four of ten requests in a benchmark, all of them probes for files that do not exist.
Why is my thumbnail read so slow?
Almost certainly ds.read(1)[::16, ::16], which reads everything before subsampling. Use out_shape with a resampling method.
Can I fix a file that has no overviews without rewriting it?
Only by adding an external .ovr sidecar next to it, hosted on the same server and permitted by CPL_VSIL_CURL_ALLOWED_EXTENSIONS.
Why is sampling points so slow?
Each point can cost a whole tile — 470 kB here. Read one block covering the points and index into it instead.
Does increasing the cache help?
Only if reads are clustered or repetitive. In a four-scattered-window benchmark, enabling VSI_CACHE changed nothing.