⚙️ How-To Web maps & tiles Cloud-native formats Publishing Web Maps · step 11 of 15 Updated

How to Render Raster Tiles from a COG in Python

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 — measured, a whole-scene thumbnail from a COG costs 0.115 MB in one request against 24.672 MB from the same data without overviews.

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 calculate_default_transform, reproject, Resampling
from PIL import Image


def render_tile(src_path, tile, size=256, vmin=None, vmax=None):
    bounds = mercantile.xy_bounds(tile)          # Web Mercator metres
    dst_transform = rasterio.transform.from_bounds(*bounds, size, size)

    with rasterio.open(src_path) 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]), mode="RGBA")

reproject into the tile's own transform does the resampling, the reprojection and the overview selection in one call.

A COG reprojected into a Web Mercator tile transform, stretched to bytes, with nodata becoming a transparent alpha channel.
One reprojection call per tile. GDAL picks the overview level from the output resolution.

Step-by-step solution

1. Compute the tile's Web Mercator bounds 

bounds = mercantile.xy_bounds(tile)     # metres in EPSG:3857

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, at the resolution the output needs — which means it uses an overview automatically.

3. Rely on the overviews

The saving is the whole reason to serve from a COG:

whole scene at 1/16 resolution
  with overviews     1 request    0.115 MB
  without           3 requests   24.672 MB

A tile at zoom 10 over a 10 m raster is asking for roughly 1/9 resolution. Without overviews, every one of those tiles reads full-resolution pixels 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.

Adjacent tiles stretched independently showing a visible seam, against a global stretch producing continuous tone.
A per-tile stretch is the most common cause of a checkerboard-looking raster layer.

Code examples

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)

    with rasterio.open(src_path) 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, mode="RGBA")

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 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, int(np.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)

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.

Example 3 — a minimal dynamic tile server 

import io
import mercantile
from fastapi import FastAPI, HTTPException, Response

app = FastAPI()
SOURCE = "https://example.com/scene.tif"
STRETCH = None


@app.on_event("startup")
def startup():
    global STRETCH
    STRETCH = global_stretch(SOURCE)
    print(f"  stretch fixed at {STRETCH}")


@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"})

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.

Explanation

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, from the source, at exactly the output resolution. GDAL selects the overview whose resolution best matches, so the read is small.

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.

A tile without an alpha channel showing nodata as an opaque rectangle against one with alpha letting the basemap through.
Two lines of alpha handling, and every edge tile stops covering the basemap.

Edge cases or notes

  • mercantile.xy_bounds for Web Mercator metres, bounds for 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 tile cost full resolution — 215× more bytes in a measured case.
  • Set Cache-Control; tiles are conventionally immutable.

FAQ

How do I serve a COG as web map tiles?

Reproject each 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?

Effectively yes. Without them every zoomed-out tile reads full-resolution pixels — 24.7 MB against 0.115 MB for a measured thumbnail.

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.