Skip to content

McHenrys Peak Terrain Relief

This example drapes aligned NAIP orthoimagery over USGS 3DEP elevation.

Preview

McHenrys Peak Terrain Relief

Open the live WebGPU example.

Run And Adapt

Commands below assume a Datoviz source checkout and start at the repository root. Use your configured build environment; Python routes additionally require local bindings.

Route Availability Command or action
C Canonical native source just example-c showcases/terrain_relief (build and run), or rerun ./build/examples/c/showcases/terrain_relief
Python Available python3 -m examples.python.gallery.showcases.terrain_relief
Data preparation Additional integration source; check optional dependencies python3 -m tools.data.prepare_terrain_relief
Browser Live WebGPU route Open live example

Prepared data required

This example intentionally fails when its prepared input is absent; it does not substitute synthetic data. Expected input: data/examples/terrain_relief/prepared/terrain.bin; data/examples/terrain_relief/prepared/terrain.jpg. Prepare it from the repository root with uv run tools/data/prepare_terrain_relief.py.

Use this example as capability or integration evidence, not as a minimal copy-paste template. Start from the nearest supported, copy-safe example and add this feature after verifying the linked API reference.

What To Look For

The real bare-earth DEM drives a lit surface-grid mesh while the matching natural-color aerial image is sampled through the grid UVs. Rotate the arcball view to inspect McHenrys Peak, Glacier Gorge, alpine lakes, and the Continental Divide.

Prepare the cache before running:

Uv run tools/data/prepare_terrain_relief.py

Source

#!/usr/bin/env python3
"""Prepared USGS elevation with aligned NAIP imagery."""

from __future__ import annotations

import ctypes
import struct
from pathlib import Path

import numpy as np
from PIL import Image

import datoviz as dvz
from examples.python.gallery import common as ex


BUNDLE_DIRS = (Path("data/examples/terrain_relief/prepared"), Path(".cache/datoviz/examples/terrain_relief/prepared"))
HEADER = struct.Struct("<8sIIIffff")
DISPLAY_WIDTH = 5.0
EXAGGERATION = 1.35
INITIAL_ANGLES = (ctypes.c_float * 3)(+0.050954, -0.163512, -0.032719)


def _load():
    directory = next((path for path in BUNDLE_DIRS if (path / "terrain.bin").exists() and (path / "terrain.jpg").exists()), None)
    if directory is None:
        raise FileNotFoundError("missing prepared USGS terrain; run `uv run tools/data/prepare_terrain_relief.py`")
    payload = (directory / "terrain.bin").read_bytes()
    magic, version, rows, cols, width_m, depth_m, elevation_min, elevation_max = HEADER.unpack_from(payload)
    if magic[:7] != b"DVZTRN1" or version != 1 or rows < 2 or cols < 2:
        raise ValueError("invalid prepared terrain header")
    if HEADER.size + rows * cols * np.dtype("<f4").itemsize != len(payload):
        raise ValueError("unexpected prepared terrain size")
    heights = np.frombuffer(payload, "<f4", rows * cols, HEADER.size).astype(np.float64)
    heights -= elevation_min
    texture = np.asarray(Image.open(directory / "terrain.jpg").convert("RGBA"), dtype=np.uint8)
    return rows, cols, width_m, depth_m, heights, texture


def _build_scene():
    rows, cols, width_m, depth_m, heights, texture = _load()
    horizontal_scale = DISPLAY_WIDTH / width_m
    display_depth = depth_m * horizontal_scale
    desc = dvz.dvz_geometry_surface_grid_desc()
    desc.rows, desc.cols = rows, cols
    desc.heights = heights.ctypes.data_as(ctypes.POINTER(ctypes.c_double))
    desc.color = dvz.DvzColor(255, 255, 255, 255)
    desc.origin[:] = (-0.5 * DISPLAY_WIDTH, 0.0, +0.5 * display_depth)
    desc.col_basis[:] = (DISPLAY_WIDTH / (cols - 1), 0.0, 0.0)
    desc.row_basis[:] = (0.0, 0.0, -display_depth / (rows - 1))
    desc.height_axis[:] = (0.0, 1.0, 0.0)
    desc.height_scale = horizontal_scale * EXAGGERATION
    geometry = dvz.dvz_geometry_surface_grid(ctypes.byref(desc))
    if not geometry:
        raise RuntimeError("dvz_geometry_surface_grid() failed")

    scene, figure, panel = ex.scene_panel()
    dvz.dvz_panel_set_background_color(panel, dvz.DvzColor(6, 9, 10, 255))
    camera = dvz.dvz_camera_desc()
    camera.view.eye[:] = (+4.95, +3.95, +5.70)
    camera.view.target[:] = (0.0, +0.38, 0.0)
    camera.view.up[:] = (0.0, 1.0, 0.0)
    camera.projection.fov_y = 0.56
    camera.projection.near_clip = 0.03
    camera.projection.far_clip = 100.0
    dvz.dvz_panel_set_camera_desc(panel, ctypes.byref(camera))

    mesh = dvz.dvz_mesh(scene, 0)
    try:
        if dvz.dvz_mesh_set_geometry(mesh, geometry) != 0:
            raise RuntimeError("dvz_mesh_set_geometry() failed")
    finally:
        dvz.dvz_geometry_destroy(geometry)
    material = dvz.dvz_phong_material_desc()
    material.light_direction[:] = (-0.38, +0.76, +0.52)
    material.phong.ambient = 0.734
    material.phong.diffuse = 0.454
    material.phong.specular = 0.049
    material.phong.shininess = 7.594
    dvz.dvz_visual_set_material(mesh, ctypes.byref(material))
    field = dvz.dvz_sampled_field_from_array(scene, texture)
    if not field or dvz.dvz_visual_set_field(mesh, b"texture", field) != 0:
        raise RuntimeError("terrain texture setup failed")
    ex.add_visual(panel, mesh)
    return scene, figure, panel, rows, cols


def main() -> None:
    scene, figure, panel, rows, cols = _build_scene()
    print(f"terrain_relief: {cols} x {rows} prepared elevation grid")

    def configure(view) -> None:
        arcball = dvz.dvz_view_arcball(view, panel, None)
        if not arcball or dvz.dvz_arcball_set(arcball, INITIAL_ANGLES) != 0:
            raise RuntimeError("arcball setup failed")
        dvz.dvz_arcball_zoom(arcball, 0.606531)
        pan = (ctypes.c_float * 2)(-0.191645, +0.730556)
        dvz.dvz_arcball_pan(arcball, pan)

    ex.run_with_view(scene, figure, "McHenrys Peak Terrain Relief", configure)


