Generating Web Mercator Tile Pyramids

Reproject at exactly a zoom level’s cell size, with the grid snapped to tile edges:

import morecantile

TMS = morecantile.tms.get("WebMercatorQuad")
zoom = 14
cell = TMS.matrix(zoom).cellSize                 # ~9.55 m at the equator

# resolution=(cell, cell) plus bounds snapped to multiples of the tile size
# makes every pixel edge coincide with a tile edge at this zoom

Get the resolution and the origin right and every zoom level below it lands exactly on a stored overview. This page belongs to preparing rasters for the web in Visualization, Tiling & Web Delivery.

Zoom Levels Are Fixed Resolutions

Which zoom matches which sensor Web Mercator cell size halves with each zoom level: about 153 metres at zoom 10, 38 at zoom 12, 9.55 at zoom 14 and 2.4 at zoom 16. Landsat's 30 metres sits nearest zoom 12, Sentinel-2's 10 metres nearest zoom 14, and sub-metre commercial imagery at zoom 18 or beyond. Cell size halves with every zoom level z10 · 153 m z12 · 38 m Landsat, 30 m z14 · 9.6 m Sentinel-2, 10 m z16 · 2.4 m high-resolution commercial Choose the zoom whose cell size is just finer than the source — not coarser, not much finer. Cell sizes are at the equator; they shrink by the cosine of latitude on the ground.

Picking a maximum zoom much finer than the source resolution is the commonest waste in pyramid generation. Each extra level quadruples the base raster’s pixel count while adding no information, because the pixels are merely upsampled copies of the source.

Environment & Setup

Package Version pin Used for
rasterio >=1.3.0 Warping and writing the aligned raster
morecantile >=5.0 Tile matrix definitions and cell sizes
rio-cogeo >=5.0 Overviews and COG validation 
pip install "rasterio>=1.3.0" "morecantile>=5.0" "rio-cogeo>=5.0"

Complete Working Example

import math

import morecantile
import rasterio
from rasterio.enums import Resampling
from rasterio.transform import from_origin
from rasterio.warp import reproject, transform_bounds

TMS = morecantile.tms.get("WebMercatorQuad")


def best_zoom(src_res_m: float) -> int:
    """Smallest zoom whose cell size is at least as fine as the source."""
    for z in range(0, 24):
        if TMS.matrix(z).cellSize <= src_res_m:
            return z
    return 23


def aligned_web_mercator(src_path: str, dst_path: str, *,
                         resampling: Resampling = Resampling.bilinear) -> int:
    with rasterio.open(src_path) as src:
        src_res = abs(src.transform.a)                  # assumes metres
        zoom = best_zoom(src_res)
        cell = TMS.matrix(zoom).cellSize
        tile_m = cell * 256

        left, bottom, right, top = transform_bounds(src.crs, "EPSG:3857",
                                                    *src.bounds, densify_pts=21)
        # Snap outward to whole tiles so pixel and tile edges coincide
        origin_x = TMS.xy_bounds(morecantile.Tile(0, 0, 0)).left
        origin_y = TMS.xy_bounds(morecantile.Tile(0, 0, 0)).top
        left = origin_x + math.floor((left - origin_x) / tile_m) * tile_m
        right = origin_x + math.ceil((right - origin_x) / tile_m) * tile_m
        top = origin_y - math.floor((origin_y - top) / tile_m) * tile_m
        bottom = origin_y - math.ceil((origin_y - bottom) / tile_m) * tile_m

        width = int(round((right - left) / cell))
        height = int(round((top - bottom) / cell))
        transform = from_origin(left, top, cell, cell)

        profile = src.profile | {
            "crs": "EPSG:3857", "transform": transform,
            "width": width, "height": height,
            "tiled": True, "blockxsize": 256, "blockysize": 256,
            "compress": "deflate",
        }
        with rasterio.open(dst_path, "w", **profile) as dst:
            for b in range(1, src.count + 1):
                reproject(rasterio.band(src, b), rasterio.band(dst, b),
                          src_transform=src.transform, src_crs=src.crs,
                          dst_transform=transform, dst_crs="EPSG:3857",
                          src_nodata=src.nodata, dst_nodata=src.nodata,
                          resampling=resampling, num_threads=4)
            factors = [2 ** i for i in range(1, 7)]
            dst.build_overviews(factors, Resampling.average)
    return zoom

A block size of 256 here is deliberate and differs from the 512 used elsewhere on this site. With the grid snapped to 256-pixel tiles at the maximum zoom, each internal block corresponds exactly to one map tile at that zoom, so serving a tile is a single block read with no resampling at all.

Why Alignment Matters

