Chapter 4: Web Mercator and Tile Coordinates

Introduction

The Mercator projection and Slippy Map tile coordinates are familiar concepts for GIS professionals, but while 2D maps simply arrange square tiles on a plane, 3D rendering requires drawing them as curves on a spherical surface.

In this chapter, you can observe in 3D how the tile grid that looks familiar on a flat map appears on a sphere. Switching zoom levels lets you intuitively understand how the tile subdivision pattern changes on the spherical surface.

What You Will Learn

• The mathematics behind the Mercator projection -- deriving the forward and inverse transformation formulas
• Characteristics and limitations of Web Mercator (EPSG:3857) -- differences from the strict Mercator projection
• The mathematical origin of the valid latitude range (+/-85.05 degrees) -- the constraint required to form square tiles
• How Slippy Map Tiles (XYZ tiles) work -- the z/x/y coordinate conversion algorithm
• Conversion between longitude/latitude and tile coordinates -- deriving X coordinates and deriving Y coordinates via Mercator
• Tile bounding box calculation and tile URL conventions -- practical tile operations

Zoom Level

Tile count: 16 (4 x 4)

What Is a Map Projection?

The Earth is a three-dimensional ellipsoid, but map tiles are two-dimensional square images. The operation of converting a 3D curved surface to a 2D plane is called a map projection. There are many types of projections, but web map tile systems almost exclusively use the Mercator projection.

Mathematics of the Mercator Projection

The Mercator projection is a conformal cylindrical projection devised by Gerardus Mercator in 1569. "Conformal" means that local angles are preserved at any point on the map (i.e., locally small shapes retain their form).

The formula to convert longitude lam (radians) and latitude phi (radians) to Mercator coordinates (x, y):

Here, a is the WGS84 semi-major axis (6,378,137 m). The longitude direction maps linearly, but the latitude direction uses a nonlinear transformation. This transformation stretches higher latitudes, causing the well-known distortion where Greenland appears larger than Africa.

export function lonLatToMercator(
  lon: number, lat: number
): MercatorCoord {
  const lam = degToRad(lon);
  const phi = degToRad(clampLatitude(lat));
  return {
    x: lam * WGS84.a,
    y: WGS84.a * Math.log(
      Math.tan(Math.PI/4 + phi/2)
    )
  };
}

Inverse Transformation

The latitude calculation for the inverse transformation is based on the following identity:

export function mercatorToLonLat(
  mx: number, my: number
): { lon: number; lat: number } {
  const lam = mx / WGS84.a;
  const phi = 2 * Math.atan(
    Math.exp(my / WGS84.a)
  ) - Math.PI / 2;
  return {
    lon: radToDeg(lam),
    lat: radToDeg(phi)
  };
}

Why Is +/-85.05 Degrees the Limit?

In the Mercator Y coordinate formula y = ln(tan(pi/4 + phi/2)), as phi approaches +/-90 degrees, y approaches +/-infinity. At the poles, the projection diverges to infinity.

In tile map systems, the entire map must fit within a square. The X range in Mercator coordinates is [-pi*a, pi*a]. To find the latitude where the Y range also fits within the same extent:

In other words, +/-85.05 degrees is the boundary latitude where the Mercator map forms a square. The interior of Antarctica falls outside the rendering area, but the practical impact is limited.

export const MAX_MERCATOR_LATITUDE =
  85.051_128_78;

export function clampLatitude(
  lat: number
): number {
  return Math.max(
    -MAX_MERCATOR_LATITUDE,
    Math.min(MAX_MERCATOR_LATITUDE, lat)
  );
}

Normalized Mercator Y

For latitude interpolation in the tile coordinate system (used in Chapter 5), we provide a "normalized Mercator Y" function without the WGS84 scaling. The difference from lonLatToMercator is that it does not multiply by a (the equatorial radius).

export function latToMercatorY(
  lat: number
): number {
  const phi = degToRad(clampLatitude(lat));
  return Math.log(
    Math.tan(Math.PI/4 + phi/2)
  );
}

When generating spherical meshes for tiles (Chapter 5), we need to "linearly interpolate latitude in Mercator Y space within the tile extent," and scaling is unnecessary for this purpose. Tile images are Mercator-projected, so the pixel Y coordinate of a tile image is linear with respect to Mercator Y. If we were to subdivide at equal intervals in latitude space, the texture would stretch at higher latitudes, causing misalignment with the tile image.

Tile Coordinate System (Slippy Map Tiles)

The tile coordinate system widely used in web maps is called Slippy Map Tiles and was standardized by OpenStreetMap. The basic mechanism works as follows:

• Treat the Mercator projection of the entire Earth as a square
• At zoom level z, divide this square into a 2^z x 2^z grid
• Assign (x, y, z) coordinates to each grid cell
• Serve a 256x256 pixel image for each cell

The tile coordinate origin (0,0) is at the top-left, with x increasing to the right (east) and y increasing downward (south).

