Web Mercator Projection in a WGSL Vertex Shader

Projecting on the GPU means the vertex buffer holds longitude and latitude and the shader turns them into map coordinates per vertex, per frame. It costs a log and a tan per vertex, buys the ability to change projection without touching a buffer, and has two failure modes that are severe enough to be worth knowing before writing the two lines: a latitude at the pole produces an infinity that makes an entire draw call disappear, and the projected output is a full-magnitude coordinate that undoes any precision work done upstream. This page covers the formula, the clamp, the inverse, and the decision about where projection belongs. It is one stage of coordinate precision and projection on the GPU.

What happens to a vertex between the buffer and clip space Four stages left to right. The vertex buffer holds longitude and latitude in degrees as a pair of single-precision floats. The shader converts to radians and clamps the latitude to plus or minus 85.051129 degrees, which is what keeps the logarithm finite. The Mercator formula produces a full-magnitude coordinate in metres, up to about twenty million. The tile origin is subtracted immediately, bringing the value back to a few hundred metres before the matrix multiply takes it to clip space. ONE VERTEX · FOUR STAGES lon, lat degrees ±180, ±90 from the buffer Clamp the latitude ±85.051129° keeps log finite Mercator metres up to 20 037 508 full magnitude Subtract the origin back to a few hundred then the matrix The clamp is not a safety net — it is part of the definition that makes the world square.
The third box is the dangerous one and it exists for only one instruction. Anything that reads the projected value before the subtraction is reading a number with two metres of resolution.

Runnable reference implementation

The forward projection is two expressions. R is the WGS 84 semi-major axis, which is what the Web Mercator definition uses regardless of the fact that it then treats the ellipsoid as a sphere.

wgsl
const PI  : f32 = 3.14159265359;
const R   : f32 = 6378137.0;              // WGS 84 semi-major axis, metres
// ±85.051129° in radians: the latitude whose projection is ±R·π, which is
// what makes the Web Mercator world square.
const LAT_MAX : f32 = 1.4844222297;

fn web_mercator(lon_deg : f32, lat_deg : f32) -> vec2<f32> {
  let lon = radians(lon_deg);
  let lat = clamp(radians(lat_deg), -LAT_MAX, LAT_MAX);
  let y   = R * log(tan(PI * 0.25 + lat * 0.5));
  return vec2<f32>(R * lon, y);
}

struct Frame {
  world_to_clip : mat4x4<f32>,
  origin        : vec2<f32>,     // subtracted after projection, see below
  _pad          : vec2<f32>,     // vec2 pairs keep the 16-byte alignment
};
@group(0) @binding(0) var<uniform> frame : Frame;

@vertex
fn vs(@location(0) lonlat : vec2<f32>) -> @builtin(position) vec4<f32> {
  let projected = web_mercator(lonlat.x, lonlat.y);
  // The subtraction has to happen here, immediately, because `projected`
  // is a full-magnitude value and everything downstream wants a small one.
  return frame.world_to_clip * vec4<f32>(projected - frame.origin, 0.0, 1.0);
}

The inverse is worth having in the same file, because picking and hit-testing need it and reimplementing it on the host is how the two drift apart.

wgsl
fn web_mercator_inverse(x : f32, y : f32) -> vec2<f32> {
  let lon = degrees(x / R);
  let lat = degrees(2.0 * atan(exp(y / R)) - PI * 0.5);
  return vec2<f32>(lon, lat);
}

Parameter reference

Value Setting Why
R 6 378 137 m The WGS 84 semi-major axis. Using the mean radius instead produces coordinates about 0.3 per cent wrong — enough to misalign with every standard tile scheme.
LAT_MAX 1.4844222297 rad 85.051129°, the latitude that makes the projected world square. Not an approximation of 90°; a definition.
Input type f32 degrees Degrees are at most ±180, so f32 carries about six decimal places — roughly 10 cm at the equator. Adequate for a basemap, marginal for survey data.
Origin subtraction immediately after The projected value is full-magnitude; it must be reduced before anything else touches it.
Uniform padding vec2 pairs A vec2<f32> after a mat4x4<f32> needs the next 8-byte boundary; pairing two keeps the struct at a clean 16-byte multiple.
What using the wrong Earth radius costs A bar chart in metres of the horizontal error introduced by three common mistakes, evaluated at latitude 45 degrees. Using the mean Earth radius of 6 371 000 metres instead of the semi-major axis introduces about 21 000 metres of error. Applying an ellipsoidal Mercator rather than the spherical form that EPSG 3857 defines introduces about 14 000 metres. Omitting the pole clamp introduces no error at all until the latitude reaches the pole, at which point the draw call is lost entirely. ERROR AT LATITUDE 45° · METRES Mean radius ≈21 km Ellipsoidal form ≈14 km No pole clamp 0 m — until the pole 0 6000 12000 18000 24000 m wrong by kilometres correct, then catastrophic EPSG:3857 is deliberately a spherical projection of ellipsoidal coordinates — matching it means copying that inconsistency.
The first two are quiet and enormous: a basemap that is fourteen kilometres out still looks like a map. The third is silent and then total, which is arguably the easier failure to debug.

When projecting on the GPU is wrong

The technique has a narrow sweet spot, and it is worth being explicit about the cases outside it.

