Fixing Misaligned Tiles and Wrong Web Mercator Extents

The tile set loads, the imagery appears, and it is not quite where it should be: a building offset from its footprint on the basemap, a road running parallel to the one in the orthomosaic. Or everything is in the right place and visibly softer than the source. Both are geometry faults in the warp rather than problems with the tiles, and there are four of them.

Four ways a tile set ends up wrong

Pixel misalignment. Without target-aligned pixels the warped raster’s origin falls wherever the reprojected corner landed, which is not on a tile boundary. Every tile then resamples the source at a fractional offset, and the whole set is softer than it should be. Position is correct; sharpness is lost.

Resampling choice. Nearest-neighbour on visual imagery produces aliasing at every zoom; cubic on a categorical raster invents values. Position is correct; content is wrong.

Sphere versus ellipsoid. Web Mercator is defined on a sphere of the WGS84 semi-major axis while the coordinates it consumes are ellipsoidal. Treating EPSG:3857 as an ordinary ellipsoidal Mercator introduces a latitude-dependent northing error reaching about 20 km at mid-latitudes. Modern GDAL handles this correctly; hand-rolled conversion formulas frequently do not.

The latitude limit. The Web Mercator grid is square, which requires truncating at about ±85.051°. A survey above that cannot be tiled in this scheme at all, and a naive transformation returns a northing that grows without bound rather than failing.

Four faults, and which of position and sharpness each affects A two-by-two classification of tile faults against their symptoms. Pixel misalignment leaves position correct and sharpness degraded. Wrong resampling leaves position correct and content wrong, with aliasing on imagery or invented values on a classification. Treating Web Mercator as ellipsoidal shifts position by up to twenty kilometres at mid-latitudes while sharpness is unaffected. Exceeding the latitude limit produces coordinates that grow without bound, which fails outright. A note states that separating the position faults from the sharpness faults is the first useful question. position correct position wrong pixel misalignment no -tap on the warp every tile softer than the source wrong resampling nearest on imagery → aliasing cubic on classes → invented values ellipsoidal Mercator assumed 3857 is defined on a sphere up to ≈ 20 km northing error beyond ±85.051° the grid does not extend there northing grows without bound "Is it in the wrong place, or is it soft?" separates the left column from the right before anything is opened.

Figure 1 — Two symptoms, four causes. The first question narrows it to two candidates, and each of those has a one-line check.

Minimal reproducible solution

Check the warped raster’s geometry directly. All four faults are visible in its metadata and bounds, before a single tile is generated.

import math

import rasterio
from rasterio.warp import transform_bounds

WEB_MERCATOR_LIMIT_DEG = 85.05112877980659
EARTH_CIRCUM_M = 2 * math.pi * 6378137.0


def check_warped_for_tiling(path: str) -> list[str]:
    """Every geometry fault a tile pyramid can inherit, from the warped raster."""
    problems: list[str] = []
    with rasterio.open(path) as ds:
        if ds.crs.to_epsg() != 3857:
            problems.append(f"CRS is {ds.crs}, not EPSG:3857")
            return problems

        res = abs(ds.transform.a)
        ox, oy = ds.transform.c, ds.transform.f

        # -tap alignment: the origin must be an integer multiple of the
        # resolution measured from the grid origin at the top-left corner.
        half = EARTH_CIRCUM_M / 2
        for name, origin in (("x", ox + half), ("y", half - oy)):
            frac = (origin / res) % 1.0
            if min(frac, 1 - frac) > 1e-6:
                problems.append(
                    f"{name} origin is {frac:.4f} of a pixel off the grid — "
                    "warp without -tap; every tile resamples at an offset")

        west, south, east, north = transform_bounds(ds.crs, "EPSG:4326", *ds.bounds)
        if max(abs(north), abs(south)) > WEB_MERCATOR_LIMIT_DEG:
            problems.append(
                f"extent reaches {max(abs(north), abs(south)):.3f}° — beyond the "
                f"Web Mercator limit of {WEB_MERCATOR_LIMIT_DEG:.3f}°")

        if ds.count < 4 and ds.nodata is None:
            problems.append("no alpha band and no NoData — the footprint will "
                            "render as opaque")

        tags = ds.tags()
        if tags.get("RESAMPLING", "").lower() in {"near", "nearest"} and ds.count >= 3:
            problems.append("nearest-neighbour resampling on multi-band imagery "
                            "— expect aliasing at every zoom")
    return problems

The alignment test is the one that repays writing. It converts “the tiles look soft” — a subjective judgement nobody can act on — into a number: the fraction of a pixel by which the grid is offset. Zero means aligned; anything else means the warp needs -tap and nothing about the tile generation will help.

The sphere-versus-ellipsoid fault does not appear in this list because GDAL handles it correctly. It appears when someone converts coordinates by hand:

import math


def lonlat_to_3857(lon_deg: float, lat_deg: float) -> tuple[float, float]:
    """The correct spherical formula. Note there is no eccentricity term.

    Web Mercator is defined on a sphere of the WGS84 semi-major axis, so any
    formula containing the ellipsoid's eccentricity is computing a different
    projection — one that differs by kilometres at mid-latitudes.
    """
    if abs(lat_deg) > WEB_MERCATOR_LIMIT_DEG:
        raise ValueError(f"latitude {lat_deg} is outside the Web Mercator grid")
    r = 6378137.0
    x = r * math.radians(lon_deg)
    y = r * math.log(math.tan(math.pi / 4 + math.radians(lat_deg) / 2))
    return x, y