if __name__ == "__main__":
    main()
/*
 * Copyright (c) 2021 Cyrille Rossant and contributors. All rights reserved.
 * Licensed under the MIT license. See LICENSE file in the project root for details.
 * SPDX-License-Identifier: MIT
 */

/* terrain_relief - This example drapes aligned NAIP orthoimagery over USGS 3DEP elevation.
 *
 * What to look for: the real bare-earth DEM drives a lit surface-grid mesh while the matching
 * natural-color aerial image is sampled through the grid UVs. Rotate the arcball view to inspect
 * McHenrys Peak, Glacier Gorge, alpine lakes, and the Continental Divide.
 *
 * Prepare the cache before running:
 *
 *   uv run tools/data/prepare_terrain_relief.py
 *
 * Scenario: showcases_terrain_relief
 * Style: showcase, graphite_cyan, 1280x720 window target
 *
 * Build:  just example-c showcases/terrain_relief
 * Run:    ./build/examples/c/showcases/terrain_relief --live
 * Smoke:  ./build/examples/c/showcases/terrain_relief --png
 */



/*************************************************************************************************/
/*  Includes                                                                                     */
/*************************************************************************************************/

#include <math.h>
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
#include <string.h>

#include "_alloc.h"
#include "_assertions.h"
#include "datoviz/fileio.h"
#include "datoviz/geom.h"
#include "datoviz/scene.h"
#include "example_common.h"
#include "example_style.h"
#include "example_tuner.h"
#include "runner/scenario_runner.h"



/*************************************************************************************************/
/*  Forward declarations                                                                         */
/*************************************************************************************************/

DvzScenarioSpec dvz_showcase_terrain_relief_scenario(void);



/*************************************************************************************************/
/*  Constants                                                                                    */
/*************************************************************************************************/

#define WIDTH  EXAMPLE_WINDOW_WIDTH
#define HEIGHT EXAMPLE_WINDOW_HEIGHT

#define TERRAIN_DATA_DIR  "data/examples/terrain_relief/prepared"
#define TERRAIN_CACHE_DIR ".cache/datoviz/examples/terrain_relief/prepared"
#define TERRAIN_BIN_NAME  "terrain.bin"
#define TERRAIN_JPEG_NAME "terrain.jpg"

#define TERRAIN_MAGIC       "DVZTRN1"
#define TERRAIN_VERSION     1u
#define TERRAIN_HEADER_SIZE 36u
#define TERRAIN_MAX_DIM     2048u

#define TERRAIN_DISPLAY_WIDTH 5.0
#define TERRAIN_EXAGGERATION  1.35
#define TERRAIN_BASE_HEIGHT   0.18f
#define TERRAIN_PREVIEW_TAU   6.28318530718f



/*************************************************************************************************/
/*  Structs                                                                                      */
/*************************************************************************************************/

typedef struct TerrainData
{
    uint32_t rows;
    uint32_t cols;
    float width_m;
    float depth_m;
    float elevation_min_m;
    float elevation_max_m;
    double* heights_m;
} TerrainData;


typedef struct TerrainReliefState
{
    ExampleTuner tuner;
    DvzArcball* arcball;
    DvzVisual* visual;
    DvzMaterialDesc material;
    DvzExampleGuiMsaaControls msaa;
} TerrainReliefState;



/*************************************************************************************************/
/*  Binary loading                                                                               */
/*************************************************************************************************/

/**
 * Read one little-endian uint32.
 *
 * @param bytes four input bytes
 * @return decoded value
 */
static uint32_t _u32_le(const uint8_t bytes[4])
{
    ANN(bytes);
    return (uint32_t)bytes[0] | ((uint32_t)bytes[1] << 8u) | ((uint32_t)bytes[2] << 16u) |
           ((uint32_t)bytes[3] << 24u);
}



/**
 * Read one little-endian float32.
 *
 * @param bytes four input bytes
 * @return decoded value
 */
static float _f32_le(const uint8_t bytes[4])
{
    const uint32_t value = _u32_le(bytes);
    float out = 0.0f;
    dvz_memcpy(&out, sizeof(out), &value, sizeof(value));
    return out;
}



/**
 * Join a prepared bundle directory and filename.
 *
 * @param directory bundle directory
 * @param filename prepared filename
 * @param out output path
 * @param out_size output buffer size
 * @return whether the path fit in the output buffer
 */
static bool _prepared_path(const char* directory, const char* filename, char* out, size_t out_size)
{
    ANN(directory);
    ANN(filename);
    ANN(out);
    if (out_size == 0)
        return false;
    const int count = snprintf(out, out_size, "%s/%s", directory, filename);
    return count > 0 && (size_t)count < out_size;
}



/**
 * Return whether a prepared file is readable without logging a failed candidate path.
 *
 * @param path file path
 * @return whether the file can be opened
 */
static bool _file_exists(const char* path)
{
    ANN(path);
    FILE* file = fopen(path, "rb");
    if (file == NULL)
        return false;
    fclose(file);
    return true;
}



/**
 * Resolve a complete prepared terrain bundle without mixing data and cache files.
 *
 * @param terrain_path output terrain binary path
 * @param texture_path output texture path
 * @param path_size size of both output buffers
 * @return whether a complete bundle was found
 */
static bool _resolve_bundle(char* terrain_path, char* texture_path, size_t path_size)
{
    ANN(terrain_path);
    ANN(texture_path);

    const char* directories[] = {TERRAIN_DATA_DIR, TERRAIN_CACHE_DIR};
    for (uint32_t i = 0; i < 2; i++)
    {
        if (!_prepared_path(directories[i], TERRAIN_BIN_NAME, terrain_path, path_size) ||
            !_prepared_path(directories[i], TERRAIN_JPEG_NAME, texture_path, path_size))
        {
            return false;
        }
        if (_file_exists(terrain_path) && _file_exists(texture_path))
            return true;
    }
    return false;
}



/**
 * Load and validate a prepared terrain binary.
 *
 * @param path terrain binary path
 * @param out decoded terrain data
 * @return whether the file was valid
 */
