From a DEM to Estrat

Once the selected region is chosen in the app, it is enclosed in a bounding box that becomes the output grid. For each cell of that grid, Estrat computes its position in the Web Mercator tile grid. It then uses each of those geographic positions to sample the elevation and style tiles and the rasterized boundary.

The elevation tiles use Terrarium RGB encoding, which the CPU decodes with the following formula per each pixel and normalises it

height = R × 256 + G + B / 256 − 32768

CPU acquisitionshared destination grid

The boundary polygon is rasterised with an even-odd fill which gives us a binary coverage texture. A binary edge doesn't give good results because it creates an abrubt jumpt at the boundary which would result in rendering artifacts, so we turn it into a signed distance field by running a two-pass chamfer transform twice. On each transform, feature pixels start at distance 0 and every other pixel starts at Infinity. A forward, left-to-right/top-to-bottom, scan and a backward, right-to-left/bottom-to-top, scan then relaxes those values at cost 1 on cardinal neighbours and cost √2 at diagonal ones.

The first transform allows us to compute the distance to the nearest inside pixel. The second one measures the distance to the nearest outside pixel. For interior pixels we then keep the distance to the outside and for exterior pixels the negative distance to the inside:

SDF(p) =  distance to outside, p inside
       −distance to inside,  p outside

This results in an approximate Euclidean distance, not a perfect geometric SDF. Its main job is to provide a stable sign and a useful contour that can be tested during ray traversal.

SDF: exterior < 0zero contour = boundaryinterior > 0

After that Estrat finds sea-level pixels that touch the raster edge or the outside of the selected polygon and uses a flood algorithm to remove ocean-connected sea-level cells without removing enclosed lakes.

Estrat builds two elevation textures from the elevation source: a macro field and a detail field.

For the macro field we repeat a mask-aware separable blur pair three times to remove isolated DEM noise and preserve larger landforms. We then mix that with 15% of the source to avoid the result becoming completely flat. This becomes the authoritative surface for silhouettes, intersections and long shadows.

For the detail field we use two separable blur pairs using a 1-2-1 filter, creating a lightly smoothed detail field used for normals, short horizons, lowland cavities and elevation tint.

two height texturesone surface, two uses

Separating these frequencies keeps a noisy one-pixel DEM feature from becoming a spike or a kilometre-long shadow while preserving small structures in flat areas.

The CPU maps the selected area ground dimensions into a centered right-handed world, computing the aspect ratio from physical kilometers to avoid the view being distorted by longitude shrinking with latitude. Estrat then converts the world x,y to UV coordinates and samples the macro field as H(x,y).

Instead of just sampling the nearest pixel, Estrat loads the four surrounding values, interpolates across the top and bottom pairs and then interpolates between those two results. This bilinear lookup gives the ray marcher a continuous heigh function instead of pixel values that look like an staircase. The same kind of lookup is performed when sampling the signed-distance coverage field.

The macro field and the signed-distance coverage field answer different questions. The height texture says how high the surface is at a given position while the signed distance field says wether that position belongs to the selected area or not and how close it is to the edge. Having those two allows Estrat to stop treating the terrain as valid on boundaries during ray traversal to avoid creating vertical faces.

For an open relief, the implicit surface is simply:

F(P) = P.z − H(P.x, P.y)

Estrat draws a full-screen triangle and, for each fragment, it builds an orthographic ray from yaw, pitch and distance values. Ray origins are backed away from the exaggerated relief to avoid tall peaks being behind the camera plane.

Given the previously defined surface function F(P) and the parametrical equation of a ray, P(t), what we need to do condenses to:

find t where F(P(t)) = 0

Before sampling the DEM, the shader intersects the ray with the terrain’s axis-aligned bounding box using the slab method. This produces an entry and exit distance, tmin and tmax. Ray marching starts at max(tmin, 0) and stops at tmax with rays that miss the box returning immediately.

On ordinary rays, Estrat uses an adaptive advance. Given the current point, the shader measures the vertical separation |F(P)|, multiplies it by a conservative heuristic constant, and clamps the result between a minimum and maximum step. When the ray is far from the surface, this moves quickly. When the ray is near the surface, it takes shorter steps. The minimum prevents a nearly tangent ray from not moving while the maximum limits how many terrain cells the ray can cross in one advance.