Edge-case matrix

Fault Position Sharpness Check
No -tap correct degraded Origin modulo resolution
Nearest on imagery correct aliased Resampling tag, band count
Cubic on a classification correct invented values Band count, dtype
Ellipsoidal formula up to 20 km off unaffected Compare against pyproj
Beyond ±85.051° fails or wraps — Bounds check
Source CRS undeclared arbitrary — CRS check before warping
Alpha absent correct correct Band count and NoData
Warped twice correct doubly degraded Provenance in the manifest

The warped-twice row is a quiet one. A raster reprojected to Web Mercator by an earlier step and then handed to a tiler that warps again is resampled twice, and the softness compounds. Recording the CRS in the run manifest at each step is what makes it visible; the raster itself looks fine.

Verification snippet

Position is best verified against something external — a known coordinate, checked through the tile grid arithmetic rather than through the same library that produced the raster.

import math

import rasterio


def assert_known_point_lands_correctly(tile_dir, warped_path: str,
                                       lon: float, lat: float,
                                       zoom: int, tol_px: float = 1.0) -> None:
    """A surveyed point must fall at the expected pixel of the expected tile."""
    n = 2 ** zoom
    x = (lon + 180.0) / 360.0 * n
    y = (1.0 - math.asinh(math.tan(math.radians(lat))) / math.pi) / 2.0 * n
    tx, ty = int(x), int(y)
    px, py = (x - tx) * 256, (y - ty) * 256

    tile = tile_dir / str(zoom) / str(tx) / f"{ty}.png"
    assert tile.exists(), f"expected tile {zoom}/{tx}/{ty} does not exist"

    # Independently, where does the warped raster put the same coordinate?
    with rasterio.open(warped_path) as ds:
        from pyproj import Transformer
        tr = Transformer.from_crs("EPSG:4326", ds.crs, always_xy=True)
        mx, my = tr.transform(lon, lat)
        row, col = ds.index(mx, my)

    # Convert the raster position into the same tile-pixel frame.
    res = abs(ds.transform.a)
    world_per_tile = 2 * math.pi * 6378137.0 / n
    exp_px = ((mx + math.pi * 6378137.0) % world_per_tile) / world_per_tile * 256
    assert abs(exp_px - px) <= tol_px, (
        f"the raster places the point at pixel {exp_px:.1f} of its tile while "
        f"the tile grid expects {px:.1f} — the warp is not aligned to the grid")

Computing the expected tile from the tile-grid formula and the actual position from the raster’s own transform is what makes this a real check: two independent paths to the same answer, which only agree if the warp respected the grid.

What pixel misalignment does at a tile boundary Two tile boundaries drawn over a source pixel grid. In the aligned case the tile edge coincides exactly with a source pixel edge, so each output pixel takes the value of one source pixel and the tile is as sharp as the source. In the misaligned case the tile edge falls halfway through a source pixel, so every output pixel is a weighted blend of two source pixels and the whole tile is blurred by half a pixel. A note observes that the effect is uniform across the tile rather than confined to its edge, which is why the whole pyramid looks soft rather than showing visible seams. aligned — tile edge on a pixel edge tile boundary one source pixel per output pixel as sharp as the source misaligned — edge mid-pixel tile boundary, half a pixel in every output pixel blends two source pixels uniformly soft, with no visible seam The blur is uniform across the tile, not confined to its edge. Which is why it reads as "the imagery is a bit soft" rather than as an alignment fault.

Figure 2 — Why misalignment is easy to live with and worth fixing. It produces no artefact, no seam and no error — only a uniform loss of the resolution the survey paid for.

When to escalate

  • The alignment check passes and the imagery is still soft. The source was warped more than once, or the top zoom exceeds the source resolution so the client is looking at upsampled tiles. Both are visible in the run manifest if the CRS and the zoom range were recorded.
  • A survey extends beyond ±85°. Web Mercator cannot represent it. Use a polar projection and a tile scheme that supports it; there is no adjustment to the warp that makes this work.
  • The tiles align and the client’s basemap does not. Their basemap may be in a different scheme or a different datum realisation. Compare a surveyed coordinate against both rather than assuming the newer product is the wrong one.

Run the geometry checks on the warped raster rather than on the tiles, and run them before generation. Every fault on this page is present in the warped raster and takes milliseconds to detect there; the same fault found after a four-hour pyramid build costs the build.

← Tiling and Serving Orthomosaics as XYZ Tiles

Northing error from treating Web Mercator as ellipsoidal A curve of northing error against latitude for coordinates converted with an ellipsoidal Mercator formula instead of the spherical definition Web Mercator uses. The error is zero at the equator, reaches about eleven kilometres at thirty degrees, about twenty kilometres at forty-five degrees, and falls back toward the poles. A note observes that the error is a smooth function of latitude, so a survey converted this way is displaced uniformly and looks internally correct while sitting kilometres from its true position. 0° 20° 45° 65° 85° latitude northing error ≈ 20 km at 45° Smooth in latitude, so the survey is displaced uniformly and looks internally perfect while being kilometres out.

Figure 3 — The error from one missing assumption. GDAL never makes it; hand-rolled conversion code makes it regularly, and the result is a survey that is internally consistent and in the wrong country’s worth of offset.