static bool _terrain_load(const char* path, TerrainData* out)
{
    ANN(path);
    ANN(out);

    DvzSize size = 0;
    uint8_t* bytes = (uint8_t*)dvz_read_file(path, &size);
    if (bytes == NULL || size < TERRAIN_HEADER_SIZE)
        goto error;

    if (memcmp(bytes, TERRAIN_MAGIC, strlen(TERRAIN_MAGIC)) != 0)
        goto error;
    const uint32_t version = _u32_le(&bytes[8]);
    out->rows = _u32_le(&bytes[12]);
    out->cols = _u32_le(&bytes[16]);
    out->width_m = _f32_le(&bytes[20]);
    out->depth_m = _f32_le(&bytes[24]);
    out->elevation_min_m = _f32_le(&bytes[28]);
    out->elevation_max_m = _f32_le(&bytes[32]);

    if (version != TERRAIN_VERSION || out->rows < 2 || out->cols < 2 ||
        out->rows > TERRAIN_MAX_DIM || out->cols > TERRAIN_MAX_DIM || !isfinite(out->width_m) ||
        !isfinite(out->depth_m) || out->width_m <= 0.0f || out->depth_m <= 0.0f ||
        !isfinite(out->elevation_min_m) || !isfinite(out->elevation_max_m) ||
        out->elevation_max_m <= out->elevation_min_m)
    {
        goto error;
    }

    const uint64_t count = (uint64_t)out->rows * out->cols;
    if (count > (UINT64_MAX - TERRAIN_HEADER_SIZE) / sizeof(float))
        goto error;
    const uint64_t expected_size = TERRAIN_HEADER_SIZE + count * sizeof(float);
    if (size != expected_size || count > SIZE_MAX / sizeof(double))
        goto error;

    out->heights_m = (double*)dvz_calloc((size_t)count, sizeof(double));
    if (out->heights_m == NULL)
        goto error;
    for (uint64_t i = 0; i < count; i++)
    {
        const float elevation = _f32_le(&bytes[TERRAIN_HEADER_SIZE + i * sizeof(float)]);
        if (!isfinite(elevation))
            goto error;
        out->heights_m[i] = (double)elevation - out->elevation_min_m;
    }

    dvz_free(bytes);
    return true;

error:
    dvz_free(bytes);
    dvz_free(out->heights_m);
    out->heights_m = NULL;
    return false;
}



/**
 * Release decoded terrain data.
 *
 * @param terrain terrain data
 */
static void _terrain_destroy(TerrainData* terrain)
{
    if (terrain == NULL)
        return;
    dvz_free(terrain->heights_m);
    terrain->heights_m = NULL;
}



/*************************************************************************************************/
/*  Scene construction                                                                           */
/*************************************************************************************************/

/**
 * Create the scene-owned sampled field for the prepared JPEG texture.
 *
 * @param scene owning scene
 * @param path texture path
 * @return sampled texture field, or NULL on failure
 */
static DvzSampledField* _create_texture(DvzScene* scene, const char* path)
{
    ANN(scene);
    ANN(path);

    uint32_t width = 0;
    uint32_t height = 0;
    uint8_t* pixels = dvz_read_jpeg(path, &width, &height);
    if (pixels == NULL || width == 0 || height == 0)
    {
        dvz_free(pixels);
        return NULL;
    }

    DvzSampledField* texture = dvz_sampled_field(
        scene, &(DvzSampledFieldDesc){
                   DVZ_STRUCT_INIT_FIELDS(DvzSampledFieldDesc), .dim = DVZ_FIELD_DIM_2D,
                   .format = DVZ_FIELD_FORMAT_RGBA8_UNORM, .semantic = DVZ_FIELD_SEMANTIC_COLOR,
                   .width = width, .height = height, .depth = 1});
    DvzResult result = DVZ_ERROR;
    if (texture != NULL)
    {
        result = dvz_sampled_field_set_data(
            texture, &(DvzFieldDataView){
                         DVZ_STRUCT_INIT_FIELDS(DvzFieldDataView), .data = pixels,
                         .bytes_per_row = width * 4u, .rows_per_image = height});
    }
    dvz_free(pixels);
    return result == DVZ_OK ? texture : NULL;
}



/**
 * Create the physically proportioned surface-grid geometry.
 *
 * @param terrain decoded terrain data
 * @param out_depth display depth output
 * @return surface geometry, or NULL on failure
 */
static DvzGeometry* _create_terrain_geometry(const TerrainData* terrain, float* out_depth)
{
    ANN(terrain);
    ANN(out_depth);

    const double horizontal_scale = TERRAIN_DISPLAY_WIDTH / terrain->width_m;
    const double display_depth = terrain->depth_m * horizontal_scale;
    *out_depth = (float)display_depth;

    DvzGeometrySurfaceGridDesc desc = dvz_geometry_surface_grid_desc();
    desc.rows = terrain->rows;
    desc.cols = terrain->cols;
    desc.heights = terrain->heights_m;
    desc.color = (DvzColor){255, 255, 255, 255};
    desc.origin[0] = -0.5 * TERRAIN_DISPLAY_WIDTH;
    desc.origin[2] = +0.5 * display_depth;
    desc.col_basis[0] = TERRAIN_DISPLAY_WIDTH / (double)(terrain->cols - 1u);
    desc.col_basis[1] = 0.0;
    desc.col_basis[2] = 0.0;
    desc.row_basis[0] = 0.0;
    desc.row_basis[1] = 0.0;
    desc.row_basis[2] = -display_depth / (double)(terrain->rows - 1u);
    desc.height_axis[0] = 0.0;
    desc.height_axis[1] = 1.0;
    desc.height_axis[2] = 0.0;
    desc.height_scale = horizontal_scale * TERRAIN_EXAGGERATION;
    return dvz_geometry_surface_grid(&desc);
}



/**
 * Append one outward-facing vertical terrain-skirt quad.
 *
 * @param skirt destination geometry
 * @param quad quad index
 * @param top0 first top edge position
 * @param top1 second top edge position
 */
static void _set_skirt_quad(DvzGeometry* skirt, uint32_t quad, const dvec3 top0, const dvec3 top1)
{
    ANN(skirt);
    ANN(top0);
    ANN(top1);

    const uint32_t vertex = 4u * quad;
    const uint32_t index = 6u * quad;
    const DvzColor color = {24, 32, 35, 255};
    for (uint32_t i = 0; i < 4; i++)
        skirt->colors[vertex + i] = color;

    for (uint32_t axis = 0; axis < 3; axis++)
    {
        skirt->positions[vertex + 0u][axis] = top0[axis];
        skirt->positions[vertex + 1u][axis] = top1[axis];
        skirt->positions[vertex + 2u][axis] = top1[axis];
        skirt->positions[vertex + 3u][axis] = top0[axis];
    }
    skirt->positions[vertex + 2u][1] = 0.0;
    skirt->positions[vertex + 3u][1] = 0.0;

    skirt->indices[index + 0u] = vertex + 0u;
    skirt->indices[index + 1u] = vertex + 1u;
    skirt->indices[index + 2u] = vertex + 2u;
    skirt->indices[index + 3u] = vertex + 0u;
    skirt->indices[index + 4u] = vertex + 2u;
    skirt->indices[index + 5u] = vertex + 3u;
}



