How to Render Raster Tiles from a COG in Python
On this page
Problem statement
Serving a raster on a web map means producing 256-pixel PNG tiles on the XYZ grid. Doing that from a Cloud-Optimised GeoTIFF is efficient only if each low-zoom tile is read from an overview, and the obvious code does not do that. Measured on the Copernicus 30 m DEM over the Alps, one zoom-8 tile cost 30.69 MB reprojected straight from the band — the whole file — and 0.606 MB once the renderer opens the matching overview level (step 3).
Four things have to be right, and each has a characteristic failure:
- reproject into Web Mercator, or the tiles are in the wrong place
- read from the right overview, or every tile costs full resolution
- stretch the values to 0–255, or the tile is black
- handle nodata, or the edges are opaque black rectangles
Quick answer
import mercantile
import numpy as np
import rasterio
from rasterio.warp import reproject, Resampling
from PIL import Image
def overview_level(src_path, tile, size=256):
"""The coarsest overview still finer than the tile's pixels; -1 means full resolution."""
west, _, east, _ = mercantile.bounds(tile)
wanted = (east - west) / size # degrees per output pixel (source is EPSG:4326)
with rasterio.open(src_path) as src:
levels = [i for i, f in enumerate(src.overviews(1)) if src.res[0] * f <= wanted]
return levels[-1] if levels else -1
def render_tile(src_path, tile, size=256, vmin=None, vmax=None, use_overviews=True):
bounds = mercantile.xy_bounds(tile) # Web Mercator metres
dst_transform = rasterio.transform.from_bounds(*bounds, size, size)
level = overview_level(src_path, tile, size) if use_overviews else -1
with rasterio.open(src_path, **({"overview_level": level} if level >= 0 else {})) as src:
data = np.zeros((size, size), dtype="float32")
reproject(rasterio.band(src, 1), data,
dst_transform=dst_transform, dst_crs="EPSG:3857",
resampling=Resampling.bilinear,
src_nodata=src.nodata, dst_nodata=np.nan)
alpha = np.where(np.isfinite(data), 255, 0).astype("uint8")
lo = vmin if vmin is not None else np.nanpercentile(data, 2)
hi = vmax if vmax is not None else np.nanpercentile(data, 98)
scaled = np.clip((data - lo) / max(hi - lo, 1e-9), 0, 1)
rgb = (np.nan_to_num(scaled) * 255).astype("uint8")
return Image.fromarray(np.dstack([rgb, rgb, rgb, alpha]))image = render_tile("dem_cog.tif", mercantile.tile(10.5, 46.5, 12))
image.save("tile_12.png")
print(image.size, image.mode, f"{os.path.getsize('tile_12.png') / 1e3:.1f} kB")(256, 256) RGBA 26.6 kB
reproject into the tile's own transform does the resampling and the reprojection in one call. It does not choose an overview: from rasterio.band(src, 1) it reads the band at full resolution. Opening the file at the right overview_level first is what makes a low-zoom tile cheap. overview_level compares degrees because this source is in EPSG:4326; for a projected source, compare the tile's metres per pixel instead.
Tile bounds in metres, reproject into the tile, use overviews
1. Compute the tile's Web Mercator bounds
tile = mercantile.tile(10.5, 46.5, 12)
print(mercantile.bounds(tile)) # degrees
print(mercantile.xy_bounds(tile)) # metres in EPSG:3857LngLatBbox(west=10.458984375, south=46.498392258597626, east=10.546875, north=46.55886030311718)
Bbox(left=1164288.8148398064, bottom=5860579.832681034, right=1174072.754460309, top=5870363.772301536)
mercantile.bounds gives longitude and latitude; xy_bounds gives Web Mercator metres, which is what the tile transform needs.
2. Reproject directly into the tile
Rather than reprojecting the whole raster and then cutting tiles, reproject each tile's window straight into its 256 × 256 output grid. GDAL then reads only the source pixels that contribute — from whichever dataset you opened. Open the full-resolution band and it reads full-resolution pixels even for a tile covering a hundred kilometres; open the overview level that matches the tile and it reads that.
3. Rely on the overviews
The saving is the whole reason to serve from a COG. One tile at each of three zooms, rendered from the files served over HTTP, counting what GDAL fetched:
for zoom in (8, 10, 12):
tile = mercantile.tile(10.5, 46.5, zoom)
for label, path, use in (("COG, band read directly", "dem_cog.tif", False),
("COG, overview level", "dem_cog.tif", True),
("no overviews", "dem_no_overviews.tif", True)):
n, size = fetched(lambda: render_tile(url(path), tile, use_overviews=use))
print(f"z{zoom:<3} {label:25} {n:3} requests {size / 1e6:7.3f} MB")z8 COG, band read directly 4 requests 30.693 MB
z8 COG, overview level 2 requests 0.606 MB
z8 no overviews 4 requests 30.675 MB
z10 COG, band read directly 3 requests 3.704 MB
z10 COG, overview level 3 requests 1.250 MB
z10 no overviews 3 requests 3.704 MB
z12 COG, band read directly 2 requests 0.655 MB
z12 COG, overview level 2 requests 0.655 MB
z12 no overviews 2 requests 0.655 MB
A zoom-8 tile over a 30 m raster at 46° N is asking for roughly 1/14 resolution. Without overviews, the tile reads every full-resolution pixel it covers and averages them down.
4. Stretch consistently across tiles
Computing percentiles per tile makes each tile stretch to its own range, so adjacent tiles have visibly different contrast and the seams show.
Compute the stretch once from the whole raster's statistics, and pass fixed vmin and vmax to every tile.
5. Use alpha for nodata
A tile is a rectangle; the data usually is not. Without an alpha channel the area outside the data is black — an opaque rectangle over the basemap.
alpha = np.where(np.isfinite(data), 255, 0).astype("uint8")Pass dst_nodata=np.nan to reproject so the fill is distinguishable, and build the alpha from it.
A tile renderer and a global stretch
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 — rewritten with GDAL's COG driver as a COG with 512-pixel tiles and overviews, and again as a tiled GeoTIFF without overviews. A small local HTTP server serves both and logs every request, so the bytes behind each tile can be counted.
import atexit
import itertools
import os
import subprocess
import sys
import rasterio
import requests
from rasterio.shutil import copy as rio_copy
SOURCE_URL = ("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_URL, 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 = r"""
import http.server, os
class Counting(http.server.SimpleHTTPRequestHandler):
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(); 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))
with open("requests.log", "a") as log:
log.write(f"{last - first + 1}\n")
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 = itertools.count(1)
def url(path):
"""A URL GDAL has not cached; the server ignores the query."""
return f"{BASE}/{path}?run={next(RUNS)}"
def fetched(action):
open("requests.log", "w").close()
with rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR"):
action()
sizes = [int(line) for line in open("requests.log")]
return len(sizes), sum(sizes)
with rasterio.open("dem_cog.tif") as ds:
print(ds.width, "x", ds.height, ds.dtypes[0], ds.crs, "| overviews", ds.overviews(1), "| nodata", ds.nodata)3600 x 3600 float32 EPSG:4326 | overviews [2, 4, 8] | nodata None
Example 1 — a tile renderer with a colour map
import mercantile
import numpy as np
import rasterio
from matplotlib import colormaps
from rasterio.warp import reproject, Resampling
from PIL import Image
def render_tile(src_path, tile, size=256, vmin=None, vmax=None,
cmap="viridis", resampling=Resampling.bilinear):
"""One RGBA tile, with nodata transparent and a fixed stretch."""
bounds = mercantile.xy_bounds(tile)
transform = rasterio.transform.from_bounds(*bounds, size, size)
level = overview_level(src_path, tile, size)
with rasterio.open(src_path, **({"overview_level": level} if level >= 0 else {})) as src:
data = np.full((size, size), np.nan, dtype="float32")
reproject(rasterio.band(src, 1), data,
dst_transform=transform, dst_crs="EPSG:3857",
resampling=resampling,
src_nodata=src.nodata, dst_nodata=np.nan)
valid = np.isfinite(data)
if not valid.any():
return None # empty tile: do not write it
lo = vmin if vmin is not None else float(np.nanpercentile(data, 2))
hi = vmax if vmax is not None else float(np.nanpercentile(data, 98))
normalised = np.clip((data - lo) / max(hi - lo, 1e-9), 0, 1)
mapper = colormaps[cmap]
rgba = (mapper(np.nan_to_num(normalised)) * 255).astype("uint8")
rgba[..., 3] = np.where(valid, 255, 0)
return Image.fromarray(rgba)edge = render_tile("dem_cog.tif", mercantile.Tile(270, 181, 9)) # straddles the raster's west edge
sea = render_tile("dem_cog.tif", mercantile.Tile(273, 181, 9)) # east of the raster
alpha = np.asarray(edge)[..., 3]
print(f"edge tile: {np.mean(alpha == 0):.0%} transparent; tile beyond the data: {sea}")edge tile: 22% transparent; tile beyond the data: None
Returning None for an empty tile matters. Writing a fully transparent PNG for every tile over the sea multiplies the pyramid size for no benefit, and clients treat a 404 as "nothing here" perfectly well.
Example 2 — computing the stretch once
import math
import numpy as np
import rasterio
def global_stretch(src_path, band=1, percentiles=(2, 98), max_pixels=4_000_000):
"""Percentile stretch from a downsampled read of the whole raster."""
with rasterio.open(src_path) as src:
factor = max(1, math.ceil(math.sqrt(src.width * src.height / max_pixels)))
data = src.read(band,
out_shape=(src.height // factor, src.width // factor),
resampling=rasterio.enums.Resampling.average,
masked=True)
values = data.compressed()
lo, hi = np.percentile(values, percentiles)
print(f" sampled {values.size:,} pixels at 1/{factor}")
print(f" p{percentiles[0]} {lo:.3f}, p{percentiles[1]} {hi:.3f}, "
f"range {values.min():.3f}..{values.max():.3f}")
return float(lo), float(hi)n, size = fetched(lambda: global_stretch(url("dem_cog.tif")))
print(f" {n} requests, {size / 1e6:.2f} MB") sampled 3,240,000 pixels at 1/2
p2 627.083, p98 3155.531, range 225.177..3870.312
2 requests, 8.18 MB
With two adjacent zoom-11 tiles, the difference between a stretch computed per tile and one computed once shows up at the shared edge — the mean jump in grey value between the last column of one tile and the first column of the next:
left, right = mercantile.Tile(1083, 724, 11), mercantile.Tile(1084, 724, 11)
lo, hi = global_stretch("dem_cog.tif")
for label, kwargs in (("per-tile stretch", {}), ("global stretch", {"vmin": lo, "vmax": hi})):
a = np.asarray(render_tile("dem_cog.tif", left, **kwargs))[..., 0].astype(int)
b = np.asarray(render_tile("dem_cog.tif", right, **kwargs))[..., 0].astype(int)
seam = np.abs(a[:, -1] - b[:, 0]).mean()
within = np.abs(a[:, -2] - a[:, -1]).mean()
print(f" {label:17} jump across the seam {seam:5.1f}, between neighbouring columns inside a tile {within:4.1f}") sampled 3,240,000 pixels at 1/2
p2 627.083, p98 3155.531, range 225.177..3870.312
per-tile stretch jump across the seam 28.3, between neighbouring columns inside a tile 3.4
global stretch jump across the seam 3.9, between neighbouring columns inside a tile 4.1
Reading a downsampled version of the whole raster costs one overview read on a COG and gives statistics representative of the whole scene — which is exactly what a consistent stretch needs. The factor has to round up: the first version of this function used int(), which for this raster's 13 million pixels gave a factor of 1 and read the whole file.
Example 3 — a minimal dynamic tile server
import io
import mercantile
from fastapi import FastAPI, HTTPException, Response
from contextlib import asynccontextmanager
SOURCE = "dem_cog.tif" # or any https:// URL to a COG
STRETCH = None
@asynccontextmanager
async def lifespan(app):
global STRETCH
STRETCH = global_stretch(SOURCE)
print(f" stretch fixed at {STRETCH}")
yield
app = FastAPI(lifespan=lifespan)
@app.get("/tiles/{z}/{x}/{y}.png")
def tile(z: int, x: int, y: int):
if not 0 <= z <= 22:
raise HTTPException(400, "zoom out of range")
image = render_tile(SOURCE, mercantile.Tile(x, y, z),
vmin=STRETCH[0], vmax=STRETCH[1])
if image is None:
raise HTTPException(404, "no data")
buffer = io.BytesIO()
image.save(buffer, format="PNG", optimize=True)
return Response(buffer.getvalue(), media_type="image/png",
headers={"Cache-Control": "public, max-age=86400"})from fastapi.testclient import TestClient
with TestClient(app) as client:
for z, x, y in ((12, 2167, 1448), (9, 273, 181), (23, 0, 0)):
response = client.get(f"/tiles/{z}/{x}/{y}.png")
print(f" {z}/{x}/{y}: {response.status_code} {response.headers['content-type']} {len(response.content):,} bytes") sampled 3,240,000 pixels at 1/2
p2 627.083, p98 3155.531, range 225.177..3870.312
stretch fixed at (627.0834741210938, 3155.53052734375)
12/2167/1448: 200 image/png 42,241 bytes
9/273/181: 404 application/json 20 bytes
23/0/0: 400 application/json 30 bytes
Rendering on demand suits data that changes, or a pyramid too large to pre-generate. The trade-off against a static pyramid is a process to run — and the COG's overviews are what make each request cheap enough for it to work.
Why per-tile reprojection is right
Why reprojecting per tile is right
The alternative — reproject the whole raster to Web Mercator once, then cut tiles — seems more efficient and is worse in two ways.
It materialises a full reprojected copy, which for a large raster is expensive in time and storage. And it resamples twice: once into Web Mercator at some chosen resolution, once again into each tile.
Reprojecting straight into the tile transform resamples once, at exactly the output resolution. Open the overview level that matches the tile and the read is small; reproject from the full-resolution band and it is not — at zoom 8 that was the whole 30 MB file.
Why a per-tile stretch produces seams
A percentile stretch computed from one tile's pixels reflects that tile's range. A tile over a mountain and a tile over a valley have different ranges, so identical values render as different greys.
The result is a visible discontinuity at every tile boundary — the classic symptom of a raster layer that looks like a checkerboard.
Computing the stretch once, from the whole raster, makes every tile consistent. It also means the stretch is a documented parameter rather than an emergent property of the tiling.
Why alpha matters more than it seems
Tiles are rectangles on a square grid, and data extents are not. Every edge tile is partly outside the data.
Without alpha, those pixels take whatever the fill value renders to — usually black, occasionally white — and the layer becomes an opaque rectangle with a ragged data area inside it, covering the basemap.
dst_nodata=np.nan plus an alpha channel from np.isfinite handles it in two lines, and also handles interior nodata such as cloud gaps.
Why resampling method matters per data type
Bilinear for continuous data — elevation, reflectance, temperature — because averaging is meaningful.
Nearest for categorical data — land cover, classifications — because averaging class codes produces codes that do not exist. A land cover tile rendered with bilinear resampling shows boundary pixels in classes that were never observed.
For downsampling by large factors, average uses every contributing source pixel where bilinear samples only four.
Bounds, transforms and stretch statistics
mercantile.xy_boundsfor Web Mercator metres,boundsfor longitude and latitude.- Reproject into the tile transform, not into an intermediate.
- Compute the stretch once, globally, or the tiles show seams.
- Use an alpha channel for nodata, or edge tiles are black rectangles.
- Return 404 for empty tiles rather than writing transparent PNGs.
- Nearest resampling for categorical rasters.
- A COG without overviews makes every low-zoom tile cost full resolution — see the zoom-8 measurement in step 3.
- Set
Cache-Control; tiles are conventionally immutable.
Internal links
- Cloud-optimised GeoTIFF explained — why overviews matter here
- How to write a cloud-optimised GeoTIFF in Python — preparing the source
- Web map tiles explained: XYZ, zoom levels and the tile pyramid — the grid
- How to export a static tile pyramid for offline use — pre-generating instead
- Raster resampling explained — choosing the method
- Rasterio output GeoTIFF is black or empty — the stretch problem in general
- How to reproject a raster in Python with rasterio — the reprojection
- My tiles are offset or in the wrong place — grid mistakes
FAQ
How do I serve a COG as web map tiles?
Open the overview level that matches the tile's resolution, reproject the tile's Web Mercator bounds directly into a 256 × 256 output grid with rasterio.warp.reproject, stretch to bytes with a fixed range, and add an alpha channel for nodata.
Why do my tiles have visible seams?
A per-tile percentile stretch. Compute the stretch once from the whole raster and pass the same vmin and vmax to every tile.
Why are my tile edges black?
No alpha channel. Areas outside the data render as the fill value; set dst_nodata=np.nan and build alpha from np.isfinite.
Does the COG need overviews?
Yes — and the renderer has to use them. A zoom-8 tile cost the whole 30 MB file both without overviews and when reprojected straight from the band of a COG that had them; opening the right overview level brought it to 0.6 MB.
Should I pre-generate tiles or render on demand?
Pre-generate for static data with a bounded extent; render on demand for changing data or very large pyramids. Overviews make on-demand rendering cheap enough to work.
Which resampling method should I use?
Bilinear for continuous data, nearest for classifications, average for large downsampling factors.
What should I return for a tile with no data?
A 404. Clients handle it correctly, and writing transparent PNGs for empty areas multiplies the pyramid for nothing.