Although I don’t work in mapping anymore, I’ll always have have a soft spot for experimental map rendering techniques. Previously I’ve built a live solar terminator map and a kinda janky 3D map renderer. Inspired by some good posts on the topic, I’ve also dabbled in heightfield shadowing. This article finally pulls the latter work together into something functional.

Our goal here is to compute terrain lighting using standard web map terrain tiles as input and producing corresponding 1:1 raster light maps as output. For terrain data, we’ll be using the excellent and open Mapterhorn tileset. We’ll start with a simple horizon visibility algorithm and extend to cover self-shadowing and ambient occlusion. It’s definitely not novel, but the devil is in the details.

Inspiration

First of all, the best resource on the topic is Tom Patterson’s shadedrelief.com. Any work in map shading could only aspire to convey topography and tell stories in the way his maps do. Alas, I’ve set myself up to fail in the shadow of such a behemoth. But at any rate, my goal is somewhat focused here, to shade map tiles for real time interactive usage.

Lighting in computer graphics is usually either complex or expensive. One approach for heightfield shading is to shoot rays. Rye Terrell’s terrain shading marches rays from each pixel along the surface and in the direction of the sun, checking for occlusion along the way. The results are beautiful, but the cost is high and it suffers from tile-boundary discontinuities which Terrell addresses by brute-force rendering 3×3 tiles (768×768) and extracting the center.

A cheaper alternative is to sweep across the grid, tracking horizon visibility in a stack. Karim Naaji’s Line-Sweep Ambient Occlusion sweeps rows of the grid, tracking visibility and applying shading in a single pass. However, the technique is limited to eight cardinal directions and still suffers from tile boundary discontinuities.

This article doesn’t offer fundamentally new insight but instead builds on and extends the line-sweep idea to address their shortcomings. In particular it extends the line-sweep to arbitrary directions, uses geometric shadowing to produce soft shadows, and avoids tile boundary discontinuities by prewarming horizon visibility with parent tile data.

Line-sweep horizon visibility

The core primitive behind this notebook is a 1-D sweep that tracks the upper convex hull of elevations seen so far along each ray. Imagine sweeping along a single slice of terrain elevations. As the sweep advances left to right, we test each new point to see if it’s a horizon. If we step uphill far enough that multiple horizons become visible, we pop horizons off of the stack until exactly one horizon remains visible.

Self-shadowing

From here, it’s only a small step farther to compute self-shadowing. At each point during the sweep, we compare the horizon angle to the sun angle. If the sun is below the horizon, the terrain is shadowed.

A particularly nice feature of this algorithm is that we can compute horizon visbility and lighting in the same pass. Neglecting the depth of the stack which generally remains quite small, the algorithm is O(n) in the resolution of terrain data, requiring that we read each terrain value only once. I can’t immediately see that we could do significantly better.

A little more work is necessary before we can apply this to terrain tiles though. The algorithm so far presumes that we’re sweeping along a well-defined grid direction. We’d like to be able to compute shadowing in arbitrary directions.

The solution is simple but requires some finesse. The diagram below illustrates grid traversal in arbitrary directions. We start our first sweep in the upper left corner, seeding subsequent sweeps with remaining pixels along the left edge. For each sweep, we step rightward exactly one pixel and upward by a fractional amount. Once we’ve exhausted all sweeps starting along the left edge, we step rightward along the bottom edge with a given stride, performing identically the same sweeps from the new seed points.

This sequence only handles sweep angles in a single octant of the possible Drag the handle below and observe that for other octants, we simply flip and rotate our canonical octant. In practical terms, this can be implemented with ndarray-style index arithmetic.

One subtlety is worth calling out. Rays sweeping at non-axis-aligned angles pass between integer pixel rows. Nearest-neighbor assignment drops a ray’s full weight onto a single row, so adjacent rays can land on the same row repeatedly and leave others starved, producing visible banding. Distributing each sample’s weight between the two neighboring rows by , where is the fractional row offset, evens out the coverage. The toggle on the figure below compares the two schemes.