Converting Longitude/Latitude to Tile Coordinates

The X coordinate calculation is intuitive. Normalize the longitude to the 0-360 range, then distribute it evenly among the 2^z tiles.

The Y coordinate is more involved because it goes through the Mercator projection. Note that Y=0 is at the north.

export function lonLatToTile(
  lon: number, lat: number, zoom: number
): TileCoord {
  const latRad = degToRad(
    clampLatitude(lat)
  );
  const n = 2 ** zoom;
  const x = Math.floor(
    ((lon + 180) / 360) * n
  );
  const y = Math.floor(
    ((1 - Math.log(
      Math.tan(latRad)
      + 1/Math.cos(latRad)
    ) / Math.PI) / 2) * n
  );
  return {
    x: Math.max(0, Math.min(n-1, x)),
    y: Math.max(0, Math.min(n-1, y)),
    z: zoom
  };
}

Inverse Conversion: Tile Coordinates to Longitude/Latitude

Returns the longitude and latitude of the top-left corner of tile (x, y, z). To get the bottom-right corner, pass (x+1, y+1, z). The Math.sinh used in the latitude inverse transformation, where sinh(x) = (exp(x) - exp(-x))/2, provides a concise expression for the Mercator Y to latitude conversion.

export function tileToLonLat(
  x: number, y: number, z: number
): { lon: number; lat: number } {
  const n = 2 ** z;
  const lon = (x / n) * 360 - 180;
  const latRad = Math.atan(
    Math.sinh(
      Math.PI * (1 - (2*y) / n)
    )
  );
  return { lon, lat: radToDeg(latRad) };
}

Tile Bounding Box

The geographic extent of a tile can be obtained by calling tileToLonLat twice. This TileBounds is used frequently in subsequent chapters: generating spherical meshes for tiles (Chapter 5), transforming MVT tile geometry, and more.

export function tileBounds(
  x: number, y: number, z: number
): TileBounds {
  const nw = tileToLonLat(x, y, z);
  const se = tileToLonLat(x+1, y+1, z);
  return {
    west: nw.lon,  east: se.lon,
    north: nw.lat, south: se.lat
  };
}

Tile URL Convention

Tile servers follow the convention of serving tile images at the z/x/y path. By defining a TileUrlGenerator type and extracting it as a function, switching between tile sources becomes straightforward.

export type TileUrlGenerator = (
  z: number, x: number, y: number
) => string;

// OpenStreetMap tile URL
export const osmTileUrl:
  TileUrlGenerator = (z, x, y) =>
  `https://tile.openstreetmap.org/$
    {z}/${ x}/${ y}.png`;

// GSI (Geospatial Information Authority of Japan) standard tile URL
export const gsiTileUrl:
  TileUrlGenerator = (z, x, y) =>
  `https://cyberjapandata.gsi.go.jp/
    xyz/std/${ z}/${ x}/${ y}.png`;

Drawing Tile Boundaries on the Sphere

When drawing tile boundaries on a sphere, connecting edges with straight lines would cause deviations due to the curvature of the sphere, so each edge is sampled at 20 points and approximated as a curve.

const bounds = tileBounds(tx, ty, z);

// North edge: sample 20 points from west to east
for (
  let i = 0; i <= EDGE_SAMPLES; i++
) {
  const lon = bounds.west
    + (bounds.east - bounds.west)
    * (i / EDGE_SAMPLES);
  points.push(
    geoToSceneVec3(
      bounds.north, lon, TILE_ALTITUDE
    )
  );
}

Color Coding and Z-Fighting Avoidance

To distinguish adjacent tiles, colors are assigned from an 8-color palette using (tx + ty) % 8. The checkerboard-like color variation makes tile boundaries visually clear.

Tile lines are drawn not on the ellipsoid surface but 50 km above it (TILE_ALTITUDE = 50000). Placing them at exactly the same position as the ellipsoid surface would result in identical depth values and unstable rendering order, a phenomenon known as Z-fighting.

Tile Labels

When z is 3 or less, a "z/x/y" label is displayed at the center of each tile. A Sprite (a billboard that always faces the camera) is generated using a Canvas texture and placed 30 km above the tile center. At z >= 4, labels are hidden because the number of tiles becomes too large.

Summary

In this chapter, we implemented the Web Mercator projection and tile coordinate system.

• Mercator transformation: x = a*lam, y = a*ln(tan(pi/4 + phi/2)) and its inverse
• Valid latitude range: +/-85.05 degrees (the boundary where the Mercator map forms a square)
• Normalized Mercator Y: an unscaled version for tile texture mapping
• Slippy Map Tiles: dividing the Earth into a 2^z x 2^z grid at zoom level z
• lonLatToTile / tileToLonLat: bidirectional conversion via the Mercator projection
• tileBounds: geographic extent of a tile (west/east/south/north)

In the next chapter, we will generate spherical meshes from these tile coordinates and implement the process of applying tile images as textures.