Aligned and unaligned pixel grids When the raster grid is snapped to the tile grid, each tile boundary falls exactly on pixel edges and a tile is served by copying one block. When the raster is offset by a fraction of a pixel, every tile boundary cuts through pixels, forcing resampling on every request and softening the rendered image. snapped to the tile grid offset by 0.4 pixel tile edges fall on pixel edges every tile edge cuts through pixels Thin lines: pixel edges. Thick lines: tile edges at the maximum zoom.

The unaligned case is not wrong, exactly — the tiler resamples and produces a plausible image. But every tile is resampled on every request, the image is subtly softer than the source, and each request reads a little more than a tile’s worth of data. Over a heavily viewed layer those costs add up, and they are the entire reason to make a Web Mercator copy in the first place.

The alternative to doing this by hand is rio cogeo create --web-optimized, which performs the same snapping and resolution choice. It is the shorter path in practice; knowing what it does is what lets you diagnose a layer that looks soft despite the flag being set — usually a source whose resolution sat awkwardly between two zooms.

Choosing Resampling for the Warp

Reprojecting into Web Mercator resamples every pixel once, and the choice of kernel matters in the usual way. Continuous imagery and indices take bilinear or cubic, which produce a smooth result without inventing extremes; class rasters take nearest, and their overviews take mode, for the reasons set out in choosing the right resampling method for Sentinel-2.

There is one consideration specific to Web Mercator: because the maximum zoom’s cell size rarely equals the source resolution exactly, the warp is also a slight resize — 10 m onto 9.55 m, for instance. That is a 5% upsample, invisible with bilinear resampling and noticeable as occasional duplicated rows and columns with nearest. For a class raster that must stay categorical, accept the duplication; for anything continuous, bilinear hides it completely.

Verification

Four properties of a correct pyramid The raster is in EPSG 3857, its resolution matches a zoom level's cell size to within a fraction of a per cent, its origin sits on a whole-tile boundary, and it carries at least six power-of-two overview levels. Assert all four EPSG:3857 the tile CRS zoom resolution within 0.1% snapped origin whole tiles from 0,0,0 overviews six or more levels A file passing all four serves every zoom from a stored level with no resampling. 
import math

import morecantile
import rasterio

TMS = morecantile.tms.get("WebMercatorQuad")

with rasterio.open("web_copy.tif") as src:
    assert src.crs.to_epsg() == 3857
    res = src.transform.a
    zoom = min(range(24), key=lambda z: abs(TMS.matrix(z).cellSize - res))
    cell = TMS.matrix(zoom).cellSize
    assert abs(res / cell - 1) < 1e-3, f"resolution off grid for z{zoom}"
    origin_x = TMS.xy_bounds(morecantile.Tile(0, 0, 0)).left
    offset = (src.transform.c - origin_x) / (cell * 256)
    assert abs(offset - round(offset)) < 1e-6, "origin not on a tile boundary"
    assert len(src.overviews(1)) >= 6

The origin check is the one that catches a hand-rolled pyramid that got the resolution right and the snapping wrong. Such a file looks nearly perfect and serves every tile with a fraction-of-a-pixel resample — the defect is real but only visible as a slight softness, which is exactly why it is worth asserting rather than eyeballing.

Common Errors

The web layer is blurrier than the source

The resolution is between zoom levels or the origin is not snapped. Check both with the assertions above, or regenerate with the web-optimized option.

The file is far larger than expected

The maximum zoom is one or two levels too fine, quadrupling the pixel count each time. Choose the zoom from the source resolution.

Reprojection fails near the poles

Web Mercator is undefined beyond about 85 degrees of latitude. Clip the source to that band first.

Overviews look blocky for a continuous product

Nearest resampling was used for the overviews. Use average for continuous data and reserve mode for class rasters.

Frequently Asked Questions

Q: What maximum zoom suits 10 metre imagery? Zoom 14, whose cell size at the equator is about 9.55 metres, is the natural match for 10 metre data. Zoom 15 would upsample every pixel by two in each direction and store four times the data for no additional detail.

Q: Why snap the bounds to the tile grid? So that pixel edges coincide with tile edges. An unaligned grid means every tile straddles pixel boundaries and must be resampled, which softens the image and makes each tile read slightly more data than it needs.

Q: Does the web_optimized flag do the same thing? Yes — rio-cogeo’s web_optimized option performs the alignment and resolution choice for you. Doing it by hand is worth understanding so you can tell whether the flag did what you expected when a layer looks soft.

Q: Should the web copy keep the source’s nodata value? Yes, carried through the warp as both source and destination nodata, so that the corners outside the rotated footprint are transparent rather than a solid colour. A missing nodata is the usual cause of black wedges around a reprojected scene.