A final note is that since sampled elevations are directly adjacent to pixels receiving shading, even small errors in sample position result in errors in computed slope angle. Any inconsistencies causes glancing sunlight to produce terrible striations on smooth hillsides.

The problems can be overcome though, and the figure below illustrates the result. Drag the handle and observe that the terrain tile correctly casts shadows on itself.

Soft shadows

We’ve been considering hard thresholds, but we’re well poised to relax this to a smooth threshold so that we can emulate the soft penumbra of a finite-sized light source like the sun. Realistically, the sun covers about degrees in the sky, making it pretty close to a point source, but we can always increase the size for visual effect.

The figure below illustrates a soft shadowing pass using a smoothstep function for horizon visibility occlusion.

Sadly though, the sun is a two-dimensional source, and we’ve only discussed softening in one dimension. Short of something fancy like mip-maps, all I was able to come up with was to naively average multiple samples of our altitude-wide analytical soft shadows across the azimuthal direction.

Samples aren’t cheap, so it’s worth thinking carefully about how to distribute samples across the sun disk. Our altitude-wise analytical soft shadows use a smoothstep function, so after some consideration, I scaled a Gaussian function to roughly approximate the smoothstep function, then distributed equally-weighted samples in the azimuthal direction using inverse-CDF sampling to approximate the desired Gaussian.

The sun isn’t really a Gaussian disk, of course. But Gaussians are mathematically separable, which makes the 1D altitude × 1D azimuth factorization work out cleanly, and the resulting penumbra looks smooth enough in practice that the approximation pays off.

Parent-tile pre-pass

A natural limitation of per-tile methods is that they have no knowledge of terrain beyond it. Mountains just outside the boundary which would cast shadows onto the tile are never considered, producing visible seams between tiles.

A natural solution is to generate a half-resolution horizon buffer from four z−1 ancestor tiles, covering the computation tile and one tile width of surrounding terrain in each direction. Before the main sweep, a pre-pass walks each ray through this coarser buffer, warming up the stack so that by the time the ray enters the computation tile, it already carries the horizon contributed by surrounding terrain.

We could also step multiple levels upward in the tile pyramid to extend the range of this pre-warming. A 2×2 block of parents physically covers comp-tile widths around the computation tile (3 tiles at , 5 at , 9 at ), so bumping by one doubles the reachable blocker radius at the cost of halving the horizon’s parent-pixel resolution. In practice, distant blockers matter less after the / smoothstep falloff, so the coarser horizon tends to be a near-free quality win.

Toggle the parent-tile pre-pass below and observe its effect in the vicinity of tile boundaries. The selector shows how the covered-area / resolution trade plays out visually — at the comp-tile red outline shrinks to of the horizon width (4× the surrounding area captured), and at to .

Line-sweep ambient occlusion

The same horizon visibility sweep works for ambient occlusion as well. Ambient occlusion proceeds similarly to shadowing, but instead of illumination by a point source, each pixel is illuminated by the portion of sky visible to it. I extended Karim Naaji’s LSAO to support parent-tile occlusion and arbitrary sweep directions. That leaves one design knob open: the per-slice visibility function that maps each swept horizon angle to a shade contribution.

Naaji uses . The function drops with slope right at , so even a shallow horizon produces an immediate darkening. The tradeoff is a long tail, so distant tall features keep contributing, and tile boundaries take extra care.

A Lambertian sky integral offers a physically motivated alternative. Incident sky radiance from a direction at zenith angle is weighted by , so the visible fraction of one azimuthal slice above a horizon at elevation α is

which normalises to : 1 at a flat horizon, 0 at a vertical wall. Its slope is flat at and steepens later, giving a softer, more forgiving response to small horizons. And because for small α, distant tall blockers fall off as , so the effect is strongly local and tile boundaries are less visible.

Honestly, I think the particular choice of function is best thought of as an artistic decision. The original function seems to produce better visual separate of features, at the cost of less locality.