It is wrong whenever the projected value has to be stored. A compute pass that bins features into a spatial grid, computes a bounding box, or writes a compacted list needs the projected coordinate as data, and a data value at full Mercator magnitude is exactly the precision problem the rest of this topic exists to avoid. In those pipelines the projection belongs on the host or the server, where the origin subtraction can follow it in f64.

It is wrong when the input precision is inadequate. Longitude in f32 degrees resolves to about 10 centimetres at the equator, which is fine for a basemap and not for anything cadastral. Projected metres in a relative frame resolve far better, so a survey-grade pipeline should pre-project and encode rather than shipping degrees.

It is right when the projection has to change. A map that animates between projections, a globe that flattens, an analysis view that offers an equal-area alternative — all of these want the buffer to hold geographic coordinates and the shader to decide what to do with them, because the alternative is re-uploading every vertex on every change.

It is also right for anything drawn directly to the screen with no intermediate storage: points, lines and polygons whose only consumer is the rasterizer. That covers most of a basemap’s overlays, which is why the technique is common despite its narrow applicability.

Zoom, scale and the tile grid

Web Mercator’s other job is to define the tile grid, and the arithmetic connecting the two is worth having in one place because it turns up in every part of a map.

At zoom level z the world is divided into 2^z tiles on each axis, and the projected world spans 2 × R × π metres — about 40 075 016. So one tile at zoom z covers 40075016 / 2^z metres: roughly 156 kilometres at zoom 8, 610 metres at zoom 16, and 38 metres at zoom 20. A 512-pixel tile at zoom 16 therefore holds about 1.2 metres per pixel at the equator.

The equator qualifier is doing real work. Mercator’s scale factor grows as 1 / cos(latitude), so the same tile at 60° north covers half the ground distance it does at the equator, and at 75° north a quarter. That is why a “metres per pixel” figure for a web map is always quoted at the equator and always wrong everywhere else, and why a scale bar has to be recomputed from the current latitude rather than from the zoom alone.

For the shader, the practical consequence is that anything sized in ground units — a road width in metres, a buffer distance, a symbol whose size should be physical rather than screen-relative — needs the latitude-dependent scale factor as a uniform, not a constant derived from zoom.

Failure modes

  • An entire draw call disappears. A latitude of exactly ±90° made tan(π/2) infinite, the multiply propagated the infinity into every clip coordinate, and the primitive was discarded. Fix: the clamp. This is the single most common bug in a hand-written Mercator shader.
  • Geometry is offset by a fraction of a percent from the basemap. R was set to a mean Earth radius (6 371 000) rather than the semi-major axis. Fix: 6 378 137.
  • The map is vertically mirrored. The tile scheme’s y axis grows downward and the projection’s grows upward. Fix: negate y in the projection or in the matrix, once, and write down which.
  • Everything looks right until the camera zooms past about z16. The input degrees in f32 have run out of resolution. Fix: pre-project on the host and encode residuals instead.
  • A polygon that crosses the antimeridian is drawn as a band across the whole world. Not a projection bug — the geometry itself needs splitting at ±180° before it reaches the GPU, because the shader has no way to know that a longitude of −179 and one of +179 are neighbours.
Deciding whether the projection belongs in the shader A decision diagram. The question asks whether the projected coordinate is consumed immediately by the matrix that takes it to clip space, or whether it has to be stored, compared or written to a buffer. On the immediate branch the projection belongs in the vertex shader: the buffer holds degrees, the projection can change between frames, and no precision is lost because the value never persists. On the stored branch the projection belongs on the host or the server, followed by an origin subtraction in double precision, because a stored full-magnitude coordinate has already lost its low bits. CONSUMED NOW, OR STORED? Is the projected value used immediately? straight into the matrix Project in the vertex shader buffer holds degrees yes Projection can change per frame no re-upload Project on the host or server then subtract in f64 no Buffer holds small residuals precision preserved Compute passes that bin or compare coordinates are always on the stored branch.
Everything about this decision follows from one fact: a stored Mercator metre value has two metres of resolution near the world edge, so the only safe place for a full-magnitude coordinate is a register that is about to be multiplied.

Backend / Python interop note

pyproj is the reference for verifying a shader implementation, and a short round-trip test is worth more than reading the formula twice.

python
from pyproj import Transformer

to_merc = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)

def check(lon: float, lat: float) -> tuple[float, float]:
    """Reference values to assert the WGSL implementation against."""
    return to_merc.transform(lon, lat)

# A handful of cases worth pinning: origin, the pole clamp, and the world edge.
CASES = [(0.0, 0.0), (180.0, 85.051129), (-180.0, -85.051129), (150.0, -33.87)]

Two notes on EPSG:3857 specifically. It is defined as a spherical Mercator applied to ellipsoidal coordinates — the latitude is used directly rather than being converted to a conformal latitude first — which is mathematically inconsistent and is precisely what makes it match every web tile scheme. Implementing the ellipsoidal Mercator instead produces coordinates that differ by up to about 20 kilometres at high latitudes, which reads as a basemap that drifts as the user moves north.

Keep the reference cases in the test suite rather than in a notebook. A shader implementation and a pyproj call agreeing to a millimetre on four well-chosen points — the origin, both poles at the clamp, and a high-longitude city — is a stronger guarantee than any amount of visual inspection, and it fails loudly when someone edits a constant.

The second is that always_xy=True matters. Without it, pyproj follows the authority’s axis order for EPSG:4326, which is latitude first — and the resulting transposed comparison will make a correct shader look broken.