Web Mercator Tiles

The internet is full of interactive tiled slippy web maps. These maps load tiles as needed in response to the user zooming and panning around, which allows users to explore extremely large datasets with minimal lagtime. This style of web map was first popularized by Google Maps in 2005, and since then the approach has been adopted by other mapping libraries, such as Leaflet.js, OpenLayers, Mapbox and many others. These interactive maps utilize an XYZ tiling scheme, where Z represents the zoom level and X and Y are denote the tile indices.

Figure 1: Overview of web mercator tiles zoom levels 0 through 5 with an example tile for each zoom level. Source: File:XYZ Tiles.png. License CC0. If you prefer an interactive 2-D map, check out MapTiler’s Tiles à la Google Maps.

These tiles use a Web (pseudo) Mercator projection (EPSG:3857), which is defined using a strange mix of ellipsoidal and spherical properties. If you care about positional accuracy, use this projection at your own risk. Also the projection strongly distorts areas. For example, compare how Greenland and Africa appear to be similar in size (2.2 million \(\text{km}^2\) and 30.4 million \(\text{km}^2\), respectively).

But if you are creating an interactive map for the web, the projection can be very useful.

Mathmatical Description

The Web Mercator projection is defined by the following formulas. The tile’s \(x\)-index is a function of longitude (\(\lambda\)) and zoom level (\(z\)): \[ x = {\frac {1}{2\pi }}\cdot 2^z\left(\pi +\lambda \right) \tag{1}\] and the tile’s \(y\)-index is a function of latitude (\(\varphi\)): \[ y = \frac{1}{2\pi } \cdot 2^z\left(\pi -\ln \left[\tan \left({\frac {\pi }{4}} + \frac{\varphi }{2}\right)\right]\right) \tag{2}\]

Rearranging Equation 1 to solve for longitude (\(\lambda\)) yields: \[ \lambda = \frac{2\pi \cdot x}{2^z} - \pi \tag{3}\]

and rearranging Equation 2 to solve for latitude (\(\varphi\)) yields: \[ \varphi = 2 \cdot \arctan(e^{\pi - \frac{2\pi \cdot y}{2^z}}) - {\frac {\pi }{2}} \tag{4}\]

Python Implementation

In Python those equations can be written as:

import numpy as np
import numpy.typing as npt

def xyz_to_lonlat(
    x: npt.NDArray[np.floating], # horizontal index
    y: npt.NDArray[np.floating], # vertical index
    z: int, # zoom level
):
    lon = 2 * np.pi * x / 2**z - np.pi
    lat = 2 * np.arctan(np.exp(np.pi - (2 * np.pi * y / 2**z))) - np.pi / 2
    return (lon * 180 / np.pi, lat * 180 / np.pi)

1