/**
 * Create vertical walls around the structured terrain boundary.
 *
 * The four edges are traversed clockwise from above so their triangle winding faces outward.
 *
 * @param surface source surface-grid geometry
 * @return skirt geometry, or NULL on failure
 */
static DvzGeometry* _create_skirt_geometry(const DvzGeometry* surface)
{
    ANN(surface);
    if (surface->grid_rows < 2 || surface->grid_cols < 2 || surface->positions == NULL)
        return NULL;

    const uint32_t rows = surface->grid_rows;
    const uint32_t cols = surface->grid_cols;
    const uint32_t quad_count = 2u * (rows - 1u) + 2u * (cols - 1u);
    DvzGeometry* skirt = dvz_geometry(4u * quad_count, 6u * quad_count);
    if (skirt == NULL)
        return NULL;

    uint32_t quad = 0;
    for (uint32_t col = 0; col < cols - 1u; col++)
    {
        const uint32_t row = rows - 1u;
        _set_skirt_quad(
            skirt, quad++, surface->positions[row * cols + col],
            surface->positions[row * cols + col + 1u]);
    }
    for (uint32_t row = rows - 1u; row > 0; row--)
    {
        _set_skirt_quad(
            skirt, quad++, surface->positions[row * cols + cols - 1u],
            surface->positions[(row - 1u) * cols + cols - 1u]);
    }
    for (uint32_t col = cols - 1u; col > 0; col--)
    {
        _set_skirt_quad(skirt, quad++, surface->positions[col], surface->positions[col - 1u]);
    }
    for (uint32_t row = 0; row < rows - 1u; row++)
    {
        _set_skirt_quad(
            skirt, quad++, surface->positions[row * cols], surface->positions[(row + 1u) * cols]);
    }

    if (quad != quad_count || dvz_geometry_compute_normals(skirt) != DVZ_OK)
    {
        dvz_geometry_destroy(skirt);
        return NULL;
    }
    return skirt;
}



/**
 * Add the vertical relief-map skirt around the terrain.
 *
 * @param scene owning scene
 * @param panel target panel
 * @param surface source surface geometry
 * @return whether the skirt was added
 */
static bool _add_skirt(DvzScene* scene, DvzPanel* panel, const DvzGeometry* surface)
{
    ANN(scene);
    ANN(panel);
    ANN(surface);

    DvzGeometry* skirt = _create_skirt_geometry(surface);
    if (skirt == NULL)
        return false;
    DvzVisual* visual = dvz_mesh(scene, 0);
    DvzResult result = visual != NULL ? dvz_mesh_set_geometry(visual, skirt) : DVZ_ERROR;
    dvz_geometry_destroy(skirt);
    if (result != DVZ_OK)
        return false;

    DvzMaterialDesc material = dvz_phong_material_desc();
    material.light_direction[0] = -0.42f;
    material.light_direction[1] = +0.72f;
    material.light_direction[2] = +0.55f;
    material.phong.ambient = 0.36f;
    material.phong.diffuse = 0.58f;
    material.phong.specular = 0.04f;
    material.phong.shininess = 16.0f;
    result = dvz_visual_set_material(visual, &material);
    if (result == DVZ_OK)
        result = dvz_panel_add_visual(panel, visual, NULL);
    return result == DVZ_OK;
}



/**
 * Add the dark relief-map base beneath the terrain.
 *
 * @param scene owning scene
 * @param panel target panel
 * @param display_depth terrain depth in display units
 * @return whether the base was added
 */
static bool _add_base(DvzScene* scene, DvzPanel* panel, float display_depth)
{
    ANN(scene);
    ANN(panel);

    DvzGeometryCubeDesc cube_desc = dvz_geometry_cube_desc();
    cube_desc.color = (DvzColor){20, 28, 32, 255};
    DvzGeometry* cube = dvz_geometry_cube(&cube_desc);
    if (cube == NULL)
        return false;

    DvzVisual* base = dvz_mesh(scene, 0);
    DvzResult result = base != NULL ? dvz_mesh_set_geometry(base, cube) : DVZ_ERROR;
    dvz_geometry_destroy(cube);
    if (result != DVZ_OK)
        return false;

    DvzMaterialDesc material = dvz_phong_material_desc();
    material.light_direction[0] = -0.42f;
    material.light_direction[1] = +0.72f;
    material.light_direction[2] = +0.55f;
    material.phong.ambient = 0.32f;
    material.phong.diffuse = 0.55f;
    material.phong.specular = 0.08f;
    material.phong.shininess = 18.0f;
    result = dvz_visual_set_material(base, &material);

    mat4 transform = {
        {TERRAIN_DISPLAY_WIDTH + 0.12f, 0.0f, 0.0f, 0.0f},
        {0.0f, TERRAIN_BASE_HEIGHT, 0.0f, 0.0f},
        {0.0f, 0.0f, display_depth + 0.12f, 0.0f},
        {0.0f, -0.5f * TERRAIN_BASE_HEIGHT, 0.0f, 1.0f},
    };
    if (result == DVZ_OK)
        result = dvz_visual_set_transform(base, transform);
    if (result == DVZ_OK)
        result = dvz_panel_add_visual(panel, base, NULL);
    return result == DVZ_OK;
}



/**
 * Add the textured terrain mesh.
 *
 * @param scene owning scene
 * @param panel target panel
 * @param geometry terrain geometry
 * @param texture aligned texture field
 * @return terrain visual, or NULL on failure
 */
static DvzVisual* _add_terrain(
    DvzScene* scene, DvzPanel* panel, const DvzGeometry* geometry, DvzSampledField* texture,
    DvzMaterialDesc* material)
{
    ANN(scene);
    ANN(panel);
    ANN(geometry);
    ANN(texture);
    ANN(material);

    DvzVisual* visual = dvz_mesh(scene, 0);
    if (visual == NULL)
        return NULL;
    DvzResult result = dvz_mesh_set_geometry(visual, geometry);

    if (result == DVZ_OK)
        result = dvz_visual_set_material(visual, material);
    if (result == DVZ_OK)
        result = dvz_visual_set_field(visual, "texture", texture);
    if (result == DVZ_OK)
        result = dvz_panel_add_visual(panel, visual, NULL);
    return result == DVZ_OK ? visual : NULL;
}



/*************************************************************************************************/
/*  Scenario callbacks                                                                           */
/*************************************************************************************************/