Readers will notice this is not sphere tracing. That's a deliberate choice. |F(P)| gives us the vertical height, not the Euclidean distance to the surface and, since an oblique ray can travel a long horizontal distance while changing very little in z, using the raw separation as the step length could jump over narrow ridges. To avoid that, we conservatively allow a small number of texels per jump. Even when the vertical gap suggests a large step, this forces the ray to evaluate the terrain regularly enough that it's less likely to skip narrow ridges or valleys between samples.

After advancing along the ray, we check for a change in value signs, which means the surface lies between those two samples. Once a bracket is found, we keep refining via a fixed number of bisection steps. The final hit is the midpoint of the remaining interval, which is then used for normals, material and lightning.

Notice in the image above, that when the ray is far above the terrain, the height difference is large, so the next evaluation can be farther away. When the ray approaches the surface, the height shrinks, evaluating the height field more frequently, allowing us to buy precision near the place that matters without spending the same cost at every step.

For any hit point, the shader samples the detail field on either side of the point and forms a scale-aware central difference to compute the normal. While the macro field defines the hit, the detail field supplies the local orientation.

normal encodingR = east/west slopeG = north/south slopeB = up

Actual settled surfaceCentral differences

Lambert's cosine law is used to compute the direct light. A slope facing the sun must receive more direct light while a slope turning away receives less. Light azimuth and altitude are used to produce the unit light vector L.

Colour + orientationNo ambient terms

After computing the direct light, Estrat uses a compute pass that fills a 384x384 world-space atlas with three values: Terrain sun visibility, which answers if the sunlight can reach the terrain surface, receiver sun visibility, which answers if the sunlight can reach the background, and a broad sky openess value, which tells how much of the surrounding sky is visibile from the terrain.

384x384 was chosen because it gives us ~145.5k cells, enough resolution for broad terrain shadows and horizons, while also being divisible by the 8x8 compute workgroup size.

In order to make real valleys receive less of the sky than open plains, Estrat probes eight azimuthal horizons around every lighting-atlas cell. For each azimuth, twelve samples are placed at fixed distances. The shader then compares the height difference on each of those distances to the horizontal distance and retains the maximum positive slope. Computing the atan of that value converts the slope to a horizon angle. After that, normalising that angle against π/2 gives a simple openness fraction. Averaging the eight directions approximates the visible sky dome.

opennessdark = enclosedwhite = open hemisphere

Broad × local opennessTerrain-space

The fragment shader adds six shorter, higher-resolution horizon probes from the detail field that are evaluated only for visible hits. While the compute atlas sees broad landforms, these higher-resolution horizon probes see local drainages without having to compute long horizons for every pixel.

Finally, the shader also evaluates a compact short-range obscurance heuristic. Six azimuths sample the detail field at three quadratically increasing distances. In each direction it first derives the slope of the tangent plane from the surface normal, then measures only terrain that rises above that plane. This conforms a view-independent contact-scale relief value that attenuates the sky term while leaving the direct sunlight alone.

The hit's color is either sampled imagery, a solid colour or an hypsometric ramp that uses real elevation stops. To increase the contrast in low-lands, a two-radius residual compares the detail field with nearby means, updating the material brightness to capture local structure after the macro field has been smoothed.

The shader then combines the following values to decide the final color:

Before output, the shader also applies an ACES-like tone curve and adds a small noise to prevent banding on smooth gradients.

Estrat uses the same macro field twice:

1. During pointer movement, a 192x192 indexed grid samples it in a vertex shader and rasterises a low-resolution mesh. This gives the browser a fixed-function triangle setup instead of running a long intersection loop for every pixel.
2. When pointer movement ends, the renderer schedules an intermediate ray-cast frame and then a full one. The shadow atlas is reused for camera-only changes while terrain, exaggeration, sun direction and shadow softness invalidate it.

The final relief is not the product of a single heroic rendering trick: aligned rasters, macro and detail fields, terrain-sky visibility, etc. are all accumulated to drive the final result.

Following one Catalonia frame makes the dependencies visible. The northern Pyrenean ridge first appears as high scalar values, survives generalisation as broad geometry, acquires orientation through derivatives, blocks the finite north-west sun, closes the sky over neighbouring valleys and finally receives the upper-end of the adaptive tint. The low central and western lands barely change their silhouette, so they depend much more on detail normals, local horizons, cavity contrast and closely spaced low tint stops.

The slider below allows you to see, on a fixed camera, how each step changes the final image.

Catalonia

Detail

1 / Albedo

The adaptive hypsometric material at the primary hit, before illumination.

←→