Skip to content

Cloud Optimized GeoTIFF (COG) cookbook#

A Cloud Optimized GeoTIFF is a regular GeoTIFF whose internal layout is organized so that HTTP clients can efficiently read just the pixels they need via HTTP range requests. pyramids provides:

  • Dataset.to_cog() — write a COG with full control over the GDAL COG driver's option surface (band subset, dtype cast, NoData, named profiles, tags/colourmap).
  • Dataset.to_cog_bytes() — encode a COG to an in-memory bytes buffer for direct object-store upload (no temp file).
  • Dataset.cog_info() — structured inspection (compression, predictor, blocksize, the overview pyramid, tags).
  • Dataset.is_cog / Dataset.validate_cog() — validate a file.
  • Dataset.read_part() / preview() / point() / read_tile() — overview-decimated partial reads, ideal over /vsicurl/.
  • Dataset.read_file(url) — transparently open COGs from S3, GCS, Azure, and HTTPS via GDAL's virtual filesystem.
  • DatasetCollection.to_cog_stack() — export a temporal stack as a COG-per-slice directory.
  • pyramids cog create|validate|info — a command-line interface.

Defaults. to_cog resolves the two pixel-affecting options per source dtype: the predictor (2 for integer rasters, 3 for float) and the overview resampling (mode for categorical sources — integer dtype or a colour table — and average for continuous float). Pass either explicitly to override. The same policy backs the write_cog(...) facade, so a direct ds.to_cog(...) and write_cog(ds, ...) produce identical output.

Behavior change — default to_cog output

The default output of Dataset.to_cog(path) changed when the write policy was unified. Previously to_cog wrote no predictor and built nearest overviews; it now writes a dtype-aware predictor (2 for integer, 3 for float) and dtype-aware overview resampling (mode for categorical, average for continuous). write_cog on an integer raster likewise switched its overviews from average to mode (which also stops a spurious "averaging corrupts categorical rasters" warning). The files are still valid COGs, but they are byte-different from older releases.

To reproduce the previous output exactly:

from pyramids.dataset import cog

ds.to_cog(
    "scene.tif",
    compression=cog.Compression(predictor=1),
    overviews=cog.Overviews(resampling="nearest"),
)

Writing a COG#

Minimal example:

import numpy as np
from pyramids.dataset import Dataset, GeoReference

arr = np.random.rand(1, 512, 512).astype("float32")
ds = Dataset.from_array(
    arr,
    geo_ref=GeoReference(top_left_corner=(0, 0), cell_size=0.001, epsg=4326),
)
ds.to_cog("scene.tif")

Compression#

Grouped options live under the cog namespace; import it once:

from pyramids.dataset import cog

ds.to_cog("scene.tif", compression=cog.Compression(compress="ZSTD", level=18))
ds.to_cog("scene.tif", compression=cog.Compression(compress="LERC"), extra={"MAX_Z_ERROR": 0.001})
ds.to_cog("scene.tif", compression=cog.Compression(compress="DEFLATE"))   # predictor auto-resolves per dtype
ds.to_cog(
    "scene.tif",
    compression=cog.Compression(compress="DEFLATE", predictor=3),   # force the float predictor
)

Supported compress values: DEFLATE (default), LZW, ZSTD, WEBP, JPEG, LERC, LERC_DEFLATE, LERC_ZSTD, NONE.

By default the predictor is chosen from the source dtype (2 for integer, 3 for float); set cog.Compression(predictor=...) to override.

Named profiles#

Instead of remembering compression + level/quality combos, pass a named profile as the compression string (deflate, zstd, lzw, packbits, jpeg, webp, lerc, lerc_deflate, lerc_zstd, raw):

ds.to_cog("scene.tif", compression="zstd")     # COMPRESS=ZSTD, LEVEL=9
ds.to_cog("scene.tif", compression="lerc_deflate")

# A full Compression object is used directly (no profile expansion):
ds.to_cog("scene.tif", compression=cog.Compression(compress="ZSTD", level=22))

