How to Read a COG from a URL Without Downloading the Whole File

Problem statement

A Cloud-Optimised GeoTIFF can be read in place: rasterio.open("https://...") and a windowed read fetches only the tiles it needs. That works out of the box, and out of the box it is slower than it should be.

Measured against a byte-counting server, four scattered 512 × 512 window reads from a COG of the Copernicus 30 m DEM over the Alps (step 1):

default settings                 7 requests (1 of them 404s)   2.572 MB
GDAL_DISABLE_READDIR_ON_OPEN     6 requests (0 of them 404s)   2.572 MB
+ VSI_CACHE                      6 requests (0 of them 404s)   2.589 MB

One of the seven requests was GDAL probing for a sidecar file that does not exist. On a local network that is invisible; against an object store 100 ms away it is 100 ms added to every file you open. VSI_CACHE changed nothing here, for a reason the explanation below gives.

Quick answer

import rasterio
from rasterio.windows import Window

SETTINGS = {
    "GDAL_DISABLE_READDIR_ON_OPEN": "EMPTY_DIR",
    "CPL_VSIL_CURL_ALLOWED_EXTENSIONS": ".tif,.TIF,.tiff",
    "VSI_CACHE": "TRUE",
    "VSI_CACHE_SIZE": "50000000",
    "GDAL_HTTP_MERGE_CONSECUTIVE_RANGES": "YES",
}

with rasterio.Env(**SETTINGS):
    with rasterio.open(f"{BASE}/dem_cog.tif") as ds:          # any https:// URL
        window = Window(1800, 1800, 512, 512)
        data = ds.read(1, window=window)
print(data.shape, f"elevations {data.min():.0f}–{data.max():.0f} m")

(512, 512) elevations 1733–3801 m

Set the environment once, at the top of the process, rather than per call. Each rasterio.Env block re-applies the configuration and clears some caches.

Stop the probing, read windows, align them to tiles

1. Turn off directory listing and sidecar probing

GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR stops GDAL listing the "directory" containing the file, which on an object store is an expensive prefix query.

CPL_VSIL_CURL_ALLOWED_EXTENSIONS restricts which sidecar files GDAL will look for — .aux.xml, .tfw, .ovr, .msk. Each missing one is a round trip.

Measured, with the four windows the problem statement describes and the cache cleared between runs:

WINDOWS = [Window(0, 0, 512, 512), Window(3072, 0, 512, 512),
           Window(0, 3072, 512, 512), Window(3072, 3072, 512, 512)]
BASELINE = {}


def four_windows(settings, path="dem_cog.tif"):
    def run():
        with rasterio.Env(**settings), rasterio.open(url(path)) as ds:
            for window in WINDOWS:
                ds.read(1, window=window)
    return fetched(run)


for label, settings in (("default settings", BASELINE),
                        ("GDAL_DISABLE_READDIR_ON_OPEN", {"GDAL_DISABLE_READDIR_ON_OPEN": "EMPTY_DIR"}),
                        ("+ VSI_CACHE", {"GDAL_DISABLE_READDIR_ON_OPEN": "EMPTY_DIR", "VSI_CACHE": "TRUE"})):
    n, size, missing = four_windows(settings)
    print(f"{label:30} {n:3} requests ({missing} of them 404s)  {size / 1e6:6.3f} MB")

default settings                 7 requests (1 of them 404s)   2.572 MB
GDAL_DISABLE_READDIR_ON_OPEN     6 requests (0 of them 404s)   2.572 MB
+ VSI_CACHE                      6 requests (0 of them 404s)   2.589 MB

2. Read windows, not whole bands

aoi_bounds = (10.40, 46.40, 10.45, 46.45)          # about 4 km by 5.5 km, in degrees
with rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR"), rasterio.open(url("dem_cog.tif")) as ds:
    window = rasterio.windows.from_bounds(*aoi_bounds, transform=ds.transform)
    n, size, _ = fetched(lambda: ds.read(1, window=window.round_offsets().round_lengths()))
