How to Write a Cloud-Optimised GeoTIFF in Python
Problem statement
Writing a COG is four settings and one extra step, and files that miss the extra step are the most common failure in the whole area — because they look correct.
A tiled, compressed GeoTIFF without overviews behaves acceptably for full-resolution windows and catastrophically for anything zoomed out. Measured on the same 4096 × 4096 raster served over HTTP:
whole scene at 1/16 resolution
with overviews 1 request 0.115 MB
without overviews 3 requests 24.672 MB
The file without overviews transferred all 24.7 MB of itself to produce a 256 × 256 thumbnail.
Quick answer
import rasterio
from rasterio.enums import Resampling
def write_cog(src_path, dst_path, blocksize=512, levels=(2, 4, 8, 16)):
with rasterio.open(src_path) as src:
profile = src.profile | {
"driver": "GTiff",
"tiled": True,
"blockxsize": blocksize,
"blockysize": blocksize,
"compress": "deflate",
"predictor": 3 if src.dtypes[0].startswith("float") else 2,
"BIGTIFF": "IF_SAFER",
}
with rasterio.open(dst_path, "w", **profile) as dst:
for band in range(1, src.count + 1):
dst.write(src.read(band), band)
dst.build_overviews(list(levels), Resampling.average)
dst.update_tags(ns="rio_overview", resampling="average")
with rasterio.open(dst_path) as check:
assert check.overviews(1), "no overviews were written"
print(f" blocks {check.block_shapes[0]}, overviews {check.overviews(1)}")The assertion is not defensive programming for its own sake. build_overviews can silently do nothing if the levels are larger than the image.
Step-by-step solution
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:
largest = max(width, height)
levels, factor = [], 2
while largest / factor >= 256:
levels.append(factor)
factor *= 2For 4096 pixels that gives [2, 4, 8, 16], whose coarsest level is 256 pixels — one tile.
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. 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. Using predictor=2 on floating-point data makes files larger.
Measured on a real Sentinel-2 band, a 4096 × 4096 uint16 crop compressed to 33.4 MB against 34.4 MB raw — satellite reflectance is noisy and barely compresses. Elevation and classifications do far better with the same settings.
5. Validate, always
with rasterio.open(dst_path) as ds:
tiled = ds.block_shapes[0][0] not in (1, ds.height)
assert tiled and ds.overviews(1)
Code examples
Example 1 — a writer that chooses its own settings
import os
import rasterio
from rasterio.enums import Resampling
RESAMPLING = {"continuous": Resampling.average,
"categorical": Resampling.mode,
"classified": Resampling.nearest}
def write_cog(src_path, dst_path, kind="continuous", blocksize=512,
compress="deflate", min_overview_px=256, nodata=None):
"""Write a validated COG, deriving overview levels from the image size."""
with rasterio.open(src_path) as src:
largest = max(src.width, src.height)
levels, factor = [], 2
while largest / factor >= min_overview_px:
levels.append(factor)
factor *= 2
is_float = src.dtypes[0].startswith("float")
profile = src.profile | {
"driver": "GTiff", "tiled": True,
"blockxsize": blocksize, "blockysize": blocksize,
"compress": compress,
"predictor": 3 if is_float else 2,
"BIGTIFF": "IF_SAFER",
}
if nodata is not None:
profile["nodata"] = nodata
with rasterio.open(dst_path, "w", **profile) as dst:
for band in range(1, src.count + 1):
dst.write(src.read(band), band)
if src.descriptions[band - 1]:
dst.set_band_description(band, src.descriptions[band - 1])
if levels:
dst.build_overviews(levels, RESAMPLING[kind])
dst.update_tags(ns="rio_overview", resampling=kind)
with rasterio.open(dst_path) as check:
tiled = check.block_shapes[0][0] not in (1, check.height)
size = os.path.getsize(dst_path)
print(f" {os.path.basename(dst_path)}: {size / 1e6:7.2f} MB, "
f"blocks {check.block_shapes[0]}, overviews {check.overviews(1)}")
if not tiled:
raise RuntimeError("output is striped, not tiled")
if levels and not check.overviews(1):
raise RuntimeError("overviews were requested but not written")
return dst_pathExample 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
with rasterio.open(dst_path, "w", **profile) as dst:
for window, block in generator:
dst.write(block.astype(dtype), 1, window=window)
written += block.size
levels, factor = [], 2
while max(height, width) / factor >= 256:
levels.append(factor); factor *= 2
dst.build_overviews(levels, Resampling.average)
print(f" wrote {written / 1e6:.1f} Mpx, overview levels {levels}")
return dst_pathOverviews are built after all the data is written, from inside the same open handle. Building them before the data exists produces empty overviews; building them from a reopened file works too but reads everything back.
Example 3 — converting a directory, in parallel
import glob
import os
from concurrent.futures import ProcessPoolExecutor
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(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, failuresSkipping 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.
Explanation
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. Using 2 on floats produces differences that are not small in the byte representation, and the file grows.
Why compression ratios vary so much
The measured Sentinel-2 crop went from 34.4 MB raw to 33.4 MB compressed — 3%. 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. A smooth elevation model at the same size and settings routinely compresses by 60–80%, and a classification raster with twenty classes by more still.
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.
Edge cases or notes
- Validate that overviews exist after writing. This is the common failure.
- Set both
blockxsizeandblockysize, andtiled=True. - Block sizes must be multiples of 16.
predictor=3for floats,2for integers. The wrong one grows the file.modeornearestoverviews for categorical data, neveraverage.- Build overviews after writing all the data, from the same handle.
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
FAQ
How do I write a COG in Python?
Write a tiled, compressed GeoTIFF and call build_overviews before closing the file. Then validate that the overviews exist.
Why does my COG behave like an ordinary GeoTIFF?
Almost always missing overviews. Without them a downsampled read costs the full resolution — 24.7 MB against 0.115 MB in a measured comparison.
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], whose coarsest level is 256 pixels.
Which compression should I use?
Deflate with predictor=2 for integers and 3 for floats. ZSTD is better where readers support it.
Why is my compressed file barely smaller?
Because the data is noisy. Satellite reflectance compressed 3% in a measured case; smooth elevation routinely compresses 60–80% with identical settings.
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.