jpeg / webp enforce their dtype/band constraints (Byte; 1-3 / 3-4 bands) up-front with a clear error.

Tile size and bigtiff#

ds.to_cog(
    "scene.tif",
    layout=cog.Layout(
        blocksize=256,        # power of 2 in [64, 4096]
        bigtiff="IF_SAFER",   # default
        num_threads="ALL_CPUS",
    ),
)

Overviews#

The COG driver generates overviews itself; the cog.Overviews object exposes the knobs:

ds.to_cog(
    "scene.tif",
    overviews=cog.Overviews(
        resampling="bilinear",   # default: dtype-aware (see below)
        count=5,                 # default: auto
        compress="DEFLATE",
    ),
)

Categorical-safe by default. When you don't pass overviews=cog.Overviews(resampling=...), to_cog picks mode for categorical sources (integer dtype or a colour table) and average for continuous float — so the default never corrupts category labels (land cover, basin IDs, classification masks). pyramids only emits a UserWarning when you explicitly request an averaging method (average, bilinear, cubic, cubicspline, lanczos) on a categorical raster.

Web-optimized COGs#

For serving to map/tile consumers:

ds.to_cog(
    "web.tif",
    tiling=cog.Tiling(scheme="GoogleMapsCompatible", resampling="bilinear"),
)
# Output is in EPSG:3857 aligned to the web-Mercator zoom grid.

Reading a COG from cloud storage#

pyramids rewrites URL-scheme paths to GDAL's /vsi* form automatically:

URL form Rewritten internally to
s3://bucket/key.tif /vsis3/bucket/key.tif
gs://bucket/key.tif /vsigs/bucket/key.tif
az://container/blob.tif /vsiaz/container/blob.tif
abfs://container/blob.tif /vsiadls/container/blob.tif
abfss://container/blob.tif /vsiadls/container/blob.tif
https://foo/x.tif /vsicurl/https://foo/x.tif
file:///C:/path/x.tif C:/path/x.tif

So you just:

from pyramids.dataset import Dataset

ds = Dataset.read_file("s3://sentinel-cogs/tile/path/B04.tif")
ds = Dataset.read_file("https://example.com/scene.tif")
ds = Dataset.read_file("gs://bucket/key.tif")

Presigned URLs (with ?sig=...) are preserved verbatim.

Credentials via environment variables#

GDAL reads the standard cloud env vars by default; pyramids does not invent its own scheme. Before running your Python process:

# AWS
export AWS_ACCESS_KEY_ID=...
export AWS_SECRET_ACCESS_KEY=...
export AWS_REGION=us-east-1

# Google Cloud Storage
export GS_OAUTH2_REFRESH_TOKEN=...

# Azure
export AZURE_STORAGE_ACCOUNT=...
export AZURE_STORAGE_ACCESS_KEY=...

For public buckets:

export AWS_NO_SIGN_REQUEST=YES

Credentials via CloudConfig#

When you need to override per-operation:

from pyramids.base.remote import CloudConfig

with CloudConfig(aws_region="us-east-1", aws_no_sign_request=True):
    ds = Dataset.read_file("s3://public-bucket/scene.tif")

with CloudConfig(
    aws_access_key_id="AK",
    aws_secret_access_key="SEC",
    aws_region="eu-west-1",
):
    ds = Dataset.read_file("s3://private-bucket/scene.tif")

CloudConfig is a context manager that applies the options for the block and restores the previous values on exit. Where it applies them depends on the option, and that matters for large reads.

Where the options land#

GDAL fetches the byte-ranges of a large read on worker threads, and those threads inherit no thread-local config. So CloudConfig writes each option to the narrowest storage that the threads doing the work can actually read. There are three, and knowing which is which explains every reported surprise:

Per bucket — SetPathSpecificOption. Credentials, endpoint, region and addressing: AWS_ACCESS_KEY_ID, AWS_SECRET_ACCESS_KEY, AWS_S3_ENDPOINT, AWS_REGION, AWS_VIRTUAL_HOSTING, and the GS_*, AZURE_*, OSS_* and SWIFT_* equivalents. This row needs a path= from which a bucket can be derived; without one these stay thread-local instead. Readable by every thread, but only for that bucket. Lasts for the block.