/**
 * Initialize the real-terrain showcase.
 *
 * @param ctx scenario context
 * @param out_user unused scenario state output
 * @return true on success
 */
static bool _scenario_init(DvzScenarioContext* ctx, void** out_user)
{
    if (ctx == NULL)
        return false;
    TerrainReliefState* state = (TerrainReliefState*)dvz_calloc(1, sizeof(*state));
    if (state == NULL)
        return false;
    state->tuner = example_tuner("Terrain relief settings");
    if (out_user != NULL)
        *out_user = state;

    bool ok = false;
    char terrain_path[1024] = {0};
    char texture_path[1024] = {0};
    TerrainData terrain = {0};
    DvzGeometry* geometry = NULL;

    if (!_resolve_bundle(terrain_path, texture_path, sizeof(terrain_path)))
    {
        dvz_fprintf(
            stderr, "terrain_relief: missing prepared USGS terrain. Run "
                    "`uv run tools/data/prepare_terrain_relief.py` from the repository root.\n");
        goto cleanup;
    }
    EXAMPLE_CHECK(_terrain_load(terrain_path, &terrain), "invalid prepared terrain binary");

    DvzSampledField* texture = _create_texture(ctx->scene, texture_path);
    EXAMPLE_CHECK(texture != NULL, "unable to load prepared terrain texture");

    float display_depth = 0.0f;
    geometry = _create_terrain_geometry(&terrain, &display_depth);
    EXAMPLE_CHECK(geometry != NULL, "dvz_geometry_surface_grid() failed");

    ctx->figure = dvz_figure(ctx->scene, ctx->width, ctx->height, 0);
    EXAMPLE_CHECK(ctx->figure != NULL, "dvz_figure() failed");
    example_tuner_figure(&state->tuner, ctx->figure);
    DvzPanel* panel = dvz_panel_full(ctx->figure);
    EXAMPLE_CHECK(panel != NULL, "dvz_panel_full() failed");
    example_graphite_cyan_set_panel_background(panel);
    dvz_panel_set_background_color(panel, dvz_color_from_unit(0.025f, 0.034f, 0.040f, 1.0f));

    DvzCameraDesc camera = dvz_camera_desc();
    camera.projection.type = DVZ_CAMERA_PERSPECTIVE;
    camera.view.eye[0] = +4.95f;
    camera.view.eye[1] = +3.95f;
    camera.view.eye[2] = +5.70f;
    camera.view.target[0] = 0.0f;
    camera.view.target[1] = +0.38f;
    camera.view.target[2] = 0.0f;
    camera.view.up[0] = 0.0f;
    camera.view.up[1] = 1.0f;
    camera.view.up[2] = 0.0f;
    camera.projection.fov_y = 0.56f;
    camera.projection.near_clip = 0.03f;
    camera.projection.far_clip = 100.0f;
    camera.projection.ortho_height = 0.0f;
    EXAMPLE_CHECK(
        dvz_panel_set_camera_desc(panel, &camera) == DVZ_OK, "dvz_panel_set_camera_desc() failed");

    DvzCamera* camera_ref = dvz_panel_camera(panel);
    EXAMPLE_CHECK(camera_ref != NULL, "dvz_panel_camera() failed");

    EXAMPLE_CHECK(_add_base(ctx->scene, panel, display_depth), "failed to add terrain base");
    EXAMPLE_CHECK(_add_skirt(ctx->scene, panel, geometry), "failed to add terrain skirt");
    state->material = dvz_material_desc();
    state->material.model = DVZ_MATERIAL_MODEL_PHONG;
    state->material.alpha_mode = DVZ_ALPHA_OPAQUE;
    state->material.opacity = 1.0f;
    state->material.base_color_factor[0] = 1.0f;
    state->material.base_color_factor[1] = 1.0f;
    state->material.base_color_factor[2] = 1.0f;
    state->material.base_color_factor[3] = 1.0f;
    state->material.light_direction[0] = -0.38f;
    state->material.light_direction[1] = +0.76f;
    state->material.light_direction[2] = +0.52f;
    state->material.phong.ambient = 0.734f;
    state->material.phong.diffuse = 0.454f;
    state->material.phong.specular = 0.049f;
    state->material.phong.shininess = 7.594f;
    state->material.standard.roughness = 0.62f;
    state->material.standard.specular = 0.34f;
    state->material.standard.metallic = 0.0f;
    state->material.standard.emissive[0] = 0.0f;
    state->material.standard.emissive[1] = 0.0f;
    state->material.standard.emissive[2] = 0.0f;
    state->material.standard.rim_strength = 0.10f;
    state->visual = _add_terrain(ctx->scene, panel, geometry, texture, &state->material);
    EXAMPLE_CHECK(state->visual != NULL, "failed to add textured terrain");
    EXAMPLE_CHECK(
        dvz_scenario_set_primary_visual(ctx, state->visual) == DVZ_OK,
        "dvz_scenario_set_primary_visual() failed");

    state->msaa = (DvzExampleGuiMsaaControls){
        .enabled = true,
        .alpha_to_coverage = false,
        .samples = 8.0f,
        .min_samples = 2.0f,
        .max_samples = 16.0f,
    };
    DvzMsaaDesc msaa_desc = dvz_msaa_desc();
    msaa_desc.enabled = true;
    msaa_desc.sample_count = 8u;
    msaa_desc.alpha_to_coverage = false;
    EXAMPLE_CHECK(dvz_panel_set_msaa(panel, &msaa_desc) == DVZ_OK, "dvz_panel_set_msaa() failed");

    DvzController* controller = dvz_arcball(ctx->scene, NULL);
    EXAMPLE_CHECK(controller != NULL, "dvz_arcball() failed");
    state->arcball = dvz_controller_arcball(controller);
    EXAMPLE_CHECK(state->arcball != NULL, "dvz_controller_arcball() failed");
    EXAMPLE_CHECK(
        dvz_scenario_bind_controller(ctx, panel, controller, DVZ_DIM_MASK_XYZ) == DVZ_OK,
        "dvz_scenario_bind_controller() failed");

    vec3 arcball_angles = {+0.050954f, -0.163512f, -0.032719f};
    vec2 arcball_pan = {-0.191645f, +0.730556f};
    EXAMPLE_CHECK(
        dvz_arcball_initial(state->arcball, arcball_angles) == DVZ_OK,
        "dvz_arcball_initial() failed");
    EXAMPLE_CHECK(
        dvz_arcball_zoom(state->arcball, 0.606531f) == DVZ_OK, "dvz_arcball_zoom() failed");
    EXAMPLE_CHECK(
        dvz_arcball_pan(state->arcball, arcball_pan) == DVZ_OK, "dvz_arcball_pan() failed");
    example_tuner_camera_ref(&state->tuner, "Camera", panel, camera_ref, &camera);
    example_tuner_arcball(
        &state->tuner, "Arcball", state->arcball, arcball_angles, 0.606531f, arcball_pan);
    example_tuner_material(&state->tuner, "Terrain material", state->visual, &state->material);
    example_tuner_msaa(&state->tuner, "MSAA", panel, &state->msaa);
    ok = true;

cleanup:
    dvz_geometry_destroy(geometry);
    _terrain_destroy(&terrain);
    return ok;
}