NumPy, the “fundamental package for scientific computing with Python” (https://numpy.org/)

2

NumPy type hinting (https://numpy.org/doc/1.21/reference/typing.html)

3

See Equation 3

4

See Equation 4

To understand how tiles differ between levels and across the Earth, we create a GeoDataFrame object containing the tile information.

import geopandas as gpd
import shapely

zoom_levels = [0, 1, 2, 3, 4]

# for layer in layer_config_list:
tile_data_list = []
for zoom in zoom_levels:
    tilerange = range(2**zoom + 1)

    # Create arrays of cell vertices.
    x_vertices = np.array(tilerange)
    y_vertices = np.array(tilerange)

    # Convert vertices to lon/lat coordinates.
    lons, lats = xyz_to_lonlat(x_vertices, y_vertices, zoom)

    for ix in tilerange[:-1]:
        for iy in tilerange[:-1]:
            poly_coords = [
                (lons[ix], lats[iy]),
                (lons[ix + 1], lats[iy]),
                (lons[ix + 1], lats[iy + 1]),
                (lons[ix], lats[iy + 1])
            ]
            # print('poly_coords', poly_coords)

            tile_data = {
                'zoom': zoom,
                'x_index': ix,
                'y_index': iy,
                'geometry': shapely.Polygon(poly_coords),
            }
            tile_data_list.append(tile_data)

gdf = gpd.GeoDataFrame(
    data=tile_data_list,
    crs="EPSG:4326"
)

Based on code described in this earlier post, we can calculate the area covered by each tile.

import math

def patch_area_single_rectangle(
    R: float, # Authalic radius, in meters
    lon1: float, # Minimum longitude, in degrees
    lat1: float, # Minimum latitude, in degrees
    lon2: float, # Maximum longitude, in degrees
    lat2: float, # Maximum latitude, in degrees
) -> float:
    # Convert from degrees to radians
    lon1_rad = lon1 * (np.pi / 180)
    lon2_rad = lon2 * (np.pi / 180)
    lat1_rad = lat1 * (np.pi / 180)
    lat2_rad = lat2 * (np.pi / 180)
    
    h = R * (math.sin(lat2_rad) - math.sin(lat1_rad))
    area_patch = R * h  * (lon2_rad - lon1_rad)
    return area_patch

gdf['area_km2'] = gdf.apply(
    lambda row:patch_area_single_rectangle(
        6371007.180918,  # Authalic radius, in meters
        *row['geometry'].bounds,  # tuple of (minx, miny, maxx, maxy)
    ) * 1e-6,
    axis=1
)

Let’s see a summary of this tile dataset.

gdf

	zoom	x_index	y_index	geometry	area_km2
0	0	0	0	POLYGON ((-180 85.05113, 180 85.05113, 180 -85...	5.081641e+08
1	1	0	0	POLYGON ((-180 85.05113, 0 85.05113, 0 0, -180...	1.270410e+08
2	1	0	1	POLYGON ((-180 0, 0 0, 0 -85.05113, -180 -85.0...	1.270410e+08
3	1	1	0	POLYGON ((0 85.05113, 180 85.05113, 180 0, 0 0...	1.270410e+08
4	1	1	1	POLYGON ((0 0, 180 0, 180 -85.05113, 0 -85.051...	1.270410e+08
...	...	...	...	...	...
336	4	15	11	POLYGON ((157.5 -55.77657, 180 -55.77657, 180 ...	1.439368e+06
337	4	15	12	POLYGON ((157.5 -66.51326, 180 -66.51326, 180 ...	7.045811e+05
338	4	15	13	POLYGON ((157.5 -74.01954, 180 -74.01954, 180 ...	3.321439e+05
339	4	15	14	POLYGON ((157.5 -79.17133, 180 -79.17133, 180 ...	1.537909e+05
340	4	15	15	POLYGON ((157.5 -82.67628, 180 -82.67628, 180 ...	7.061717e+04

341 rows × 5 columns

Among other things, the data table shows that tile areas vary by y_index and zoom level. But it is difficult to grasp how these values change around the Earth from a table.

Tile Visualization

This section creates visualizations of Web Mercator tiles for various zoom levels. We start by formatting the tile data so that can be used by PyDeck.

tile_colors = [
    [0, 0, 0],
    [255, 255, 255],
]

def add_pydeck_colums(gdf):
    """Add columns for PyDeck visualization"""
    
    # Convert geodataframe geometries into a PyDeck compatible format
    gdf['pydeck_polygon'] = [
        [
            vertex + (200000, )  # Add a third dimension
            for vertex in list(geom.exterior.coords)
        ]
        for geom in gdf['geometry'].geometry
    ]
    gdf['pydeck_centroid'] = [
        list(geom.centroid.coords[0]) + [800000] # Add a third dimension
        for geom in gdf['geometry'].geometry
    ]
    gdf['centroid_label'] = '(' + gdf['x_index'].astype(str) + ',' + gdf['y_index'].astype(str) + ')'
    gdf['bounds'] = [
        [f'{x:.2f}' for x in geom.bounds]
        for geom in gdf['geometry'].geometry
    ]

    gdf['fill_color'] = gdf.apply(
        lambda row: tile_colors[(row['x_index'] + row['y_index']) % 2],
        axis=1
    )

    gdf['area_desc'] = gdf.apply(
        lambda row: f"{row['area_km2']:.2e}", 
        axis=1
    )

    return gdf

gdf = add_pydeck_colums(gdf)

We then create a function that can generate PyDeck objects for a specified zoom level.

import pydeck as pdk
from pydeck.types import String
import pandas as pd

def get_globe_view():
    """Return a pyDeck view for the globe."""
    view = pdk.View(
        type="_GlobeView",
        controller=True,
        width='100%',
    )
    return view

def get_initial_view_state():
    """Specify the initial pyDeck view state."""
    initial_view_state = pdk.ViewState(
        latitude=45,
        longitude=-45,
        zoom=0,
        min_zoom=0,
    )
    return initial_view_state

def get_basemap_layer():
    """Get a pyDeck layer based on Natural Earth country boundaries."""
    basemap_layer = pdk.Layer(
        "GeoJsonLayer",
        id="base-map",
        data="https://d2ad6b4ur7yvpq.cloudfront.net/naturalearth-3.3.0/ne_50m_admin_0_scale_rank.geojson",
        stroked=False,
        filled=True,
        extruded=False,
        get_fill_color=[200, 200, 200],
        get_line_color=[255, 255, 255],
        opacity=0.5,
    )
    return basemap_layer

def get_tiles(zoom):
    """Create a pyDeck polygon and text layers for tiles of a specified zoom level."""
    poly_layer = pdk.Layer(
        "SolidPolygonLayer",
        gdf[gdf['zoom']==zoom],
        opacity=0.25,
        visible=True,
        get_polygon="pydeck_polygon",
        get_fill_color="fill_color",
        extruded=True,
    )

    text_layer = pdk.Layer(
        "TextLayer",
        gdf[gdf['zoom']==zoom],
        get_position="pydeck_centroid",
        get_text="centroid_label",
        get_size=12,
        # get_color=[127, 255, 255],
        get_color=[255, 255, 0],
        get_angle=0,
        get_text_anchor=String("middle"),
        get_alignment_baseline=String("center"),
        pickable=True,
        extruded=True,
    )
    return [poly_layer, text_layer]


def create_globe_with_xyz_tiles(zoom):
    
    # gdf = layer_config['gdf']
    deck = pdk.Deck(
        views=[get_globe_view()],
        initial_view_state=get_initial_view_state(),
        layers=[
            get_basemap_layer(),
            get_tiles(zoom),
        ],
        map_provider=None,
        tooltip={
            "html": (
                "<table>"
                "<tr><th>Zoom:</th><td>{zoom}</td></tr>"
                "<tr><th>X, Y:</th><td>{x_index}, {y_index}</td></tr>"
                "<tr><th>Area, km<sup>2</sup>:</th><td>{area_desc}</td></tr>"
                "<tr><th>Longitude:</th><td>{bounds.0}&deg; to {bounds.2}&deg;</td></tr>"
                "<tr><th>Latitude:</th><td>{bounds.1}&deg; to {bounds.3}&deg;</td></tr>"
                "</table>"
            ),
            "style": {
                "text-align": "right",
            }
        },
        # Note: set to True for the globe to be opaque
        parameters={"cull": True},
    )
    return deck

The following figures show how the Web Mercator tiles are distributed on the globe for zoom levels 1 through 4. Viewing Web Mercator tiles this way clearly shows that the tiles do not reach the poles, regardless of the zoom level. Rather the tiles end at latitudes \(85.051\degree\) N and \(85.051\degree\) S¹.

Tiles - Zoom Level 1

For zoom level 1, the x and y tile indices range from 0 to 1². There are a total of 4 tiles³.

create_globe_with_xyz_tiles(1)

Tiles - Zoom Level 2

For zoom level 2, the x and y tile indices range from 0 to 3⁴. There are a total of 16 tiles⁵.

create_globe_with_xyz_tiles(2)

Tiles - Zoom Level 3

For zoom level 3, the x and y tile indices range from 0 to 7⁶. There are a total of 64 tiles⁷.

create_globe_with_xyz_tiles(3)

Tiles - Zoom Level 4

For zoom level 4, the x and y tile indices range from 0 to 15⁸. There are a total of 256 tiles⁹.

create_globe_with_xyz_tiles(4)