print(f"window {window.round_lengths().width} x {window.round_lengths().height} px: "
      f"{n} requests, {size / 1e6:.2f} MB of a {os.path.getsize('dem_cog.tif') / 1e6:.1f} MB file")

window 180 x 180 px: 2 requests, 2.46 MB of a 41.6 MB file

The saving is unbounded in the file size: the window's cost depends on the tiles it touches, not on how many tiles the file holds.

3. Align windows to the tile grid

Tiles are compressed independently, so the smallest readable unit is a whole tile. Measured on the COG with 512-pixel tiles, each read from a fresh open at a tile corner:

for label, window in (("1 px", Window(1024, 1024, 1, 1)), ("64 px", Window(1024, 1024, 64, 64)),
                      ("256 px", Window(1024, 1024, 256, 256)), ("512 px, aligned", Window(1024, 1024, 512, 512)),
                      ("512 px, unaligned", Window(1280, 1280, 512, 512))):
    def read_one():
        with rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR"), rasterio.open(url("dem_cog.tif")) as ds:
            ds.read(1, window=window)
    n, size_bytes, _ = fetched(read_one)
    print(f"{label:18} {n:3} requests  {size_bytes / 1e3:9,.1f} kB")

1 px                 3 requests      639.0 kB
64 px                3 requests      639.0 kB
256 px               3 requests      639.0 kB
512 px, aligned      3 requests      639.0 kB
512 px, unaligned    4 requests    2,471.5 kB

A single pixel costs a whole tile. If you need scattered values, read the tile block containing them once rather than issuing a read per point.

4. Use out_shape for downsampled reads

from rasterio.enums import Resampling

for path in ("dem_cog.tif", "dem_no_overviews.tif"):
    def thumb():
        with rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR"), rasterio.open(url(path)) as ds:
            thumbnail = ds.read(1, out_shape=(ds.height // 16, ds.width // 16),
                                resampling=Resampling.average)
    n, size, _ = fetched(thumb)
    print(f"{path:22} {n:3} requests  {size / 1e6:7.3f} MB")

dem_cog.tif              3 requests    0.606 MB
dem_no_overviews.tif    31 requests   30.675 MB

GDAL picks the appropriate overview level automatically — when the file has one.

Never read at full resolution and downsample in NumPy. That is the code path that transfers the whole file.

5. Reuse the dataset handle

Opening a remote dataset costs two requests and the header parse. In a loop over many windows, open once and read many times — the tile cache then works across reads.

Session settings and window reads on a Copernicus DEM tile

The data used on this page

The raster is one tile of the Copernicus 30 m digital elevation model — 1° by 1° of the Alps around the Bernina and Ortler, 3,600 × 3,600 float32 elevations, published by ESA as a COG on AWS. The setup rewrites it with GDAL as a COG with 512 × 512 tiles, and again as a tiled GeoTIFF without overviews, then serves the directory from a small local HTTP server that logs every request — so each read can be counted exactly, sidecar 404s included.

import os
import subprocess
import sys
import time
import atexit

import rasterio
import requests
from rasterio.shutil import copy as rio_copy
from rasterio.windows import Window

SOURCE = ("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N46_00_E010_00_DEM/"
          "Copernicus_DSM_COG_10_N46_00_E010_00_DEM.tif")
if not os.path.exists("dem_source.tif"):
    with requests.get(SOURCE, stream=True, timeout=300) as reply:
        reply.raise_for_status()
        with open("dem_source.tif", "wb") as out:
            for chunk in reply.iter_content(1 << 20):
                out.write(chunk)
if not os.path.exists("dem_cog.tif"):
    rio_copy("dem_source.tif", "dem_cog.tif", driver="COG", BLOCKSIZE=512, COMPRESS="DEFLATE", PREDICTOR="YES")
if not os.path.exists("dem_no_overviews.tif"):
    rio_copy("dem_source.tif", "dem_no_overviews.tif", driver="GTiff", TILED="YES",
             BLOCKXSIZE=512, BLOCKYSIZE=512, COMPRESS="DEFLATE", PREDICTOR=3, COPY_SRC_OVERVIEWS="NO")

SERVER = """
import http.server, os
class Counting(http.server.SimpleHTTPRequestHandler):
    def log(self, status, size):
        with open("requests.log", "a") as log:
            log.write(f"{status} {size} {self.path}\\n")
    def do_HEAD(self):
        path = self.translate_path(self.path.split("?")[0])
        if not os.path.isfile(path):
            self.send_response(404); self.end_headers(); self.log(404, 0); return
        self.send_response(200)
        self.send_header("Content-Length", str(os.path.getsize(path)))
        self.send_header("Accept-Ranges", "bytes")
        self.end_headers(); self.log(200, 0)
    def do_GET(self):
        path = self.translate_path(self.path.split("?")[0])
        if not os.path.isfile(path):
            self.send_response(404); self.end_headers(); self.log(404, 0); return
        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))
        self.log(206 if "Range" in self.headers else 200, last - first + 1)
    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()
"""
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()}"


