TerrainRGB, Terrarium Plugins: Further Improvements

[gkjohnson]

Related to #1690

Issues

• The current implementation does not handle raycasting / levels so well or flexibly, resulting in skipped LoDs. Instead of trying to expand children into the hierarchy it would be best to create multiple sub-meshes per tile and run per-mesh bounding box checks to accelerate raycasting
• Adjusting the height scale means every bounding box must be expanded. It would be best to do this as-needed or as-traversed, instead, with a kind of dirty check.
• The height altitude should be in a consistent unit and / or scale relative to the scale of the stretched scale of the world around it (either consistent based on camera position or scaled based on locality)

Functionality / Perf

• Core bounding-volume update infrastructure in the tiles renderer (grow-only bounds, animation, parent-child consistency, as-traversed dirty updates instead of full traversal on heightScale change)
• Debug bound visualizations refreshing after heightScale changes (falls out of the above)
• Improve accounting of downloaded DEM memory (currently double counted on sub meshes with shared items)
• Avoid regenerating a mesh on raycast
• Only update bounds when traversed or enable dynamic bound checking to avoid perf hit when updating height scale

Effects

• Hillshade-style shading effect
• Concavity map revisit (dynamically generated from DEM)

Other

• SkirtedPlaneGeometry adoption in GeneratedSurfacePlugin
• Increase skirt length?
• Upstreaming the bump map chunk fixes to three.js
• Frustum culling support (likely unnecessary unless users request it since it's disabled by default - requires updating bounds)

[jo-chemla]

New Terrain visualization modes used on terrain-viewer, based on @maplibre shared on X here

Easy additions for viz modes would probably be the ones most people use and which are low-hanging fruits: the derivatives-based terrain modes:

• slope: magnitude of the gradient, color-ramped with editable bounds see eg cpt-city colorramps
• profile and planar curvatures: 2nd order derivs, representing flow acceleration or divergence via a diverging color-ramp
• and aspect: angle of the gradient, a colored representation of the surface normal azimuth (pointing east...), usually hue-coded

Formulas based on GDAL Horn approximation of the gradient in this docs page, also in details block below

Formulas

Shared inputs

Elevation decoding (lib/elevation-encoding.ts) — every mode starts by decoding an upstream raster-dem tile:

$$h_{\text{terrarium}} = \left(R \cdot 256 + G + \frac{B}{256}\right) - 32768$$

$$h_{\text{mapbox}} = -10000 + (R \cdot 256^2 + G \cdot 256 + B) \cdot 0.1$$

The Horn 3×3 gradient (hornGradient(), ported from GDAL's GDALSlopeHornAlg), the shared building block for Slope, Aspect, TRI*, Curvature, and (per-sub-cell) Blobness. For a 3×3 window $a_0..a_8$ (row-major, $a_4$ = center) with ground resolution $L$ at the tile's center latitude:

$$\frac{\partial z}{\partial x} = \frac{a_0 + 2a_3 + a_6 - a_2 - 2a_5 - a_8}{8 L \cos(\text{lat})}$$

$$\frac{\partial z}{\partial y} = \frac{a_6 + 2a_7 + a_8 - a_0 - 2a_1 - a_2}{8 L \cos(\text{lat})}$$

^{*TRI/TPI/Roughness don't use the gradient directly — they're plain 3×3 aggregations of elevation itself, see below.}

Terrain analysis

Slope

$$\text{slope}° = \arctan\left(\sqrt{\left(\tfrac{\partial z}{\partial x}\right)^2 + \left(\tfrac{\partial z}{\partial y}\right)^2}\right) \cdot \frac{180}{\pi}$$

Re-encoded via Mapbox Terrain-RGB packing (base −10000, 0.1° step). lib/slope-protocol.ts.

Aspect

Compass bearing of the gradient direction, converted from math convention (CCW from east) to compass convention (CW from north):

$$\theta_{\text{math}} = atan2\left(\tfrac{\partial z}{\partial y},, -\tfrac{\partial z}{\partial x}\right) \cdot \frac{180}{\pi}$$

$$\text{aspect}° = (90 - \theta_{\text{math}}) \bmod 360$$

(0 if $\partial z/\partial x = \partial z/\partial y = 0$, i.e. flat ground.) lib/aspect-protocol.ts.

Curvature

Four sub-modes share one Zevenbergen & Thorne (1987) second-order fit over the 3×3 window ($p, q$ = first partials; $r, t, s$ = second partials; $L$ = ground spacing):

$$p = \frac{a_5 - a_3}{2L} \qquad r = \frac{a_5 - 2a_4 + a_3}{L^2}$$

$$q = \frac{a_7 - a_1}{2L} \qquad t = \frac{a_7 - 2a_4 + a_1}{L^2}$$

$$s = \frac{a_2 - a_0 - a_8 + a_6}{4L^2}$$

Profile (curvature along steepest descent — flow acceleration) / Plan (curvature across contours — flow convergence/divergence, $\equiv \operatorname{div}(\nabla z / |\nabla z|)$), both ×100, undefined ($\to 0$) on flat ground:

$$g = p^2 + q^2$$

$$\text{profile} = \frac{100,(r p^2 + 2spq + tq^2)}{g,(1+g)^{1.5}}$$

$$\text{plan} = \frac{100,(r q^2 - 2spq + tp^2)}{g^{1.5}}$$