Preparing Rasters for the Web

An analytical raster and a web raster want different things. The analytical product wants the sensor’s native projection, its full bit depth, and nothing thrown away; the web rendition wants Web Mercator, eight bits, an aggressive stretch and a pyramid that lines up with tile zoom levels. Dynamic tiling hides much of that difference, but not all of it, and knowing which conversions still pay is the substance of this topic — part of Visualization, Tiling & Web Delivery.

The governing rule is that the web copy is a derivative. It is lossy in value and in geometry, it can be regenerated whenever the rendering decisions change, and it must never become the file people analyse.

---

Prerequisites

pip install "rasterio>=1.3.0" "rio-cogeo>=5.0" "numpy>=1.23" "morecantile>=5.0" "pillow>=10.0"

You need to know the source CRS and the display stretch, which come from extracting and parsing raster metadata and from the rendering decisions in rendering rasters with matplotlib respectively.

---

Step-by-Step Workflow

Step 1 — Decide whether a Web Mercator copy is worth it

The multi-zone case is the one people underestimate. A product covering a country that straddles three UTM zones cannot be a single analytical file without a projection choice, and Web Mercator is a reasonable choice for the display copy precisely because it is a single global grid. It is a poor choice for analysis — area is badly distorted away from the equator — so the two copies serve genuinely different purposes.

Step 2 — Reproject onto the tile grid, not merely into the CRS

import rasterio
from rasterio.enums import Resampling
from rasterio.warp import calculate_default_transform, reproject
import morecantile

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


def reproject_to_tile_grid(src_path: str, dst_path: str, *, zoom: int,
                           resampling: Resampling = Resampling.bilinear) -> None:
    """Reproject to Web Mercator at exactly the resolution of a tile zoom level."""
    target_res = TMS.matrix(zoom).cellSize      # metres per pixel at this zoom

    with rasterio.open(src_path) as src:
        transform, width, height = calculate_default_transform(
            src.crs, "EPSG:3857", src.width, src.height, *src.bounds,
            resolution=(target_res, target_res),
        )
        profile = src.profile | {
            "crs": "EPSG:3857", "transform": transform,
            "width": width, "height": height,
            "compress": "deflate", "tiled": True,
            "blockxsize": 512, "blockysize": 512,
        }
        with rasterio.open(dst_path, "w", **profile) as dst:
            for band in range(1, src.count + 1):
                reproject(
                    rasterio.band(src, band), rasterio.band(dst, band),
                    src_crs=src.crs, dst_crs="EPSG:3857",
                    src_transform=src.transform, dst_transform=transform,
                    src_nodata=src.nodata, dst_nodata=src.nodata,
                    resampling=resampling, num_threads=4,
                )

Passing an explicit resolution is what turns “in Web Mercator” into “on the tile grid”. Without it the default transform picks a resolution derived from the source, tile requests land between pyramid levels, and every tile is resampled from a level that does not match — losing most of the benefit the copy was made for. The general reprojection mechanics are in reprojecting a raster from UTM to WGS84 with rasterio; only the resolution argument is specific to this use.

Step 3 — Convert to display dtype with a recorded stretch

import json

import numpy as np
import rasterio


def to_display_8bit(src_path: str, dst_path: str,
                    limits: dict[int, tuple[float, float]]) -> None:
    """Rescale each band onto 1-255, reserving 0 for nodata, and record the stretch."""
    with rasterio.open(src_path) as src:
        profile = src.profile | {"dtype": "uint8", "nodata": 0,
                                 "compress": "deflate", "tiled": True}
        with rasterio.open(dst_path, "w", **profile) as dst:
            for band in range(1, src.count + 1):
                arr = src.read(band, masked=True).astype("float32")
                lo, hi = limits[band]
                scaled = np.clip((arr - lo) / max(hi - lo, 1e-9), 0, 1) * 254 + 1
                out = scaled.filled(0).astype("uint8")
                dst.write(out, band)
                dst.update_tags(band, stretch_min=str(lo), stretch_max=str(hi))
            dst.update_tags(display_stretch=json.dumps(
                {str(k): list(v) for k, v in limits.items()}))

The stretch must be fixed and recorded, not recomputed per scene. Percentiles derived per scene make every image internally well contrasted and mutually incomparable, so a mosaic of them has visible seams at every scene boundary for reasons that have nothing to do with the ground. Choosing one stretch for a whole collection is the same discipline as fixing a colour map, and the histogram tools in histogram matching across scenes help when the scenes genuinely differ.