Thread-local — config_options. Transport knobs: GDAL_HTTP_TIMEOUT, VSI_CACHE, retry counts, plus a copy of the routing keys. Readable by only the thread that entered. Lasts for the block.

Process-global — SetConfigOption. The routing subset only: AWS_S3_ENDPOINT, AWS_ENDPOINT_URL, AWS_REGION, AWS_DEFAULT_REGION, AWS_VIRTUAL_HOSTING and AWS_HTTPS. Readable by every thread, for every path. Lasts for one read, not the block.

The third row exists for one reason. GDAL discovers a bucket's region through CPLGetConfigOption, which takes no path, so it can never consult the per-bucket table. A custom endpoint left only in the first two rows is invisible to it: it probes real AWS instead, and a bucket whose name also exists there returns 403. Publishing the routing subset process-globally is the only way those threads can see it.

What this means for you#

  • Credentials never go process-global. They go per bucket (with path=) or thread-local (without), so another thread's read is never signed with your identity. Verified: off any other thread, a credential from this block reads back as None.
  • The routing subset is global only while a read runs — inside merge_rasters, stack_bands, read_windows, or any single read_array / read_part / to_file / info on a dataset that captured a config. Your own with CloudConfig(...) block publishes nothing by itself.
  • While such a read runs, any other store you touch inherits its routing. That includes the destination of a to_file, because the write runs inside the source dataset's own read scope. GDAL resolves path options per key and falls through to the process config for any key a bucket's scope does not set, so an explicit path= on the other store does not shield it — only giving that store its own AWS_S3_ENDPOINT does.

Do not mix two stores inside one read

with CloudConfig(extra=cdse_env):
    ds = Dataset.read_file("/vsis3/eodata/scene.jp2")
    ds.to_file("/vsis3/my-own-bucket/out.tif")   # addressed to CDSE's host

The to_file runs inside the source's read scope, so the destination bucket resolves CDSE's endpoint and its path-style addressing. Write to a local file inside the block and upload after it closes, or give the destination its own CloudConfig with an explicit AWS_S3_ENDPOINT (and AWS_VIRTUAL_HOSTING) naming its real host.

The same applies to two reads of different custom-endpoint stores that overlap in time: there is one process-global slot per key, so the later one wins for both until it finishes. GDAL offers no per-path storage that a worker thread can read, so this cannot be fixed inside pyramids.

Reading an option back

A bucket-scoped option is invisible to gdal.GetConfigOption, which only sees thread-local and process-global config. Inside a block with path=, the workers read AWS_ACCESS_KEY_ID while gdal.GetConfigOption answers None — that is the option working, not missing. To inspect one, use gdal.GetPathSpecificOption(path, key) and pass the /vsi form of the path: the options are keyed on /vsis3/bucket/, so a s3://bucket/... spelling returns None even when the option is set.

Pass path= so it knows which bucket:

with CloudConfig(aws_access_key_id="AK", aws_region="eu-west-1",
                 path="s3://private-bucket/scene.tif"):
    ds = Dataset.read_file("s3://private-bucket/scene.tif")

Without path=#

You do not have to supply it. A block with no path= keeps its options thread-local and records them, and the operation doing the reading — which does know the paths — scopes them onto its own buckets for its duration:

from pyramids.dataset.merge import merge_rasters

with CloudConfig(extra={"AWS_S3_ENDPOINT": "minio.example",
                        "AWS_ACCESS_KEY_ID": "AK"}):
    merge_rasters(["/vsis3/bucket/a.tif", "/vsis3/bucket/b.tif"], "out.tif")

This covers the operations that read inside the call — merge_rasters, stack_bands and read_windows — and any dataset opened from an object store inside the block, which carries the credentials into its later reads. A dataset opened from a local path captures nothing, so it never holds a cloud secret it cannot use.

