How to Write a Cloud-Optimised GeoTIFF in Python
On this page
Problem statement
Writing a COG looks like four settings and one extra step — tile it, compress it, pick a predictor, build overviews — and files that miss the extra step are the commonest failure in the whole area, because they look correct. A tiled, compressed GeoTIFF without overviews behaves acceptably for full-resolution windows and badly for anything zoomed out: measured over HTTP, a thumbnail of a 3,600 × 3,600 elevation model cost 0.606 MB with overviews and 30.7 MB without — the whole file.
There is a second failure that the usual recipe does not show. Writing a tiled GTiff in rasterio and calling build_overviews gives a file with overviews that GDAL's own COG validator rejects, because the overview data lands after the full-resolution image instead of before it (step 5). GDAL's COG driver writes the layout correctly, and it is one call.
Quick answer
from rasterio.shutil import copy as rio_copy
from osgeo_utils.samples.validate_cloud_optimized_geotiff import validate
def write_cog(src_path, dst_path, resampling="AVERAGE", blocksize=512):
"""GDAL's COG driver: tiles, compression, predictor, overviews and the layout, in one call."""
rio_copy(src_path, dst_path, driver="COG", BLOCKSIZE=blocksize, COMPRESS="DEFLATE",
PREDICTOR="YES", OVERVIEW_RESAMPLING=resampling, BIGTIFF="IF_SAFER")
warnings, errors, _ = validate(dst_path, full_check=True)
if errors:
raise RuntimeError(f"{dst_path} is not a valid COG: {errors}")
return dst_path
write_cog("dem.tif", "dem_cog.tif")
with rasterio.open("dem_cog.tif") as check:
print(f"blocks {check.block_shapes[0]}, overviews {check.overviews(1)}, "
f"{os.path.getsize('dem_cog.tif') / 1e6:.1f} MB, valid COG")blocks (512, 512), overviews [2, 4, 8], 41.6 MB, valid COG
PREDICTOR="YES" lets the driver pick the right predictor for the data type, and the driver chooses overview levels down to about one tile on its own. The validation is not defensive programming for its own sake: it is the only check that sees the layout.
Tiles, overview levels, resampling, compression, validation
1. Tile it
"tiled": True, "blockxsize": 512, "blockysize": 512Both block dimensions must be multiples of 16, and both must be set — setting only blockysize on a tiled write is an error, and omitting tiled makes the block settings meaningless.
512 is a good default. 256 suits small rasters and tile servers; 1024 suits very large rasters read in big blocks.
2. Choose the overview levels from the image size
Take the coarsest level down to roughly one tile, so a full-extent view is a single read:
for width, height in ((4096, 4096), (3600, 3600)):
largest = max(width, height)
levels, factor = [], 2
while largest / factor >= 256:
levels.append(factor)
factor *= 2
print(f"{width} px: levels {levels}, coarsest {largest // levels[-1]} px")4096 px: levels [2, 4, 8, 16], coarsest 256 px
3600 px: levels [2, 4, 8], coarsest 450 px
The COG driver does this itself — the elevation model above got [2, 4, 8]. The loop is for when you build overviews yourself.
3. Choose the overview resampling for the data
| data | resampling |
|---|---|
| elevation, reflectance, temperature | average |
| land cover, classes, masks | mode or nearest |
| already-binned counts | sum if supported, else average |
Averaging a classification produces class codes that do not exist. Measured on the WorldCover crop, whose legend has ten codes:
LEGEND = {10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 100}
for resampling in ("AVERAGE", "MODE"):
write_cog("landcover.tif", f"landcover_{resampling.lower()}.tif", resampling=resampling)
with rasterio.open(f"landcover_{resampling.lower()}.tif") as ds:
overview = ds.read(1, out_shape=(ds.height // 8, ds.width // 8))
codes = set(np.unique(overview).tolist())
print(f"{resampling:8} overview codes: {len(codes)}, of which {len(codes - LEGEND)} are not classes, "
f"e.g. {sorted(codes - LEGEND)[:6]}")AVERAGE overview codes: 91, of which 80 are not classes, e.g. [11, 12, 13, 14, 15, 16]
MODE overview codes: 9, of which 0 are not classes, e.g. []
This is the same rule as for any resampling, and it is applied once at write time rather than at every read — which makes getting it wrong more consequential.
4. Choose compression and predictor together
"compress": "deflate", "predictor": 2 # integers
"compress": "deflate", "predictor": 3 # floatsThe predictor stores differences between adjacent values rather than the values themselves, which compresses far better on smooth data — and not on classes:
for path, dtype in (("dem.tif", "float32"), ("landcover.tif", "uint8")):
with rasterio.open(path) as src:
data, profile = src.read(1), src.profile
raw = data.nbytes / 1e6
for label, extra in (("no compression", {}), ("deflate", {"compress": "deflate"}),
("deflate, predictor 2", {"compress": "deflate", "predictor": 2}),
("deflate, predictor 3", {"compress": "deflate", "predictor": 3})):
if extra.get("predictor") == 3 and dtype != "float32":
continue
out = f"test_{dtype}.tif"
with rasterio.open(out, "w", **(profile | {"driver": "GTiff", "tiled": True, "blockxsize": 512,
"blockysize": 512} | extra)) as dst:
dst.write(data, 1)
print(f"{path:14} {label:22} {os.path.getsize(out) / 1e6:6.2f} MB ({raw / (os.path.getsize(out) / 1e6):5.2f}x)")dem.tif no compression 67.11 MB ( 0.77x)
dem.tif deflate 44.41 MB ( 1.17x)
dem.tif deflate, predictor 2 37.59 MB ( 1.38x)
dem.tif deflate, predictor 3 30.67 MB ( 1.69x)
landcover.tif no compression 16.78 MB ( 1.00x)
landcover.tif deflate 0.97 MB (17.35x)
landcover.tif deflate, predictor 2 1.14 MB (14.73x)
Three things in that table. The predictor matters for smooth continuous data: predictor 3 took the elevation model from 44.41 MB to 30.67 MB, and predictor 2 on the same floats helped less but still helped. On the land-cover classes it did the opposite — predictor 2 made the file 18% larger than plain deflate, because differences between class codes are not small and long runs of one class compress better as they are. And "no compression" is larger than the data: 3,600 is not a multiple of 512, so the edge tiles are padded out to full size.
5. Validate, always
The check most recipes use — tiled, with overviews — passes files that are not COGs. Here is the common recipe, a tiled GTiff with build_overviews, next to the COG driver's output:
from rasterio.enums import Resampling
with rasterio.open("dem.tif") as src:
profile = src.profile | {"driver": "GTiff", "tiled": True, "blockxsize": 512, "blockysize": 512,
"compress": "deflate", "predictor": 3, "BIGTIFF": "IF_SAFER"}
with rasterio.open("dem_gtiff.tif", "w", **profile) as dst:
dst.write(src.read(1), 1)
dst.build_overviews([2, 4, 8], Resampling.average)
for path in ("dem_gtiff.tif", "dem_cog.tif"):
with rasterio.open(path) as ds:
tiled = ds.block_shapes[0][0] not in (1, ds.height)
looks_right = tiled and bool(ds.overviews(1))
warnings, errors, _ = validate(path, full_check=True)
print(f"{path:14} tiled with overviews: {looks_right}; validator errors: {len(errors)}")
for error in errors:
print(" ", error)dem_gtiff.tif tiled with overviews: True; validator errors: 3
The offset of the first block of overview of index 1 should be after the one of the overview of index 2
The offset of the first block of overview of index 0 should be after the one of the overview of index 1
The offset of the first block of the main resolution image should be after the one of the overview of index 2
dem_cog.tif tiled with overviews: True; validator errors: 0
The GTiff file has tiles and overviews and is not a COG: its overview data sits after the full-resolution image, so a reader asking for a thumbnail has to find it at the end of the file. The validator is osgeo_utils.samples.validate_cloud_optimized_geotiff, shipped with GDAL's Python bindings.
A self-configuring writer for a DEM and WorldCover
The data used on this page
Two real rasters of the same part of the Alps: one tile of the Copernicus 30 m digital elevation model (3,600 × 3,600 float32, published as a COG on AWS), and a 4,096 × 4,096 crop of ESA WorldCover 2021 at 10 m (uint8 land-cover classes), read from its public COG. Both are saved as plain striped GeoTIFFs, the way most rasters arrive before anyone thinks about the cloud.
import os
import numpy as np
import rasterio
from rasterio.windows import from_bounds
DEM = ("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")
WORLDCOVER = ("https://esa-worldcover.s3.eu-central-1.amazonaws.com/v200/2021/map/"
"ESA_WorldCover_10m_2021_v200_N45E009_Map.tif")
with rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR"):
for url, path, window_of in ((DEM, "dem.tif", None),
(WORLDCOVER, "landcover.tif", lambda ds: from_bounds(
10.2, 46.2, 10.2 + 4096 * ds.res[0], 46.2 + 4096 * ds.res[1], ds.transform))):
if os.path.exists(path):
continue
with rasterio.open(url) as src:
window = window_of(src) if window_of else None
data = src.read(1, window=window)
profile = src.profile | {"driver": "GTiff", "tiled": False, "compress": None,
"height": data.shape[0], "width": data.shape[1],
"transform": src.window_transform(window) if window else src.transform}
for key in ("blockxsize", "blockysize", "interleave", "predictor"):
profile.pop(key, None)
with rasterio.open(path, "w", **profile) as dst:
dst.write(data, 1)
for path in ("dem.tif", "landcover.tif"):
with rasterio.open(path) as ds:
print(f"{path:14} {ds.width} x {ds.height} {ds.dtypes[0]:8} striped {ds.block_shapes[0]}, "
f"{os.path.getsize(path) / 1e6:.1f} MB")dem.tif 3600 x 3600 float32 striped (1, 3600), 51.9 MB
landcover.tif 4096 x 4096 uint8 striped (2, 4096), 16.8 MB
Example 1 — a writer that chooses its own settings
The writer is a module, so Example 3 can hand it to worker processes.
%%writefile cogtools.py
import os
import rasterio
from osgeo_utils.samples.validate_cloud_optimized_geotiff import validate
from rasterio.shutil import copy as rio_copy
RESAMPLING = {"continuous": "AVERAGE", "categorical": "MODE", "classified": "NEAREST"}
PREDICTOR = {"continuous": "YES", "categorical": "NO", "classified": "NO"} # see step 4
def write_cog(src_path, dst_path, kind="continuous", blocksize=512, compress="DEFLATE"):
"""Write a validated COG with the resampling the data needs."""
with rasterio.open(src_path) as src:
descriptions = src.descriptions
rio_copy(src_path, dst_path, driver="COG", BLOCKSIZE=blocksize, COMPRESS=compress,
PREDICTOR=PREDICTOR[kind], OVERVIEW_RESAMPLING=RESAMPLING[kind], BIGTIFF="IF_SAFER")
warnings, errors, _ = validate(dst_path, full_check=True)
if errors:
raise RuntimeError(f"not a valid COG: {errors}")
with rasterio.open(dst_path) as check:
assert check.descriptions == descriptions
print(f" {os.path.basename(dst_path)}: {os.path.getsize(dst_path) / 1e6:7.2f} MB, "
f"blocks {check.block_shapes[0]}, overviews {check.overviews(1)}")
return dst_pathfrom cogtools import write_cog as write_validated_cog
write_validated_cog("dem.tif", "dem_example.tif", kind="continuous")
write_validated_cog("landcover.tif", "landcover_example.tif", kind="categorical") dem_example.tif: 41.58 MB, blocks (512, 512), overviews [2, 4, 8]
landcover_example.tif: 1.49 MB, blocks (512, 512), overviews [2, 4, 8]
Example 2 — writing a large raster in windows
import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.windows import Window
def write_cog_windowed(dst_path, shape, transform, crs, generator,
dtype="float32", blocksize=512, nodata=np.nan):
"""Write a COG larger than memory, one block at a time.
generator: yields (Window, ndarray) pairs covering the raster.
"""
height, width = shape
profile = dict(driver="GTiff", height=height, width=width, count=1,
dtype=dtype, crs=crs, transform=transform, nodata=nodata,
tiled=True, blockxsize=blocksize, blockysize=blocksize,
compress="deflate",
predictor=3 if str(dtype).startswith("float") else 2,
BIGTIFF="IF_SAFER")
written = 0
staging = dst_path + ".staging.tif"
with rasterio.open(staging, "w", **profile) as dst:
for window, block in generator:
dst.write(block.astype(dtype), 1, window=window)
written += block.size
write_cog(staging, dst_path) # overviews and the COG layout, then validated
os.remove(staging)
print(f" wrote {written / 1e6:.1f} Mpx through a tiled staging file")
return dst_pathA block-by-block write cannot produce the COG layout directly, because the layout puts the overviews first and they can only be computed once all the data exists. So the blocks go to a tiled staging file and the COG driver rewrites it — reading it back once, which is the price of the layout.
with rasterio.open("dem.tif") as src:
shape, transform, crs = (src.height, src.width), src.transform, src.crs
def feet(path="dem.tif", size=512):
"""Elevation in feet, one 512 x 512 block at a time."""
with rasterio.open(path) as src:
for row in range(0, src.height, size):
for col in range(0, src.width, size):
window = Window(col, row, min(size, src.width - col), min(size, src.height - row))
yield window, src.read(1, window=window) * 3.28084
from rasterio.windows import Window
write_cog_windowed("dem_feet.tif", shape, transform, crs, feet()) wrote 13.0 Mpx through a tiled staging file
Example 3 — converting a directory, in parallel
import glob
import os
from concurrent.futures import ProcessPoolExecutor
import cogtools
def convert_directory(pattern, out_dir, kind="continuous", workers=4):
"""Convert every raster in a directory to a validated COG."""
os.makedirs(out_dir, exist_ok=True)
paths = sorted(glob.glob(pattern))
def target(path):
return os.path.join(out_dir,
os.path.splitext(os.path.basename(path))[0] + ".tif")
todo = [(p, target(p)) for p in paths if not os.path.exists(target(p))]
print(f" {len(paths)} inputs, {len(todo)} to convert")
results, failures = [], []
with ProcessPoolExecutor(workers) as pool:
futures = {pool.submit(cogtools.write_cog, src, dst, kind=kind): src
for src, dst in todo}
for future in futures:
try:
results.append(future.result())
except Exception as exc:
failures.append((futures[future], f"{type(exc).__name__}: {exc}"))
for path, error in failures:
print(f" ! {os.path.basename(path)}: {error}")
print(f" {len(results)} converted, {len(failures)} failed")
return results, failuresFour quadrants of the elevation model as the "directory", converted, then converted again:
os.makedirs("inputs", exist_ok=True)
with rasterio.open("dem.tif") as src:
for i, (row, col) in enumerate(((0, 0), (0, 1800), (1800, 0), (1800, 1800))):
window = Window(col, row, 1800, 1800)
with rasterio.open(f"inputs/quadrant_{i}.tif", "w", **(src.profile | {
"width": 1800, "height": 1800, "transform": src.window_transform(window)})) as dst:
dst.write(src.read(1, window=window), 1)
convert_directory("inputs/*.tif", "cogs")
convert_directory("inputs/*.tif", "cogs") 4 inputs, 4 to convert
4 converted, 0 failed
4 inputs, 0 to convert
0 converted, 0 failed
Skipping outputs that already exist makes the conversion resumable, which matters because COG conversion of a large archive is measured in hours and something always interrupts it.
Overviews, predictors and compression ratios
Why overviews are a separate step
Overviews are additional image data, not a header flag. They have to be computed by reading the full-resolution data and downsampling it, which is why they cannot be part of the profile.
That separation is why they get forgotten. The profile settings are visible at the top of the function and the build_overviews call is buried after the write loop — and omitting it produces a valid, tiled, compressed file that passes every casual inspection.
The only reliable defence is validating the output, which is why every example here does.
Why the predictor matters so much
Deflate compresses repeated byte patterns. Raw pixel values in a smooth raster are all different, so there is little to repeat.
The horizontal predictor stores each value as its difference from the previous one along the row. On smooth data those differences are small and cluster near zero, which compresses very well.
predictor=2 assumes integer arithmetic; predictor=3 is the floating-point variant, which reorders bytes so that exponents and mantissas group together. Measured on the elevation model, 2 on floats gave 37.59 MB against 30.67 MB for 3 — worse, though still better than no predictor. On categorical data a predictor turns long runs of one code into runs of zeros broken by large jumps at every boundary, and the land-cover file grew.
Why compression ratios vary so much
Measured in step 4, the elevation model and the land-cover crop compressed by very different factors with the same settings. That is not a bad setting; it is the data.
Satellite reflectance at 16 bits carries sensor noise in the low bits, and noise is incompressible by definition; smooth fields and classifications with long runs of one value compress well.
If a compression ratio disappoints, check whether the data is genuinely noisy before changing the settings. Quantising reflectance to fewer bits is a lossy choice with real consequences, and worth making deliberately if at all.
Why BIGTIFF=IF_SAFER
Classic TIFF uses 32-bit offsets and cannot exceed 4 GB. BigTIFF lifts that, at the cost of support in very old readers.
IF_SAFER writes BigTIFF when the result might exceed the limit and classic TIFF otherwise. Without it, a write that grows past 4 GB fails at the end — after hours of processing.
COG validation, block sizes and categorical overviews
- Validate with GDAL's COG validator, not by checking for overviews — a tiled GTiff with overviews can still fail it.
- Set both
blockxsizeandblockysize, andtiled=True. - Block sizes must be multiples of 16.
predictor=3for continuous floats,2for smooth integers, none for classes. Predictor 2 grew the land-cover file by 18%.modeornearestoverviews for categorical data, neveraverage.- Use the COG driver for the layout; a GTiff with
build_overviewsstores overviews after the image. BIGTIFF=IF_SAFERso a large write does not fail at the end.- Skip existing outputs so a long conversion is resumable.
Internal links
- Cloud-optimised GeoTIFF explained — what each setting buys
- How to read a COG from a URL without downloading the whole file — the reading side
- Cloud-native geospatial explained — the pattern across formats
- How to choose chunk and tile sizes that actually help — sizing the blocks
- Raster resampling explained — choosing the overview method
- How to reduce GIS file size in Python without wrecking the data — compression trade-offs
- How to batch process rasters with rasterio in Python — converting an archive
- Reading a COG from a URL is slow or downloads everything — what a missing overview causes
- How to embed provenance tags in a GeoTIFF — metadata that travels in the file
FAQ
How do I write a COG in Python?
Copy it with GDAL's COG driver — rasterio.shutil.copy(src, dst, driver="COG", ...) — and validate the result with GDAL's COG validator. A tiled GTiff with build_overviews looks right and fails the validator.
Why does my COG behave like an ordinary GeoTIFF?
Missing overviews, or overviews in the wrong place. Without them a downsampled read costs the full resolution — 30.7 MB against 0.606 MB in a measured comparison — and a GTiff written with build_overviews stores them after the image, which the COG validator rejects.
What block size should I use?
512 × 512 for most rasters, 256 for small ones or tile servers. Both dimensions must be multiples of 16 and both must be set.
What overview levels should I build?
Powers of two down to roughly one tile. For a 4096-pixel image that is [2, 4, 8, 16]; the COG driver chose [2, 4, 8] for the 3,600-pixel elevation model on its own.
Which compression should I use?
Deflate, with predictor=3 for continuous floats and 2 for smooth integers — and no predictor for categorical data, where it made a measured land-cover file 18% larger. ZSTD is often better where readers support it.
Why is my compressed file barely smaller?
Because the data is noisy. Compression depends on the data more than the settings; step 4 shows the same settings giving very different ratios for elevation and land cover.
Do I need GDAL's COG driver?
It guarantees the internal layout. Writing with GTiff plus build_overviews usually produces an acceptable file — validate either way.