RUNS = iter(range(1, 10**6))


def url(path):
    """A URL GDAL has not seen before, so nothing is served from its cache; the server ignores the query."""
    return f"{BASE}/{path}?run={next(RUNS)}"


def fetched(action):
    """Requests, bytes and 404s the server saw while action() ran."""
    open("requests.log", "w").close()
    action()
    lines = [line.split() for line in open("requests.log")]
    return len(lines), sum(int(size) for _, size, _ in lines), sum(status == "404" for status, _, _ in lines)


with rasterio.open("dem_cog.tif") as ds:
    print(ds.width, "x", ds.height, ds.dtypes[0], "| tiles", ds.block_shapes[0], "| overviews", ds.overviews(1))
for name in ("dem_source.tif", "dem_cog.tif", "dem_no_overviews.tif"):
    print(f"{name:22} {os.path.getsize(name) / 1e6:6.1f} MB")

3600 x 3600 float32 | tiles (512, 512) | overviews [2, 4, 8]
dem_source.tif           41.8 MB
dem_cog.tif              41.6 MB
dem_no_overviews.tif     30.7 MB

Example 1 — a session-level configuration

import os
import rasterio

REMOTE_SETTINGS = {
    "GDAL_DISABLE_READDIR_ON_OPEN": "EMPTY_DIR",
    "CPL_VSIL_CURL_ALLOWED_EXTENSIONS": ".tif,.TIF,.tiff,.vrt",
    "VSI_CACHE": "TRUE",
    "VSI_CACHE_SIZE": str(64 * 1024 * 1024),
    "GDAL_CACHEMAX": 512,                       # MB of decoded block cache
    "GDAL_HTTP_MERGE_CONSECUTIVE_RANGES": "YES",
    "GDAL_HTTP_MULTIPLEX": "YES",
    "GDAL_HTTP_MAX_RETRY": 3,
    "GDAL_HTTP_RETRY_DELAY": 1,
    "GDAL_BAND_BLOCK_CACHE": "HASHSET",
}


def configure_remote_reads():
    """Apply once per process, before any dataset is opened."""
    for key, value in REMOTE_SETTINGS.items():
        os.environ.setdefault(key, str(value))
    print(f"  configured {len(REMOTE_SETTINGS)} GDAL options for remote reads")

configure_remote_reads()

  configured 10 GDAL options for remote reads

Setting these as process environment variables rather than inside rasterio.Env blocks means they apply to every read including those inside libraries you do not control — rioxarray, stackstac, rasterstats.

The retry settings matter in production: object stores return transient 500s and 503s, and without retries a batch job fails on the first bad response it meets.

Example 2 — reading many windows efficiently

import numpy as np
import rasterio
from rasterio.windows import Window


