Optimizing Tile Generation for Global Web Maps

Generating map tiles that span the entire globe creates a predictable computational bottleneck. The Web Mercator projection divides the Earth into a hierarchical pyramid where each zoom level quadruples the tile count. At zoom level zero, only four tiles exist. By zoom level twelve, that number exceeds sixteen million. Rendering every tile for a global raster dataset without optimization will exhaust memory, trigger redundant disk I/O, and stall production pipelines. The solution relies on spatial filtering, windowed reading, and parallel execution. This guide provides a direct, production-ready approach to generating only the tiles that contain actual data.

Efficient tile generation is a foundational component of modern Geospatial Visualization & Web Mapping workflows, bridging raw spatial data and interactive delivery. Before rendering a single pixel, the pipeline must verify whether a target tile intersects the source dataset. Global rasters rarely contain valid data across every ocean, ice sheet, or uninhabited desert. By calculating the geographic bounds of each tile and comparing them against the dataset extent, you can skip empty regions entirely. This spatial pruning eliminates millions of unnecessary read operations and transforms an intractable task into a streamlined operation.

Production-Ready Implementation

The following script uses rasterio for geospatial I/O, morecantile for Web Mercator tile mathematics, and concurrent.futures for safe multithreading. It opens the raster inside each worker to guarantee thread safety, calculates exact tile bounds, reads only intersecting windows, and skips empty or nodata regions. The per-tile pipeline, including the spatial-pruning gates that skip empty tiles, is shown below.

flowchart TD
    A["Build tile queue<br/>for dataset extent"] --> B["Worker: open raster<br/>+ compute tile bounds"]
    B --> C{Window intersects<br/>dataset?}
    C -->|No| D["Skip tile"]
    C -->|Yes| E["Read intersecting<br/>window"]
    E --> F{Any valid<br/>non-nodata pixels?}
    F -->|No| D
    F -->|Yes| G["Normalize to 8-bit"]
    G --> H["Write z/x/y.png"]

Install dependencies: pip install rasterio morecantile numpy

import os
import numpy as np
import rasterio
from rasterio.windows import from_bounds
from rasterio.enums import Resampling
from morecantile import tms as default_tms, Tile
from concurrent.futures import ThreadPoolExecutor, as_completed

INPUT_RASTER = "global_dataset.tif"
OUTPUT_DIR = "tiles"
ZOOM_LEVEL = 6
MAX_WORKERS = 8

def process_tile(tms, tile, input_path, output_dir):
    # Open raster per-thread to avoid GIL contention and handle locks safely
    with rasterio.open(input_path) as src:
        tile_bounds = tms.xy_bounds(tile)

        # Create a window from tile bounds and clip to dataset extent
        window = from_bounds(*tile_bounds, src.transform)
        window = window.intersection(src.window(*src.bounds))

        # Skip if window falls completely outside the raster
        if window.width <= 0 or window.height <= 0:
            return

        # Read only the intersecting pixels
        data = src.read(window=window, resampling=Resampling.bilinear)

        # Identify valid data (exclude NaN and explicit nodata values)
        valid_mask = ~np.isnan(data)
        if src.nodata is not None:
            valid_mask &= (data != src.nodata)

        # Skip tiles with no valid data
        if not np.any(valid_mask):
            return

        # Normalize single-band data to 8-bit for standard web display
        if data.shape[0] == 1:
            data = data[0]
            vmin, vmax = np.nanmin(data), np.nanmax(data)
            if vmax != vmin:
                data = ((data - vmin) / (vmax - vmin) * 255).astype(np.uint8)
            else:
                data = np.zeros_like(data, dtype=np.uint8)
            data = data[np.newaxis, :, :]

        # Construct standard web tile directory structure: z/x/y.png
        out_path = os.path.join(output_dir, str(ZOOM_LEVEL), str(tile.x), f"{tile.y}.png")
        os.makedirs(os.path.dirname(out_path), exist_ok=True)

        with rasterio.open(
            out_path, 'w', driver='PNG',
            height=data.shape[1], width=data.shape[2],
            count=data.shape[0], dtype=data.dtype
        ) as dst:
            dst.write(data)

def generate_optimized_tiles():
    tms = default_tms.get("WebMercatorQuad")
    os.makedirs(OUTPUT_DIR, exist_ok=True)

    # Quick bounds check to calculate tile range
    with rasterio.open(INPUT_RASTER) as src:
        left, bottom, right, top = src.bounds
        if src.crs != tms.crs:
            from rasterio.warp import transform_bounds
            left, bottom, right, top = transform_bounds(src.crs, tms.crs, left, bottom, right, top)

        # Calculate the tile indices at the dataset corners. Tile rows (y)
        # increase southward, so the bottom edge yields the larger y index.
        ll_tile = tms.tile(left, bottom, ZOOM_LEVEL)
        ur_tile = tms.tile(right, top, ZOOM_LEVEL)

        # Clamp to valid matrix dimensions, ordering each axis correctly
        matrix = tms.matrix(ZOOM_LEVEL)
        x_min = max(0, min(ll_tile.x, ur_tile.x))
        x_max = min(matrix.matrixWidth - 1, max(ll_tile.x, ur_tile.x))
        y_min = max(0, min(ll_tile.y, ur_tile.y))
        y_max = min(matrix.matrixHeight - 1, max(ll_tile.y, ur_tile.y))

        # Build tile queue
        tiles = [
            Tile(x=x, y=y, z=ZOOM_LEVEL)
            for x in range(x_min, x_max + 1)
            for y in range(y_min, y_max + 1)
        ]

    print(f"Processing {len(tiles)} tiles across {MAX_WORKERS} workers...")

    # Execute in parallel
    with ThreadPoolExecutor(max_workers=MAX_WORKERS) as executor:
        futures = {executor.submit(process_tile, tms, tile, INPUT_RASTER, OUTPUT_DIR): tile for tile in tiles}
        for future in as_completed(futures):
            tile = futures[future]
            try:
                future.result()
            except Exception as e:
                print(f"Error processing tile {tile.x}/{tile.y}: {e}")

if __name__ == "__main__":
    generate_optimized_tiles()

Debugging & Troubleshooting

Symptom Likely Cause Resolution
CRSError or TransformError Source raster CRS differs from Web Mercator and bounds aren’t transformed Ensure transform_bounds runs before tms.tile() calculation. Verify src.crs is not None.
High memory usage / OOM MAX_WORKERS too high for available RAM, or reading full-resolution windows Reduce MAX_WORKERS to os.cpu_count() // 2. Add out_shape to src.read() if downscaling is acceptable.
All output tiles are blank Nodata values not recognized, or normalization fails on uniform data Check src.nodata. If data is already 0–255, remove the normalization block. Use np.nanmin/np.nanmax safely.
PermissionError on write Output directory locked or missing parent folders Verify os.makedirs(..., exist_ok=True) runs before writing. Run script with appropriate filesystem permissions.
Slow execution despite parallelism Disk I/O bottleneck or raster stored on network drive Switch to SSD storage. Compress source raster to BLOCKSIZE=256 and TILED=YES using gdal_translate before processing.

For detailed parameter tuning, consult the official Rasterio Documentation and the morecantile repository. When integrating these optimized tiles into static basemaps, consider how they align with pre-rendered layers in workflows like Static Mapping with Matplotlib and Contextily. Adhering to the OGC Web Map Tile Service Standard ensures compatibility with downstream mapping libraries and CDN caching strategies.