/**
 * Attach the live terrain tuner after the native view exists.
 *
 * @param ctx scenario context
 * @param app native app
 * @param view native view
 * @param user scenario state
 * @return whether the tuner was attached or intentionally skipped
 */
static bool _scenario_native_view(DvzScenarioContext* ctx, DvzApp* app, DvzView* view, void* user)
{
    (void)app;
    TerrainReliefState* state = (TerrainReliefState*)user;
    if (ctx == NULL || ctx->presentation != DVZ_RUNNER_PRESENT_GLFW || state == NULL ||
        view == NULL)
    {
        return true;
    }
    return example_tuner_attach(&state->tuner, view);
}



/**
 * Rotate the complete terrain scene through one seamless preview orbit.
 *
 * @param ctx scenario context
 * @param user scenario state
 */
static void _scenario_frame(DvzScenarioContext* ctx, void* user)
{
    TerrainReliefState* state = (TerrainReliefState*)user;
    if (ctx == NULL || !ctx->preview_mode || state == NULL || state->arcball == NULL)
        return;

    const float phase =
        (float)dvz_scenario_preview_phase(ctx, DVZ_SCENARIO_PREVIEW_PHASE_SEAMLESS_LOOP);
    vec3 angles = {+0.050954f, -0.163512f + TERRAIN_PREVIEW_TAU * phase, -0.032719f};
    (void)dvz_arcball_set(state->arcball, angles);
}



/**
 * Destroy the terrain-relief showcase state.
 *
 * @param ctx scenario context
 * @param user scenario state
 */
static void _scenario_destroy(DvzScenarioContext* ctx, void* user)
{
    (void)ctx;
    TerrainReliefState* state = (TerrainReliefState*)user;
    if (state == NULL)
        return;
    example_tuner_detach(&state->tuner);
    dvz_free(state);
}



/**
 * Return the terrain-relief showcase scenario specification.
 *
 * @return scenario specification
 */
DvzScenarioSpec dvz_showcase_terrain_relief_scenario(void)
{
    return (DvzScenarioSpec){
        .id = "showcases_terrain_relief",
        .title = "McHenrys Peak Terrain Relief",
        .width = WIDTH,
        .height = HEIGHT,
        .fps = 60.0,
        .requirements =
            DVZ_SCENARIO_REQ_MESH_VISUAL | DVZ_SCENARIO_REQ_CONTROLLER | DVZ_SCENARIO_REQ_ARCBALL,
        .init = _scenario_init,
        .frame = _scenario_frame,
        .destroy = _scenario_destroy,
    };
}



/*************************************************************************************************/
/*  Functions                                                                                    */
/*************************************************************************************************/

/**
 * Run the terrain-relief showcase through the native scenario runner.
 *
 * @param argc command-line argument count
 * @param argv command-line argument vector
 * @return process exit code
 */
#ifndef DVZ_EXAMPLE_NO_MAIN
int main(int argc, char** argv)
{
    DvzScenarioSpec spec = dvz_showcase_terrain_relief_scenario();
    if (example_cli_wants_live_gui(argc, argv))
        spec.native_view = _scenario_native_view;
    return dvz_scenario_run_native_cli(&spec, argc, argv) == 0 ? 0 : 1;
}
#endif
#!/usr/bin/env -S uv run --script
# /// script
# requires-python = ">=3.10"
# dependencies = [
#   "numpy==2.2.6",
#   "pillow==11.2.1",
#   "pyproj==3.7.1",
#   "requests==2.32.4",
# ]
# ///
"""Prepare aligned USGS 3DEP elevation and NAIP imagery for the terrain showcase."""

from __future__ import annotations

import argparse
import json
import shutil
import struct
from pathlib import Path
from typing import Any

import numpy as np
import requests
from PIL import Image, ImageEnhance, ImageOps
from pyproj import Transformer

from common import (
    CACHE_ROOT,
    REPO_ROOT,
    artifact,
    command_argv,
    write_json,
    write_manifest,
    write_provenance,
)


EXAMPLE_ID = "terrain_relief"
DEFAULT_OUTPUT = CACHE_ROOT / EXAMPLE_ID

ELEVATION_SERVICE = (
    "https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer"
)
IMAGERY_SERVICE = (
    "https://imagery.nationalmap.gov/arcgis/rest/services/USGSNAIPPlus/ImageServer"
)

# McHenrys Peak, Glacier Gorge, and the Continental Divide in Rocky Mountain National Park.
WGS84_BBOX = (-105.705, 40.245, -105.625, 40.310)
OUTPUT_CRS = 26913  # NAD83 / UTM zone 13N.
ELEVATION_ROWS = 512
TEXTURE_ROWS = 2048

MAGIC = b"DVZTRN1\0"
VERSION = 1
USER_AGENT = "Datoviz-v0.4-terrain-relief-example/1.0"


def _service_metadata(session: requests.Session, service_url: str) -> dict[str, Any]:
    """Fetch the public ArcGIS image-service metadata."""
    response = session.get(service_url, params={"f": "json"}, timeout=60)
    response.raise_for_status()
    payload = response.json()
    if "error" in payload:
        raise RuntimeError(f"ArcGIS service metadata error: {payload['error']}")
    return payload


def _download_file(session: requests.Session, url: str, path: Path) -> None:
    """Download one response body atomically."""
    path.parent.mkdir(parents=True, exist_ok=True)
    tmp = path.with_suffix(path.suffix + ".tmp")
    with session.get(url, stream=True, timeout=180) as response:
        response.raise_for_status()
        with tmp.open("wb") as file:
            for chunk in response.iter_content(chunk_size=1024 * 1024):
                if chunk:
                    file.write(chunk)
    tmp.replace(path)