def read_windows(url, windows, band=1):
    """One open, many reads, so the tile cache is shared."""
    results = []
    with rasterio.open(url) as ds:
        block_h, block_w = ds.block_shapes[band - 1]
        print(f"  tiles are {block_w} x {block_h}")

        # order windows by tile so nearby reads hit the cache
        def tile_key(w):
            return (int(w.row_off // block_h), int(w.col_off // block_w))

        for window in sorted(windows, key=tile_key):
            results.append(ds.read(band, window=window))

    print(f"  read {len(results)} windows from one open handle")
    return results


def align_to_tiles(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)

import random

random.seed(1)
windows = [Window(random.randrange(0, 3400), random.randrange(0, 3400), 200, 200) for _ in range(24)]
n, size, _ = fetched(lambda: read_windows(url("dem_cog.tif"), windows))
print(f"  {n} requests, {size / 1e6:.2f} MB for 24 windows")

  tiles are 512 x 512
  read 24 windows from one open handle
  26 requests, 20.80 MB for 24 windows

Sorting windows by tile is the cheapest optimisation available. Random-order reads evict the cache constantly; tile-ordered reads hit it repeatedly.

Example 3 — sampling points without a request per point

import numpy as np
import rasterio
from rasterio.windows import Window


def sample_points(url, xs, ys, band=1, pad=0):
    """Read one block covering all the points, then index it."""
    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() - pad, 0), min(rows.max() + pad + 1, ds.height)
        col0, col1 = max(cols.min() - pad, 0), min(cols.max() + pad + 1, ds.width)
        window = Window(col0, row0, col1 - col0, row1 - row0)

        span = (row1 - row0) * (col1 - col0)
        print(f"  {len(xs)} points span {col1 - col0} x {row1 - row0} px "
              f"({span / 1e6:.1f} Mpx) — reading as one block")
        if span > 50e6:
            print("  ! the points are spread too widely for one block; "
                  "cluster them and read per cluster")

        block = ds.read(band, window=window)
    return block[rows - row0, cols - col0]

A hundred points inside one 512-pixel tile, then a hundred scattered across the whole raster, each set sampled both ways:

import numpy as np

rng = np.random.default_rng(0)
with rasterio.open("dem_cog.tif") as ds:
    POINT_SETS = {"one tile": rasterio.transform.xy(ds.transform, rng.integers(1030, 1530, 100), rng.integers(1030, 1530, 100)),
                  "whole raster": rasterio.transform.xy(ds.transform, rng.integers(0, 3600, 100), rng.integers(0, 3600, 100))}


def one_by_one(xs, ys):
    with rasterio.open(url("dem_cog.tif")) as ds:
        list(ds.sample(zip(xs, ys)))


for label, (xs, ys) in POINT_SETS.items():
    n, size, _ = fetched(lambda: sample_points(url("dem_cog.tif"), xs, ys))
    print(f"  {label:12} one block: {n:3} requests, {size / 1e6:6.2f} MB")
    n, size, _ = fetched(lambda: one_by_one(xs, ys))
    print(f"  {label:12} ds.sample: {n:3} requests, {size / 1e6:6.2f} MB")

  100 points span 492 x 498 px (0.2 Mpx) — reading as one block
  one tile     one block:   3 requests,   0.66 MB
  one tile     ds.sample:   3 requests,   0.66 MB
  100 points span 3565 x 3439 px (12.3 Mpx) — reading as one block
  whole raster one block:   5 requests,  30.71 MB
  whole raster ds.sample:  50 requests,  29.44 MB

ds.sample() is convenient and issues a read per point. Measured, for a hundred points inside one tile the two were identical — GDAL's block cache absorbed every repeat. Scattered across the raster, both pulled about 30 MB, because a hundred points touched most of the tiles, but ds.sample made 50 requests to the block read's 5 — on a link with 100 ms of latency, four and a half seconds more.

Where points are spread across the raster, cluster them spatially and read one block per cluster.

Sidecar probes, tiles and caching

Why the sidecar probes cost so much

When GDAL opens a dataset it looks for auxiliary files: .aux.xml for statistics, .ovr for external overviews, .msk for masks, world files for georeferencing.

On a local filesystem these are cheap stat calls. Over HTTP each is a request, and a 404 costs a full round trip.

GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR tells GDAL to assume the directory is empty, so it skips them. CPL_VSIL_CURL_ALLOWED_EXTENSIONS narrows which it will consider even when listing is allowed.

Measured with GDAL 3.12: one request — a 404 — removed from a seven-request read, before any data was transferred.

Why the tile is the unit

COG tiles are compressed independently. A client cannot decompress half a tile, so the minimum transfer for any pixel is that pixel's whole compressed tile.

That is why the benchmark shows 1 pixel, 64 × 64, 256 × 256 and an aligned 512 × 512 all costing 639.0 kB — all four touch exactly one tile plus the header — while the same 512 × 512 window shifted off the grid touched four tiles and cost 2,471.5 kB.

The practical rule: make your access pattern tile-shaped. Sampling a raster at scattered points is close to the worst case; reading blocks is close to the best.

Why caching helps less than you expect

VSI_CACHE caches raw byte ranges; GDAL_CACHEMAX caches decoded blocks. Both help when the same tiles are read repeatedly.

In the benchmark, adding VSI_CACHE changed nothing — six requests either way — because each of the four windows touched different tiles, so there was nothing to reuse.

Caching pays when access is repetitive or clustered: a tile server rendering adjacent tiles, a time series reading the same window from many files, a zoom-out followed by a zoom-in. Ordering reads to create that repetition is worth more than increasing the cache size.

Why to open once

Opening a remote dataset costs two requests and a header parse. Inside a loop that is two requests per iteration for nothing.

Worse, a fresh open resets the block cache, so windows that would have hit it miss instead. Keeping the handle open across a loop is often a larger saving than any environment setting.

Settings per process, sorted windows and retries

• Set the environment once per process, not per read.
• GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR is the single highest-value setting.
• Use out_shape for downsampled reads, never full-resolution plus NumPy.
• A 1-pixel read costs a whole tile — 639 kB here.
• Sort windows by tile so the cache is used.
• Open once, read many. A fresh open costs two requests and clears the cache.
• Set retries. Object stores return transient errors, and one 503 should not fail a batch.
• Requester-pays buckets need credentials and AWS_REQUEST_PAYER=requester.

Internal links

• Cloud-optimised GeoTIFF explained — why the tile is the unit
• Reading a COG from a URL is slow or downloads everything — diagnosing a bad read
• Cloud-native geospatial explained — the wider pattern
• How to read spatial data from S3 and other object storage — credentials and requester-pays
• GDAL cannot open an S3 or HTTPS path — when the open itself fails
• How to load Sentinel-2 bands into Python as an analysis-ready array — windowed reads in a real pipeline
• How to choose chunk and tile sizes that actually help — matching access to layout
• How to extract raster values at point locations with rasterio — the point-sampling case

FAQ

How do I read a COG from a URL in Python?

rasterio.open("https://...") and read a window. Set GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR first — with GDAL 3.12 it removed a 404 sidecar probe, one of seven requests, in the benchmark here.

How much data does a windowed read transfer?

Only the tiles the window touches. A 180 × 180 window cost 2.46 MB from a 41.6 MB file, and would cost the same from a much larger one.

Why does reading a single pixel take so long?

Because tiles are compressed independently, so the smallest transfer is one whole tile — 639 kB on the 512-tile COG here.

How do I read a thumbnail cheaply?

ds.read(1, out_shape=(h // 16, w // 16), resampling=Resampling.average). GDAL selects the right overview: 3 requests and 0.606 MB, against 31 requests and 30.7 MB from the same data without overviews.

Should I set VSI_CACHE?

Yes, but expect little from it unless your reads are clustered or repetitive. In a four-scattered-window benchmark it changed nothing.

Why open the file once?

A fresh open costs two requests and clears the block cache. In a loop over windows, opening once is often the largest single saving.

What about authentication?

For S3, set credentials in the environment and AWS_REQUEST_PAYER=requester for requester-pays buckets. Public HTTPS needs nothing.