Operations that fan out their own reads without that wrapper, such as crop, align, to_crs and the gdal.Warp-backed ones, still need an explicit path=.

A captured config has no expiry

A dataset that captured credentials keeps re-installing them on every read, and carries them into its pickle recipe and any dask graph it enters. Reopen the dataset when the credentials rotate.

Validating a COG#

from pyramids.dataset.cog import validate

report = validate("scene.tif")
if report:
    print("valid COG; blocksize=", report.details.get("blocksize"))
else:
    for err in report.errors:
        print("ERROR:", err)
    for warn in report.warnings:
        print("WARN :", warn)

# Strict: warnings (e.g., "no overviews") become errors
strict = validate("scene.tif", strict=True)

Or from a Dataset instance:

ds = Dataset.read_file("scene.tif")
if ds.is_cog:
    print("is a COG")

report = ds.validate_cog(strict=True)

Exporting a DatasetCollection as a COG stack#

One COG per time slice — the typical static-STAC pattern:

from pyramids.dataset import DatasetCollection, cog

dc = DatasetCollection.read_multiple_files("path/to/folder")
dc.open_multi_dataset(band=0)

paths = dc.to_cog_stack(
    "out_cogs/",
    pattern="B04_{i:04d}.tif",   # filename template
    name="B04",
    compression=cog.Compression(compress="ZSTD"),
)
# -> [PosixPath('out_cogs/B04_0000.tif'), ..., PosixPath('out_cogs/B04_NNNN.tif')]

Selecting bands, casting dtype, and NoData on write#

to_cog can pre-process the source in one call — band subset/reorder (indexes, 0-based), dtype cast (out_dtype), and an explicit nodata, grouped in a cog.BandSelection — routed through an in-memory gdal.Translate so your open dataset is never mutated:

# RGB-from-multiband, cast to uint8, set NoData=0:
ds.to_cog(
    "rgb.tif",
    bands=cog.BandSelection(indexes=[3, 2, 1], out_dtype="uint8", nodata=0),
)

The dtype-aware predictor is resolved from the post-cast dtype (so a float→int16 cast correctly yields PREDICTOR=2).

Attaching tags, a colour table, and metadata#

Useful when the source is a bare array that carries no band descriptions or palette:

ds.to_cog(
    "landcover.tif",
    tags=cog.Tags(
        band_tags={0: {"name": "landcover", "units": "class"}},
        colormap={0: (0, 0, 0, 255), 1: (34, 139, 34, 255), 2: (70, 130, 180, 255)},
        metadata={"source": "my-pipeline", "version": "1.0"},
    ),
)

These are stamped onto an in-memory copy before the write, so the original dataset is left untouched.

Encoding to bytes (in-memory)#

To upload a COG straight to an object store without a temp file, use to_cog_bytes — it accepts the same keywords as to_cog:

blob = ds.to_cog_bytes(compression=cog.Compression(compress="ZSTD"))
# e.g. boto3:  s3.put_object(Bucket="b", Key="scene.tif", Body=blob)

Inspecting a COG#

cog_info() reads only headers/metadata (no pixels), so it is cheap even for a large remote COG:

info = ds.cog_info()              # or: from pyramids.dataset.cog import cog_info
info.compression                  # 'DEFLATE'
info.predictor                    # '3'
info.blocksize                    # (512, 512)
info.dtype, info.band_count       # ('Float32', 1)
info.crs_epsg, info.bounds        # (4326, (...))
[o.decimation for o in info.overviews]   # [2, 4, 8]

Partial / overview-decimated reads#

The whole point of a COG is reading only what you need. These methods ask GDAL for a smaller output size, which serves the data from the nearest overview — over /vsicurl/ only the relevant byte ranges are fetched:

# Whole-image thumbnail (long edge <= 1024 px):
thumb = ds.preview(max_size=1024)