def _export_image(
    session: requests.Session,
    service_url: str,
    path: Path,
    metadata_path: Path,
    *,
    cols: int,
    rows: int,
    image_format: str,
    pixel_type: str | None = None,
    refresh: bool,
) -> dict[str, Any]:
    """Export and cache one projected ArcGIS image-service raster."""
    if path.exists() and metadata_path.exists() and not refresh:
        return json.loads(metadata_path.read_text(encoding="utf8"))["export"]

    params: dict[str, str] = {
        "bbox": ",".join(f"{value:.9f}" for value in WGS84_BBOX),
        "bboxSR": "4326",
        "imageSR": str(OUTPUT_CRS),
        "size": f"{cols},{rows}",
        "format": image_format,
        "interpolation": "RSP_BilinearInterpolation",
        "f": "json",
    }
    if pixel_type is not None:
        params["pixelType"] = pixel_type

    endpoint = f"{service_url}/exportImage"
    response = session.get(endpoint, params=params, timeout=180)
    response.raise_for_status()
    export = response.json()
    if "error" in export:
        raise RuntimeError(f"ArcGIS export error: {export['error']}")
    href = export.get("href")
    if not isinstance(href, str) or not href.startswith("https://"):
        raise RuntimeError(f"ArcGIS export did not return an HTTPS download: {export}")

    _download_file(session, href, path)
    write_json(
        metadata_path,
        {
            "endpoint": endpoint,
            "parameters": params,
            "export": export,
        },
    )
    return export


def _projected_dimensions() -> tuple[float, float]:
    """Return the requested WGS84 crop dimensions in UTM meters."""
    west, south, east, north = WGS84_BBOX
    transformer = Transformer.from_crs(4326, OUTPUT_CRS, always_xy=True)
    xmin, ymin = transformer.transform(west, south)
    xmax, ymax = transformer.transform(east, north)
    return float(xmax - xmin), float(ymax - ymin)


def _cols_for_rows(rows: int, width_m: float, depth_m: float) -> int:
    """Choose a column count that keeps projected pixels approximately square."""
    return max(2, int(round(rows * width_m / depth_m)))


def _extent(export: dict[str, Any]) -> tuple[float, float, float, float]:
    """Read one projected extent from an ArcGIS export response."""
    extent = export.get("extent")
    if not isinstance(extent, dict):
        raise RuntimeError("ArcGIS export response is missing its projected extent")
    try:
        keys = ("xmin", "ymin", "xmax", "ymax")
        return tuple(float(extent[key]) for key in keys)  # type: ignore[return-value]
    except (KeyError, TypeError, ValueError) as exc:
        raise RuntimeError(f"invalid ArcGIS export extent: {extent}") from exc


def _assert_aligned(elevation_export: dict[str, Any], imagery_export: dict[str, Any]) -> None:
    """Verify that elevation and imagery cover the same projected rectangle."""
    elevation_extent = np.asarray(_extent(elevation_export), dtype=np.float64)
    imagery_extent = np.asarray(_extent(imagery_export), dtype=np.float64)
    if not np.allclose(elevation_extent, imagery_extent, rtol=0.0, atol=0.02):
        raise RuntimeError(
            "USGS elevation and imagery exports are not aligned: "
            f"{elevation_extent.tolist()} != {imagery_extent.tolist()}"
        )


def _load_elevation(path: Path, rows: int, cols: int) -> np.ndarray:
    """Load and validate the exported float32 elevation grid."""
    with Image.open(path) as image:
        elevation = np.asarray(image, dtype=np.float32)
    if elevation.shape != (rows, cols):
        raise RuntimeError(f"unexpected elevation shape {elevation.shape}; expected {(rows, cols)}")
    if not np.all(np.isfinite(elevation)):
        raise RuntimeError("elevation export contains non-finite samples")
    return np.ascontiguousarray(elevation, dtype="<f4")


def _prepare_texture(source: Path, output: Path, rows: int, cols: int) -> None:
    """Apply a restrained gallery grade while preserving the aligned NAIP pixels."""
    with Image.open(source) as image:
        texture = ImageOps.exif_transpose(image).convert("RGB")
        if texture.size != (cols, rows):
            raise RuntimeError(
                f"unexpected imagery size {texture.size}; expected {(cols, rows)}"
            )
        texture = ImageOps.autocontrast(texture, cutoff=(0.25, 0.25))
        texture = ImageEnhance.Color(texture).enhance(1.08)
        texture = ImageEnhance.Contrast(texture).enhance(1.06)
        texture = ImageEnhance.Brightness(texture).enhance(0.96)
        output.parent.mkdir(parents=True, exist_ok=True)
        texture.save(output, format="JPEG", quality=92, subsampling=0, optimize=True)


def _write_elevation(
    path: Path,
    elevation: np.ndarray,
    *,
    width_m: float,
    depth_m: float,
) -> None:
    """Write the little-endian terrain-grid binary consumed by the C example."""
    rows, cols = elevation.shape
    header = struct.pack(
        "<8sIIIffff",
        MAGIC,
        VERSION,
        rows,
        cols,
        width_m,
        depth_m,
        float(np.min(elevation)),
        float(np.max(elevation)),
    )
    path.parent.mkdir(parents=True, exist_ok=True)
    with path.open("wb") as file:
        file.write(header)
        file.write(elevation.tobytes(order="C"))