Finally, a contrast slider tone-maps the response. For pixel let and ; an tonemap squashes the scaled shade into :

At the slider emits solid white; as the scaled shade runs to infinity and the tonemap saturates toward full darkness. Naaji’s reference fires eight axis-aligned passes; because our canonical reparameterization handles arbitrary azimuths, we can use any number of directions, and each pass is independent so we run them in the same Web Worker pool as the soft shadows.

Relief shading

The shadow map captures shadows, but even fully lit surfaces vary in brightness depending on their orientation. The most natural lighting choice is diffuse Lambertian shading, but this tends to darken the surface when lit at oblique angles, which isn’t necessarily what we want for a map.

Relief shading separates the terrain slope from its aspect (direction it faces relative to the sun), so the sun’s altitude doesn’t darken flat ground. Flat terrain stays white; slopes facing the sun brighten, slopes facing away darken, and the strength of the effect is proportional to the slope angle.

Concretely, I’ve slightly adapted GDAL’s relief shading function, coputing the slope angle and the aspect alignment , then a darkness term where is the strength parameter. The final relief value is . The z-factor (where is the sun’s altitude) dampens gentle slopes following the GDAL convention, and modulates the perceived slope by sun altitude: at overhead sun the terrain appears flat and the image is white; as the sun lowers, slopes become more visually pronounced, reaching full contrast at the horizon.

WebGPU implementation

The CPU path’s two modes (soft shadow with stratified azimuth integration, and LSAO) port cleanly to WGSL compute via dynamic shader specialization. Each (mode, prewarm) permutation is built once into its own GPUComputePipeline and cached, so the inner hull-scan loop has no per-pixel mode branch. The parent-tile prewarm is itself a compile-time flag: when off, the shader skips the warmup block entirely and the horizon binding shrinks to a 16-byte stub. Multi-sample averaging (azimuth offsets for soft, uniform directions for LSAO) is dispatched as N successive compute passes that each accumulate weight · bit · fpScale into the shared atomic shadow buffer.

A caveat from building this out is that the GPU path turned out to be less useful than hoped. It’s measurably faster per tile than the CPU worker pool, but only modestly so, and the same device also drives the map’s own render passes. Every heavy compute submit stalls the tile viewer’s compositor and makes pans and zooms feel gummy. The CPU worker pool runs slightly slower per tile but leaves the GPU free to stay responsive to pointer events, which matters more than raw throughput for an interactive viewer. Both paths are kept here for comparison, and the live viewer at the end of the notebook defaults to the CPU bake.

Zoom-dependent gradient attenuation

One annoyance on a slippy map is that terrain visibly flattens at low zoom. Peaks survive downsampling well, but the gradient field between them does not. Decimation fails to commute with the nonlinearities in between (the resampling kernel, the tonemap, the slope calculation itself), so low-zoom tiles systematically under-report slope. Cast shadows read and hold up fine, while relief and LSAO read and go flat.

A quick probe across a handful of Mapterhorn summits suggests a rough power-law attenuation with trending toward . We apply a hand-wavy corrective factor only to the passes, contracting pxSize rather than inflating heights so that the elevation buffer stays untouched:

A theoretical value of matches the asymptote but reads visually overeager, so the viewer’s zoom-compensation slider defaults to . Setting it to disables the correction.

Live WebGPU tile viewer

Finally, we have all the pieces in place to produce a map. In the map below, each tile goes through a two-phase bake (CPU or GPU). The static bake runs once when the tile’s elevation data arrives: compute surface normals, 8-pass LSAO, and pack the results into an rgba8unorm texture:

The shadow re-bake runs whenever the sun direction changes, packing shadow values into the alpha channel for use in the final shading step.

The map viewer below is thrown together hastily with the help of Claude since I was not able to find another viewer which satisfied my requirements. Although it’s not hard to get tiles onto a map, it’s tricky if you really care about tile lifecycle and smooth updates. Pan with click+drag; zoom with the mouse wheel.