# A window, decimated to an explicit output size. `bbox_crs` names the CRS of
# the bbox; omit it when the bbox is already in the raster's own coordinates.
part = ds.read_part(
    (12.4, 41.8, 12.6, 42.0), dst_width=256, dst_height=256, bbox_crs=4326
)

# A single coordinate value, reprojected from `point_crs` when it differs:
value = ds.point(12.5, 41.9, point_crs=4326)

# A Web-Mercator XYZ/slippy-map tile (no morecantile dependency):
tile = ds.read_tile(z=12, x=2191, y=1539, tilesize=256)

read_part / read_tile accept a resampling argument (nearest, bilinear, cubic, cubicspline, lanczos, average, mode) and a band (0-based) to read a single band.

Tuning remote reads from Python#

to_cog, validate, and cog_info accept a config dict applied via gdal.config_options for the duration of the call. Remote reads with no explicit config automatically get cloud-friendly defaults (GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR, curl range tuning):

from pyramids.dataset.cog import validate, cog_info

report = validate("s3://bucket/scene.tif")          # remote defaults applied
info = cog_info("/vsicurl/https://host/scene.tif",  # or override explicitly
                config={"GDAL_HTTP_MAX_RETRY": "5"})

Command line#

The pyramids cog command group exposes the common workflow from the shell:

pyramids cog create input.tif scene_cog.tif --profile zstd
pyramids cog validate scene_cog.tif --strict
pyramids cog info scene_cog.tif

create validates the output by default (--no-validate to skip); validate exits non-zero for a non-COG; info prints the structured metadata and the overview pyramid. The entry point is registered as the pyramids console script (available after pip install pyramids-core).

Performance tuning#

Recommended GDAL environment variables for cloud reads:

export GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR
export GDAL_HTTP_MERGE_CONSECUTIVE_RANGES=YES
export VSI_CACHE=TRUE
export VSI_CACHE_SIZE=26214400       # 25 MiB
export GDAL_CACHEMAX=512             # MiB

GDAL_DISABLE_READDIR_ON_OPEN is the biggest win — it stops GDAL from probing for sidecar files (.aux.xml, .ovr) that usually don't exist on public buckets, saving a directory listing and an extra HTTP request per open. These same options are the COG_READ_DEFAULTS preset that validate() and cog_info() apply automatically to remote paths.

Don't set CPL_VSIL_CURL_ALLOWED_EXTENSIONS

Earlier versions of this guide recommended CPL_VSIL_CURL_ALLOWED_EXTENSIONS=.tif,.TIF,.vrt. Don't. It makes GDAL refuse any URL whose path does not end in one of those extensions — extensionless object keys, presigned S3 links carrying a query string, and most STAC / Planetary Computer asset hrefs — so an option meant to speed reads up instead makes them impossible. It was removed from COG_READ_DEFAULTS for the same reason; GDAL_DISABLE_READDIR_ON_OPEN already delivers the latency win on its own.

Troubleshooting#

"Range downloading not supported by this server!" The remote HTTP server doesn't support HTTP Range headers. Python's stdlib http.server is one such case. Real cloud storage (S3, GCS, Azure) always supports Range.

/vsis3/ returns opaque error for a private bucket Check that AWS_ACCESS_KEY_ID, AWS_SECRET_ACCESS_KEY, and AWS_REGION are set, or use CloudConfig(...) explicitly. For public buckets use AWS_NO_SIGN_REQUEST=YES.

GoogleMapsCompatible output has unexpected bounds The tiling scheme reprojects to EPSG:3857 and aligns to a zoom grid. Your output will always be larger than the input extent — that's by design; it lets a tile server serve the file without further reprojection.

overviews=cog.Overviews(resampling="average") warning on my categorical raster This is the documented safety check. Either switch to cog.Overviews(resampling="nearest") / "mode", or suppress the warning if you really want averaged overviews:

import warnings

from pyramids.dataset import cog

with warnings.catch_warnings():
    warnings.simplefilter("ignore", UserWarning)
    ds.to_cog("x.tif", overviews=cog.Overviews(resampling="average"))