def prepare(args: argparse.Namespace) -> None:
    """Prepare the cache-only McHenrys Peak terrain bundle."""
    bundle_root = args.output.resolve()
    source = bundle_root / "source"
    prepared = bundle_root / "prepared"
    if args.refresh_source and source.exists():
        shutil.rmtree(source)
    if args.force and prepared.exists():
        shutil.rmtree(prepared)
    source.mkdir(parents=True, exist_ok=True)
    prepared.mkdir(parents=True, exist_ok=True)

    approximate_width_m, approximate_depth_m = _projected_dimensions()
    elevation_cols = _cols_for_rows(
        ELEVATION_ROWS, approximate_width_m, approximate_depth_m
    )
    texture_cols = _cols_for_rows(TEXTURE_ROWS, approximate_width_m, approximate_depth_m)

    session = requests.Session()
    session.headers.update({"User-Agent": USER_AGENT})
    elevation_service = _service_metadata(session, ELEVATION_SERVICE)
    imagery_service = _service_metadata(session, IMAGERY_SERVICE)

    elevation_source = source / "mchenrys_peak_3dep.tif"
    imagery_source = source / "mchenrys_peak_naip.jpg"
    elevation_export = _export_image(
        session,
        ELEVATION_SERVICE,
        elevation_source,
        source / "mchenrys_peak_3dep.export.json",
        cols=elevation_cols,
        rows=ELEVATION_ROWS,
        image_format="tiff",
        pixel_type="F32",
        refresh=args.refresh_source,
    )
    imagery_export = _export_image(
        session,
        IMAGERY_SERVICE,
        imagery_source,
        source / "mchenrys_peak_naip.export.json",
        cols=texture_cols,
        rows=TEXTURE_ROWS,
        image_format="jpgpng",
        refresh=args.refresh_source,
    )
    _assert_aligned(elevation_export, imagery_export)

    xmin, ymin, xmax, ymax = _extent(elevation_export)
    width_m = xmax - xmin
    depth_m = ymax - ymin
    elevation = _load_elevation(elevation_source, ELEVATION_ROWS, elevation_cols)

    elevation_output = prepared / "terrain.bin"
    texture_output = prepared / "terrain.jpg"
    _write_elevation(elevation_output, elevation, width_m=width_m, depth_m=depth_m)
    _prepare_texture(imagery_source, texture_output, TEXTURE_ROWS, texture_cols)

    source_metadata = {
        "name": "USGS McHenrys Peak terrain and orthoimagery",
        "bbox_wgs84": list(WGS84_BBOX),
        "output_crs": f"EPSG:{OUTPUT_CRS}",
        "elevation": {
            "service": ELEVATION_SERVICE,
            "description": elevation_service.get("serviceDescription", ""),
            "copyright": elevation_service.get("copyrightText", ""),
        },
        "imagery": {
            "service": IMAGERY_SERVICE,
            "description": imagery_service.get("serviceDescription", ""),
            "copyright": imagery_service.get("copyrightText", ""),
        },
    }
    command = command_argv(
        "tools/data/prepare_terrain_relief.py",
        ["--output", str(bundle_root.relative_to(REPO_ROOT))]
        if bundle_root.is_relative_to(REPO_ROOT)
        else ["--output", str(bundle_root)],
    )
    write_manifest(
        bundle_root,
        example_id=EXAMPLE_ID,
        title="McHenrys Peak Terrain Relief",
        status="cache-only",
        script="tools/data/prepare_terrain_relief.py",
        command=command,
        source=source_metadata,
        artifacts=[
            artifact(
                elevation_output,
                bundle_root,
                "elevation-grid",
                "DVZTRN1 little-endian float32",
                rows=ELEVATION_ROWS,
                cols=elevation_cols,
                elevation_min_m=float(np.min(elevation)),
                elevation_max_m=float(np.max(elevation)),
            ),
            artifact(
                texture_output,
                bundle_root,
                "surface-texture",
                "JPEG RGB8",
                width=texture_cols,
                height=TEXTURE_ROWS,
            ),
        ],
        validation={
            "aligned_projected_extent": [xmin, ymin, xmax, ymax],
            "width_m": width_m,
            "depth_m": depth_m,
            "finite_elevation": True,
        },
        extra={
            "encoding": {
                "terrain_header": "<8sIIIffff",
                "row_order": "north-to-south",
                "column_order": "west-to-east",
                "elevation_units": "meters NAVD88",
                "horizontal_crs": f"EPSG:{OUTPUT_CRS}",
            }
        },
    )
    write_provenance(
        bundle_root,
        title="McHenrys Peak Terrain Relief",
        source_lines=[
            "Bare-earth elevation: USGS 3D Elevation Program dynamic image service.",
            "Natural-color texture: USGS The National Map NAIP Plus image service.",
            f"WGS84 crop: {WGS84_BBOX} around McHenrys Peak and Glacier Gorge, Colorado.",
        ],
        processing_lines=[
            f"Both rasters were exported into EPSG:{OUTPUT_CRS} over the same projected extent.",
            f"Elevation was sampled to {elevation_cols}x{ELEVATION_ROWS} float32 meters.",
            f"Imagery was sampled to {texture_cols}x{TEXTURE_ROWS} RGB pixels.",
            "The imagery received restrained autocontrast, saturation, contrast, and "
            "brightness grading.",
        ],
        license_lines=[
            "USGS 3DEP products are public domain and available without use restrictions.",
            "NAIP imagery acquired by USDA has been placed in the public domain.",
            "Credit: U.S. Geological Survey 3DEP and USDA National Agriculture Imagery Program.",
        ],
        notes=[
            "The ArcGIS service descriptions and exact export requests are retained under source/.",
            "Service-backed source mosaics evolve; the manifest records checksums for this "
            "prepared snapshot.",
        ],
    )

    print(f"prepared terrain relief in {bundle_root}")
    print(
        f"grid {elevation_cols}x{ELEVATION_ROWS}, texture {texture_cols}x{TEXTURE_ROWS}, "
        f"elevation {float(np.min(elevation)):.1f}-{float(np.max(elevation)):.1f} m"
    )


def main() -> int:
    """Parse command-line arguments and prepare the terrain cache."""
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--output", type=Path, default=DEFAULT_OUTPUT)
    parser.add_argument("--force", action="store_true", help="replace prepared outputs")
    parser.add_argument(
        "--refresh-source", action="store_true", help="redownload both USGS service exports"
    )
    args = parser.parse_args()
    prepare(args)
    return 0


if __name__ == "__main__":
    raise SystemExit(main())
Example details

Tags

scientific, real-data, terrain, dem, orthoimagery, surface-grid, mesh, sampled-field, texture, material, arcball, camera, msaa

Data

Field Value
kind prepared

Dataset

Field Value
name McHenrys Peak USGS 3DEP terrain and NAIP orthoimagery
source USGS 3D Elevation Program and USDA National Agriculture Imagery Program
source_url https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer; https://imagery.nationalmap.gov/arcgis/rest/services/USGSNAIPPlus/ImageServer
license USGS 3DEP products and USDA NAIP imagery are public domain; credit the U.S. Geological Survey and USDA National Agriculture Imagery Program as the data sources.
citation U.S. Geological Survey 3D Elevation Program; USDA National Agriculture Imagery Program.
preprocessing uv run tools/data/prepare_terrain_relief.py
prepared_source data/examples/terrain_relief/prepared/terrain.bin; data/examples/terrain_relief/prepared/terrain.jpg
cache_prepared_path .cache/datoviz/examples/terrain_relief/prepared
provenance An aligned 6.95 by 7.27 km crop around McHenrys Peak and Glacier Gorge is exported from the public USGS image services into NAD83 / UTM zone 13N. The bare-earth DEM drives the mesh and the natural-color orthoimage supplies its sampled texture. Exact export requests, service descriptions, extents, and artifact checksums are retained in the cache manifest.

Encoding

Field Value
elevation 490 by 512 little-endian float32 meters, row-major north to south
texture 1960 by 2048 RGB JPEG over the identical projected extent
geometry physical UTM aspect ratio with 1.35x vertical exaggeration