Step 4 — Build the pyramid and validate

import rasterio
from rasterio.enums import Resampling
from rio_cogeo.cogeo import cog_translate, cog_validate
from rio_cogeo.profiles import cog_profiles


def finalise_web_cog(src_path: str, dst_path: str, *, categorical: bool = False) -> None:
    resampling = "mode" if categorical else "average"
    cog_translate(
        src_path, dst_path, cog_profiles.get("deflate"),
        overview_level=6, overview_resampling=resampling,
        web_optimized=True,                  # aligns the grid to WebMercatorQuad
        in_memory=False, quiet=True,
    )
    valid, errors, _ = cog_validate(dst_path)
    if not valid:
        raise RuntimeError(errors)

web_optimized=True does the alignment work described in step two automatically, which makes it the shorter route when the stretch and dtype conversion have already happened. It is worth knowing what it does rather than treating it as a magic flag, because when the output looks subtly resampled the reason is almost always that the source resolution did not correspond to a zoom level.

---

Static Tiles and When They Still Win

The object count is the part that surprises people. Each additional zoom level quadruples the tile count, so extending a pyramid from zoom 12 to zoom 14 multiplies it by sixteen. Storage cost is usually fine; the per-object overhead, the listing time and the sheer awkwardness of managing millions of small files are what make static pyramids unattractive for anything that changes.

Where they remain the right answer is genuinely fixed content — a basemap, a one-off published map, an offline deployment on a vessel or in the field with no server — and there the build is a one-time cost paid to remove all serving infrastructure.

import morecantile
from rio_tiler.io import Reader

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


def write_static_tiles(cog_path: str, out_dir: str, *, min_zoom: int = 8,
                       max_zoom: int = 13, rescale: tuple[float, float] = (0, 3000)):
    """Pre-render a pyramid to disk in the conventional z/x/y layout."""
    from pathlib import Path

    with Reader(cog_path) as cog:
        for z in range(min_zoom, max_zoom + 1):
            for tile in TMS.tiles(*cog.geographic_bounds, [z]):
                try:
                    img = cog.tile(tile.x, tile.y, tile.z, tilesize=256)
                except Exception:
                    continue                      # tile outside the data footprint
                path = Path(out_dir) / str(z) / str(tile.x)
                path.mkdir(parents=True, exist_ok=True)
                (path / f"{tile.y}.png").write_bytes(
                    img.rescale(in_range=((rescale[0], rescale[1]),))
                       .render(img_format="PNG"))

---

Thumbnails and Quicklooks

Every published product should carry a small rendered image next to it. It costs a second of compute, makes the archive browsable with no tooling, and catches defects that no schema check will.

import numpy as np
import rasterio
from PIL import Image


def write_thumbnail(src_path: str, png_path: str, *, width: int = 512,
                    limits: tuple[float, float] | None = None) -> None:
    with rasterio.open(src_path) as src:
        factor = max(1, src.width // width)
        arr = src.read(1, out_shape=(max(1, src.height // factor),
                                     max(1, src.width // factor)), masked=True)
    lo, hi = limits or (float(np.percentile(arr.compressed(), 2)),
                        float(np.percentile(arr.compressed(), 98)))
    scaled = np.clip((arr.astype("float32") - lo) / max(hi - lo, 1e-9), 0, 1)
    rgba = np.zeros((*scaled.shape, 4), dtype="uint8")
    rgba[..., :3] = (scaled.filled(0) * 255).astype("uint8")[..., None]
    rgba[..., 3] = np.where(np.ma.getmaskarray(scaled), 0, 255)
    Image.fromarray(rgba, mode="RGBA").save(png_path, optimize=True)

Writing the alpha channel from the mask is what makes a thumbnail honest: nodata is transparent rather than black, so a scene that is half missing looks half missing rather than half dark. Registering the thumbnail as an asset on the product’s catalog record — the thumbnail role described in attaching model metadata to STAC items — means every catalog browser shows it without further work.

---

What Web Mercator Costs You

Reprojecting for display is not free in information, and it is worth being explicit about what is given up so nobody analyses the web copy by mistake.

Three specific losses come with the display copy. Area is distorted by the square of the secant of the latitude, so any hectare figure computed from a Web Mercator raster at mid-latitudes is roughly double the truth — the area work in estimating area with stratified random sampling must be done on the analytical file. Values are quantised to 256 levels by the 8-bit conversion, so subtle gradients become visible steps and no meaningful difference can be computed. And every pixel has been resampled once, so the crisp edges of a class raster are softened unless mode resampling was used throughout.

None of those matter for a map somebody looks at, and all of them matter for a number somebody quotes. Keeping the two files clearly named and clearly documented — *_analysis.tif and *_web.tif, with the web copy’s tags saying what was done to it — prevents the failure where an analyst picks up the display rendition because it was the one that loaded quickly.

---

Parameter Reference

---

Verification & Testing

import morecantile
import rasterio

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

with rasterio.open("web_copy.tif") as src:
    assert src.crs.to_epsg() == 3857, "not in Web Mercator"
    res = abs(src.transform.a)
    zooms = {z: TMS.matrix(z).cellSize for z in range(6, 18)}
    nearest = min(zooms, key=lambda z: abs(zooms[z] - res))
    ratio = res / zooms[nearest]
    assert 0.99 < ratio < 1.01, f"resolution is off the tile grid (z{nearest}, {ratio:.3f})"
    assert src.overviews(1), "no overviews"
    assert src.profile["tiled"], "not internally tiled"

The resolution-to-zoom check is the one worth automating, because a file that is nominally correct but sits between zoom levels produces tiles that are always slightly resampled — a soft, subtly blurred layer that people describe as “looking worse than the original” without being able to say why.

For the 8-bit conversion, verify that the stretch tags round-trip and that the nodata footprint survived:

with rasterio.open("web_display.tif") as src:
    assert src.tags().get("display_stretch"), "stretch not recorded"
    assert src.nodata == 0
    arr = src.read(1)
    assert arr.min() == 0 or (arr > 0).all(), "0 must mean nodata, not a valid value"

---

Troubleshooting

The web layer looks blurrier than the source

The file’s resolution does not match a zoom level, so every tile is resampled. Reproject with an explicit resolution from the tile matrix, or use web_optimized=True.

Scene boundaries are visible in a mosaic of web copies

Each scene was stretched on its own percentiles. Fix one stretch for the collection and apply it everywhere.

The conversion runs out of memory

in_memory defaulted to true for a file just under the threshold. Pass in_memory=False explicitly for anything above a few hundred megabytes.

Static tiles are missing around the edges of the data

Tiles outside the footprint raise rather than returning empty, and the exception handler skipped them — which is correct — but the footprint used for the tile list came from the file’s bounds rather than its valid-data mask. Use the data footprint if the scene has large nodata regions.

The thumbnail is black

Nodata dominated the percentile calculation, or the read was not masked. Read with masked=True and derive limits from compressed() values only.

---

Frequently Asked Questions

Q: Do I have to reproject to Web Mercator? No — a dynamic tiler reprojects per request, and for most layers that is cheap enough. Make a Web Mercator copy when a layer is viewed constantly, when the source spans several UTM zones, or when tiles must be pre-rendered.

Q: Should the web copy replace the analytical file? Never. The 8-bit Web Mercator rendition is a derivative for display and is lossy in both value and geometry. Keep the analytical product as the authoritative file and treat the web copy as something that can be regenerated at any time.

Q: Are static tiles ever the right answer now? For a fixed basemap served at very high rates, or for an offline or air-gapped deployment where no server can run, yes. For anything that changes or needs adjustable rendering, dynamic tiling is simpler and cheaper.

Q: What about JPEG compression for the web copy? It is reasonable for a true-colour display rendition and wrong for anything else. JPEG is lossy in a way that alters pixel values, so a classification or an index copy must stay lossless; for imagery intended only to be looked at, the size saving is substantial and the artefacts are invisible at normal zoom.

Q: How do I keep the web copy from going stale when the product is reprocessed? Generate it in the same pipeline run that produces the analytical file, never as a separate manual step, and write the source product’s identifier into the web copy’s tags. A nightly job that regenerates display copies independently will eventually run against a product that has since been superseded, and nothing in the file itself will say so. Deriving both in one run also means a single failure leaves neither published, which is far easier to reason about than a display layer that quietly describes last month’s data.

---

Related

• Generating Web Mercator Tile Pyramids — grid alignment and zoom-level resolutions in detail.
• Converting Rasters to 8-bit for Display — fixed stretches, nodata handling and recording the transformation.
• Producing PNG Thumbnails and Quicklooks — batch quicklooks as a pipeline quality gate.
• Publishing Static Tiles to Object Storage — layout, upload strategy and cache headers.
• Serving Raster Tiles with TiTiler — the dynamic alternative most of this exists to compare against.