sudo apt-get install gdal-bin \
  libgdal-dev
  python3-gdal

---

DEM terrain from OpenTopography

This page shows how to download a real Digital Elevation Model (DEM) from OpenTopography and use it as terrain in Gazebo Harmonic.

The default example downloads a 2 km x 2 km area using USGS10m, if the selected area is covered by USGS 3DEP. For a lighter Gazebo world, use the 0.5 km x 0.5 km option shown below. A 10 m DEM gives about 50 x 50 height samples for a 0.5 km square.

Note

USGS10m is a USGS 3DEP dataset, so it covers the United States, Alaska, Hawaii, and U.S. territories. For other countries, use a global dataset such as COP30, NASADEM, AW3D30, or SRTMGL1. Those are usually 30 m resolution, not 10 m.

request size

A 2 km x 2 km world is only 4 km2. OpenTopography documents a much larger request limit for USGS10m, so the simulation size here is chosen for Gazebo performance, not because the API requires such a small area.

Install tools

Gazebo uses GDAL internally for DEM files. Use GDAL on the command line to inspect and reproject the DEM before loading it in simulation.

sudo apt update
sudo apt install gdal-bin python3-gdal python3-numpy curl

Check:

gdalinfo --version
python3 -c "from osgeo import gdal; print(gdal.VersionInfo())"
gz sim --versions

Get an OpenTopography API key

1. Open OpenTopography.
2. Create an account or sign in.
3. Open your MyOpenTopo dashboard.
4. Click Get an API Key or Request API key.
5. Store it in an environment variable:

export OPENTOPOGRAPHY_API_KEY="your_key_here"

Do not commit the API key into the repository.

OpenTopography also documents daily API limits. Keep generated examples and notebooks using OPENTOPOGRAPHY_API_KEY from the environment instead of hard-coding a key.

Pick an area

Choose the center of the terrain as latitude and longitude.

Example near Mount Baker, Washington:

export CENTER_LAT=48.7769
export CENTER_LON=-121.8144

The download script converts the center point into an approximate square latitude/longitude bounding box. By default it downloads 2000 m x 2000 m.

Download script

./download_usgs10m_2km.sh "$CENTER_LAT" "$CENTER_LON"

Output:

raw/ot_usgs10m_2000m.tif

For a smaller 0.5 km x 0.5 km area, pass 500 as the third argument:

./download_usgs10m_2km.sh "$CENTER_LAT" "$CENTER_LON" 500

Output:

raw/ot_usgs10m_500m.tif

The script calls:

https://portal.opentopography.org/API/usgsdem

with:

Parameter	Value
datasetName	USGS10m
south north west east	bounding box around the center point
outputFormat	GTiff
API_Key	from OPENTOPOGRAPHY_API_KEY

Prepare the DEM for Gazebo

The downloaded DEM is usually in geographic coordinates, meaning latitude and longitude degrees. For simulation, reproject it to a meter-based coordinate system.

The script below chooses the UTM zone from the longitude and resamples to 10 m cells:

Prepare script

./prepare_dem_utm.sh "$CENTER_LAT" "$CENTER_LON"

Output:

dem_2000m_10m_utm.tif
dem_2000m_10m_utm.info.txt
textures/dem_diffuse.png
textures/flat_normal.png

For the 0.5 km x 0.5 km DEM, pass the same size:

./prepare_dem_utm.sh "$CENTER_LAT" "$CENTER_LON" 500

Output:

dem_500m_10m_utm.tif
dem_500m_10m_utm.info.txt

Inspect the file:

gdalinfo -stats dem_500m_10m_utm.tif

Look for:

• Size is ...
• Pixel Size = (10.000..., -10.000...)
• Minimum=...
• Maximum=...

For a 2 km x 2 km area at 10 m resolution, expect roughly 200 x 200 pixels. For a 0.5 km x 0.5 km area at 10 m resolution, expect roughly 50 x 50 pixels. It may not be exact because the first crop is defined in latitude/longitude.

Generate the Gazebo world

Generate the SDF from the actual DEM metadata:

World generator

./make_dem_world.py dem_500m_10m_utm.tif dem_world.sdf

Then run:

gz sim -r -v4 dem_world.sdf

The generated SDF uses the DEM as both collision and visual geometry:

<heightmap>
  <uri>dem_500m_10m_utm.tif</uri>
  <size>500 500 120</size>
  <pos>0 0 -650</pos>
  <sampling>1</sampling>
  <texture>
    <diffuse>textures/dem_diffuse.png</diffuse>
    <normal>textures/flat_normal.png</normal>
    <size>50</size>
  </texture>
</heightmap>

The exact size and pos values are generated from gdalinfo:

• x = raster_width_pixels * pixel_width_m
• y = raster_height_pixels * pixel_height_m
• z = Maximum elevation - Minimum elevation
• pos z = -Maximum elevation

Using pos z = -Maximum elevation moves the highest point near z = 0, which makes the terrain easier to view from the default Gazebo camera. If you want the lowest point near z = 0, change pos z to -Minimum elevation.

Fuel-style mesh model

The Gazebo moon examples show two useful patterns:

• dem_moon.sdf is a world that configures moon gravity, GUI, physics, a falling box, and then includes a terrain model.
• The Moon 60S 30S Fuel model is a regular model folder with model.config, model.sdf, and a .dae mesh. It is not loaded as a DEM heightmap at runtime.

This mesh approach can be more stable for rendering because Gazebo uses the normal mesh pipeline instead of the OGRE heightmap terrain shader.

Mesh model generator

./dem_to_mesh_model.py dem_500m_10m_utm.tif .

Output:

models/local_dem_terrain/model.config
models/local_dem_terrain/model.sdf
models/local_dem_terrain/meshes/dem_terrain.dae
dem_mesh_world.sdf

The generated model is structured like a Fuel model:

local_dem_terrain/
├── model.config
├── model.sdf
└── meshes/
    └── dem_terrain.dae

The world includes it with a model URI:

<include>
  <uri>model://local_dem_terrain</uri>
</include>

Before running the world, point Gazebo at the local model folder:

export GZ_SIM_RESOURCE_PATH="$PWD/models"
gz sim -r -v4 dem_mesh_world.sdf

The mesh generator converts the DEM grid into a centered terrain mesh:

• x and y are local meter coordinates from the UTM raster.
• z = elevation - minimum_elevation, so the lowest point is near z = 0.
• The generated .dae is used for both collision and visual geometry.
• The world uses DART with Bullet collision detection, matching the moon DEM example's heightmap-friendly physics setting.

Files

Download script

#!/usr/bin/env bash
set -euo pipefail

if [ "$#" -lt 2 ] || [ "$#" -gt 3 ]; then
    echo "Usage: $0 CENTER_LAT CENTER_LON [SIZE_M]" >&2
    echo "Example 500 m square: $0 48.7769 -121.8144 500" >&2
    exit 1
fi

if [ -z "${OPENTOPOGRAPHY_API_KEY:-}" ]; then
    echo "Set OPENTOPOGRAPHY_API_KEY before running this script." >&2
    exit 1
fi

CENTER_LAT="$1"
CENTER_LON="$2"
SIZE_M="${3:-2000}"
HALF_SIZE_M="$(
python3 - "$SIZE_M" <<'PY'
import sys

size_m = float(sys.argv[1])
if size_m <= 0:
    raise SystemExit("SIZE_M must be positive.")
print(size_m / 2.0)
PY
)"
SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
OUT_DIR="${SCRIPT_DIR}/raw"
SIZE_LABEL="$(python3 - "$SIZE_M" <<'PY'
import sys

size_m = float(sys.argv[1])
if size_m.is_integer():
    print(str(int(size_m)))
else:
    print(str(size_m).replace(".", "p"))
PY
)"
OUT_FILE="${OUT_DIR}/ot_usgs10m_${SIZE_LABEL}m.tif"

mkdir -p "$OUT_DIR"

read -r SOUTH NORTH WEST EAST < <(
python3 - "$CENTER_LAT" "$CENTER_LON" "$HALF_SIZE_M" <<'PY'
import math
import sys

lat = float(sys.argv[1])
lon = float(sys.argv[2])
half_size_m = float(sys.argv[3])

meters_per_degree_lat = 111_320.0
meters_per_degree_lon = meters_per_degree_lat * math.cos(math.radians(lat))

if abs(meters_per_degree_lon) < 1:
    raise SystemExit("Longitude scale is too small near the pole.")

dlat = half_size_m / meters_per_degree_lat
dlon = half_size_m / meters_per_degree_lon

print(f"{lat - dlat:.8f} {lat + dlat:.8f} {lon - dlon:.8f} {lon + dlon:.8f}")
PY
)

URL="https://portal.opentopography.org/API/usgsdem"

curl --fail --location --get "$URL" \
    --data-urlencode "datasetName=USGS10m" \
    --data-urlencode "south=${SOUTH}" \
    --data-urlencode "north=${NORTH}" \
    --data-urlencode "west=${WEST}" \
    --data-urlencode "east=${EAST}" \
    --data-urlencode "outputFormat=GTiff" \
    --data-urlencode "API_Key=${OPENTOPOGRAPHY_API_KEY}" \
    --output "$OUT_FILE"

echo "Downloaded: $OUT_FILE"
gdalinfo "$OUT_FILE" | sed -n '1,35p'

Prepare script

#!/usr/bin/env bash
set -euo pipefail

if [ "$#" -lt 2 ] || [ "$#" -gt 3 ]; then
    echo "Usage: $0 CENTER_LAT CENTER_LON [SIZE_M]" >&2
    echo "Example 500 m square: $0 48.7769 -121.8144 500" >&2
    exit 1
fi

CENTER_LAT="$1"
CENTER_LON="$2"
SIZE_M="${3:-2000}"
SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
SIZE_LABEL="$(python3 - "$SIZE_M" <<'PY'
import sys

size_m = float(sys.argv[1])
if size_m.is_integer():
    print(str(int(size_m)))
else:
    print(str(size_m).replace(".", "p"))
PY
)"
IN_FILE="${SCRIPT_DIR}/raw/ot_usgs10m_${SIZE_LABEL}m.tif"
OUT_FILE="${SCRIPT_DIR}/dem_${SIZE_LABEL}m_10m_utm.tif"

if [ ! -f "$IN_FILE" ]; then
    echo "Missing $IN_FILE. Run download_usgs10m_2km.sh with the same SIZE_M first." >&2
    exit 1
fi

UTM_EPSG="$(
python3 - "$CENTER_LAT" "$CENTER_LON" <<'PY'
import math
import sys

lat = float(sys.argv[1])
lon = float(sys.argv[2])
zone = int(math.floor((lon + 180.0) / 6.0) + 1)
epsg = 32600 + zone if lat >= 0 else 32700 + zone
print(epsg)
PY
)"

gdalwarp \
    -overwrite \
    -t_srs "EPSG:${UTM_EPSG}" \
    -tr 10 10 \
    -r bilinear \
    -of GTiff \
    -co COMPRESS=DEFLATE \
    "$IN_FILE" \
    "$OUT_FILE"

gdalinfo -stats "$OUT_FILE" > "${OUT_FILE%.tif}.info.txt"

echo "Prepared: $OUT_FILE"
echo "Projection: EPSG:${UTM_EPSG}"
sed -n '1,80p' "${OUT_FILE%.tif}.info.txt"

Mesh model generator

#!/usr/bin/python3
import math
import sys
import struct
import zlib
from pathlib import Path
from xml.sax.saxutils import escape

try:
    import numpy as np
    from osgeo import gdal
except ModuleNotFoundError as exc:
    raise SystemExit(
        "Missing Python dependency. Install it on Ubuntu with:\n"
        "  sudo apt install python3-gdal python3-numpy\n"
        "If you use a virtual environment, install GDAL bindings that match your system GDAL."
    ) from exc

gdal.UseExceptions()


def load_dem(path: Path) -> tuple[np.ndarray, float, float, float, float, float, float]:
    dataset = gdal.Open(str(path))
    if dataset is None:
        raise RuntimeError(f"Could not open DEM: {path}")

    band = dataset.GetRasterBand(1)
    data = band.ReadAsArray().astype(float)
    nodata = band.GetNoDataValue()
    valid = np.isfinite(data)
    if nodata is not None:
        valid &= data != nodata

    if not np.any(valid):
        raise RuntimeError("DEM has no valid elevation samples.")

    min_elevation = float(np.min(data[valid]))
    max_elevation = float(np.max(data[valid]))
    data[~valid] = min_elevation

    transform = dataset.GetGeoTransform()
    pixel_width = abs(float(transform[1]))
    pixel_height = abs(float(transform[5]))
    width_m = data.shape[1] * pixel_width
    depth_m = data.shape[0] * pixel_height

    return data, pixel_width, pixel_height, width_m, depth_m, min_elevation, max_elevation


def normal_for_triangle(a: np.ndarray, b: np.ndarray, c: np.ndarray) -> np.ndarray:
    normal = np.cross(b - a, c - a)
    length = np.linalg.norm(normal)
    if length == 0:
        return np.array([0.0, 0.0, 1.0])
    return normal / length


def vertex_normals(vertices: np.ndarray, triangles: list[tuple[int, int, int]]) -> np.ndarray:
    normals = np.zeros_like(vertices)
    for i0, i1, i2 in triangles:
        normal = normal_for_triangle(vertices[i0], vertices[i1], vertices[i2])
        normals[i0] += normal
        normals[i1] += normal
        normals[i2] += normal

    lengths = np.linalg.norm(normals, axis=1)
    zero = lengths == 0
    lengths[zero] = 1.0
    normals = normals / lengths[:, None]
    normals[zero] = np.array([0.0, 0.0, 1.0])
    return normals


def write_png(path: Path, width: int, height: int, rows: list[list[int]]) -> None:
    raw = b"".join(b"\x00" + bytes(row) for row in rows)

    def chunk(kind: bytes, data: bytes) -> bytes:
        crc = zlib.crc32(kind + data) & 0xFFFFFFFF
        return struct.pack(">I", len(data)) + kind + data + struct.pack(">I", crc)

    png = b"\x89PNG\r\n\x1a\n"
    png += chunk(b"IHDR", struct.pack(">IIBBBBB", width, height, 8, 2, 0, 0, 0))
    png += chunk(b"IDAT", zlib.compress(raw, 9))
    png += chunk(b"IEND", b"")
    path.write_bytes(png)


def terrain_color(value: float) -> tuple[int, int, int]:
    # A compact elevation color ramp: low green, mid tan, high light rock.
    stops = [
        (0.00, (62, 86, 55)),
        (0.35, (97, 122, 72)),
        (0.65, (146, 132, 92)),
        (1.00, (218, 214, 190)),
    ]
    for (left_pos, left), (right_pos, right) in zip(stops, stops[1:]):
        if value <= right_pos:
            t = (value - left_pos) / (right_pos - left_pos)
            return tuple(round(left[i] + t * (right[i] - left[i])) for i in range(3))
    return stops[-1][1]


def write_texture(path: Path, elevation: np.ndarray, min_elevation: float, max_elevation: float) -> None:
    span = max(max_elevation - min_elevation, 1.0)
    normalized = np.clip((elevation - min_elevation) / span, 0.0, 1.0)
    rows = []
    for row in normalized:
        png_row = []
        for value in row:
            png_row.extend(terrain_color(float(value)))
        rows.append(png_row)
    write_png(path, elevation.shape[1], elevation.shape[0], rows)


def write_dae(
    path: Path,
    vertices: np.ndarray,
    texcoords: np.ndarray,
    triangles: list[tuple[int, int, int]],
    texture_file: str,
) -> None:
    normals = vertex_normals(vertices, triangles)
    position_values = " ".join(f"{value:.6f}" for value in vertices.reshape(-1))
    normal_values = " ".join(f"{value:.6f}" for value in normals.reshape(-1))
    texcoord_values = " ".join(f"{value:.6f}" for value in texcoords.reshape(-1))

    # Each triangle uses the same index for VERTEX, NORMAL, and TEXCOORD.
    triangle_indices = " ".join(
        f"{i0} {i0} {i0} {i1} {i1} {i1} {i2} {i2} {i2}" for i0, i1, i2 in triangles
    )

    path.write_text(f'''<?xml version="1.0" encoding="utf-8"?>
<COLLADA xmlns="http://www.collada.org/2005/11/COLLADASchema" version="1.4.1">
  <asset>
    <contributor>
      <authoring_tool>dem_to_mesh_model.py</authoring_tool>
    </contributor>
    <unit name="meter" meter="1"/>
    <up_axis>Z_UP</up_axis>
  </asset>

  <library_effects>
    <effect id="terrain_effect">
      <profile_COMMON>
        <newparam sid="terrain_texture_surface">
          <surface type="2D">
            <init_from>terrain_texture_image</init_from>
          </surface>
        </newparam>
        <newparam sid="terrain_texture_sampler">
          <sampler2D>
            <source>terrain_texture_surface</source>
          </sampler2D>
        </newparam>
        <technique sid="common">
          <lambert>
            <diffuse>
              <texture texture="terrain_texture_sampler" texcoord="UVMap"/>
            </diffuse>
          </lambert>
        </technique>
      </profile_COMMON>
    </effect>
  </library_effects>

  <library_images>
    <image id="terrain_texture_image" name="terrain_texture_image">
      <init_from>{escape(texture_file)}</init_from>
    </image>
  </library_images>

  <library_materials>
    <material id="terrain_material" name="terrain_material">
      <instance_effect url="#terrain_effect"/>
    </material>
  </library_materials>

  <library_geometries>
    <geometry id="dem_terrain_mesh" name="dem_terrain_mesh">
      <mesh>
        <source id="dem_terrain_positions">
          <float_array id="dem_terrain_positions_array" count="{vertices.size}">{position_values}</float_array>
          <technique_common>
            <accessor source="#dem_terrain_positions_array" count="{len(vertices)}" stride="3">
              <param name="X" type="float"/>
              <param name="Y" type="float"/>
              <param name="Z" type="float"/>
            </accessor>
          </technique_common>
        </source>
        <source id="dem_terrain_normals">
          <float_array id="dem_terrain_normals_array" count="{normals.size}">{normal_values}</float_array>
          <technique_common>
            <accessor source="#dem_terrain_normals_array" count="{len(normals)}" stride="3">
              <param name="X" type="float"/>
              <param name="Y" type="float"/>
              <param name="Z" type="float"/>
            </accessor>
          </technique_common>
        </source>
        <source id="dem_terrain_uvmap">
          <float_array id="dem_terrain_uvmap_array" count="{texcoords.size}">{texcoord_values}</float_array>
          <technique_common>
            <accessor source="#dem_terrain_uvmap_array" count="{len(texcoords)}" stride="2">
              <param name="S" type="float"/>
              <param name="T" type="float"/>
            </accessor>
          </technique_common>
        </source>
        <vertices id="dem_terrain_vertices">
          <input semantic="POSITION" source="#dem_terrain_positions"/>
        </vertices>
        <triangles material="terrain_material" count="{len(triangles)}">
          <input semantic="VERTEX" source="#dem_terrain_vertices" offset="0"/>
          <input semantic="NORMAL" source="#dem_terrain_normals" offset="1"/>
          <input semantic="TEXCOORD" source="#dem_terrain_uvmap" offset="2" set="0"/>
          <p>{triangle_indices}</p>
        </triangles>
      </mesh>
    </geometry>
  </library_geometries>

  <library_visual_scenes>
    <visual_scene id="Scene" name="Scene">
      <node id="dem_terrain" name="dem_terrain">
        <instance_geometry url="#dem_terrain_mesh">
          <bind_material>
            <technique_common>
              <instance_material symbol="terrain_material" target="#terrain_material">
                <bind_vertex_input semantic="UVMap" input_semantic="TEXCOORD" input_set="0"/>
              </instance_material>
            </technique_common>
          </bind_material>
        </instance_geometry>
      </node>
    </visual_scene>
  </library_visual_scenes>

  <scene>
    <instance_visual_scene url="#Scene"/>
  </scene>
</COLLADA>
''', encoding="utf-8")


def write_model_config(path: Path) -> None:
    path.write_text('''<?xml version="1.0" ?>
<model>
  <name>Local DEM Terrain Mesh</name>
  <version>1.0</version>
  <sdf version="1.9">model.sdf</sdf>

  <author>
    <name>new_blog demo</name>
  </author>

  <description>
    Mesh terrain generated from a local DEM GeoTIFF.
  </description>
</model>
''', encoding="utf-8")


def write_model_sdf(path: Path) -> None:
    path.write_text('''<?xml version="1.0" ?>
<sdf version="1.9">
  <model name="local_dem_terrain_mesh">
    <static>true</static>
    <link name="link">
      <collision name="collision">
        <geometry>
          <mesh>
            <uri>meshes/dem_terrain.dae</uri>
          </mesh>
        </geometry>
      </collision>
      <visual name="visual">
        <cast_shadows>true</cast_shadows>
        <geometry>
          <mesh>
            <uri>meshes/dem_terrain.dae</uri>
          </mesh>
        </geometry>
      </visual>
    </link>
  </model>
</sdf>
''', encoding="utf-8")


def write_world(path: Path, width_m: float, depth_m: float, max_z: float) -> None:
    camera_x = -0.65 * width_m
    camera_y = -0.75 * depth_m
    camera_z = max(max_z + 0.9 * max(width_m, depth_m), 500.0)
    box_z = max_z + 50.0
    path.write_text(f'''<?xml version="1.0" ?>
<sdf version="1.9">
  <world name="dem_mesh_world">
    <gravity>0 0 -9.8</gravity>

    <physics name="1ms" type="ignored">
      <max_step_size>0.001</max_step_size>
      <real_time_factor>1.0</real_time_factor>
      <dart>
        <collision_detector>bullet</collision_detector>
      </dart>
    </physics>

    <plugin
      filename="gz-sim-physics-system"
      name="gz::sim::systems::Physics">
    </plugin>
    <plugin
      filename="gz-sim-user-commands-system"
      name="gz::sim::systems::UserCommands">
    </plugin>
    <plugin
      filename="gz-sim-scene-broadcaster-system"
      name="gz::sim::systems::SceneBroadcaster">
    </plugin>

    <gui fullscreen="0">
      <plugin filename="MinimalScene" name="3D View">
        <gz-gui>
          <title>3D View</title>
          <property type="bool" key="showTitleBar">false</property>
          <property type="string" key="state">docked</property>
        </gz-gui>
        <engine>ogre2</engine>
        <scene>scene</scene>
        <ambient_light>0.35 0.35 0.35</ambient_light>
        <background_color>0.78 0.82 0.86</background_color>
        <camera_pose>{camera_x:.3f} {camera_y:.3f} {camera_z:.3f} 0 0.95 0.72</camera_pose>
      </plugin>
    </gui>

    <light type="directional" name="sun">
      <cast_shadows>true</cast_shadows>
      <pose>0 0 1200 0 0 0</pose>
      <diffuse>0.9 0.86 0.78 1</diffuse>
      <specular>0.25 0.25 0.25 1</specular>
      <direction>-0.45 0.35 -0.82</direction>
    </light>

    <include>
      <uri>model://local_dem_terrain</uri>
    </include>

    <model name="falling_box">
      <pose>0 0 {box_z:.3f} 0 0 0</pose>
      <link name="box_link">
        <inertial>
          <mass>1.0</mass>
          <inertia>
            <ixx>0.16666</ixx>
            <iyy>0.16666</iyy>
            <izz>0.16666</izz>
          </inertia>
        </inertial>
        <collision name="box_collision">
          <geometry>
            <box>
              <size>10 10 10</size>
            </box>
          </geometry>
        </collision>
        <visual name="box_visual">
          <geometry>
            <box>
              <size>10 10 10</size>
            </box>
          </geometry>
          <material>
            <ambient>1 0.1 0.1 1</ambient>
            <diffuse>1 0.1 0.1 1</diffuse>
            <specular>0.4 0.4 0.4 1</specular>
          </material>
        </visual>
      </link>
    </model>
  </world>
</sdf>
''', encoding="utf-8")


def main() -> int:
    if len(sys.argv) != 3:
        print(f"Usage: {sys.argv[0]} DEM_TIF OUT_CODE_DIR", file=sys.stderr)
        return 1

    dem_path = Path(sys.argv[1]).resolve()
    out_dir = Path(sys.argv[2]).resolve()
    model_dir = out_dir / "models" / "local_dem_terrain"
    mesh_dir = model_dir / "meshes"
    mesh_dir.mkdir(parents=True, exist_ok=True)

    data, pixel_width, pixel_height, width_m, depth_m, min_elevation, max_elevation = load_dem(dem_path)
    rows, cols = data.shape
    z = data - min_elevation
    max_z = max_elevation - min_elevation

    vertices = []
    texcoords = []
    for row in range(rows):
        y = depth_m / 2.0 - row * pixel_height
        for col in range(cols):
            x = col * pixel_width - width_m / 2.0
            vertices.append((x, y, z[row, col]))
            u = col / (cols - 1) if cols > 1 else 0.0
            v = 1.0 - row / (rows - 1) if rows > 1 else 0.0
            texcoords.append((u, v))
    vertices_np = np.array(vertices, dtype=float)
    texcoords_np = np.array(texcoords, dtype=float)

    triangles = []
    for row in range(rows - 1):
        for col in range(cols - 1):
            i00 = row * cols + col
            i10 = row * cols + col + 1
            i01 = (row + 1) * cols + col
            i11 = (row + 1) * cols + col + 1
            triangles.append((i00, i01, i10))
            triangles.append((i10, i01, i11))

    texture_file = "dem_texture.png"
    write_texture(mesh_dir / texture_file, data, min_elevation, max_elevation)
    write_dae(mesh_dir / "dem_terrain.dae", vertices_np, texcoords_np, triangles, texture_file)
    write_model_config(model_dir / "model.config")
    write_model_sdf(model_dir / "model.sdf")
    write_world(out_dir / "dem_mesh_world.sdf", width_m, depth_m, max_z)

    print(f"Model written: {model_dir}")
    print(f"World written: {out_dir / 'dem_mesh_world.sdf'}")
    print(f"DEM samples: {cols} x {rows}")
    print(f"Mesh vertices: {len(vertices_np)}")
    print(f"Mesh triangles: {len(triangles)}")
    print(f"World size: {width_m:.3f} x {depth_m:.3f} x {max_z:.3f} m")
    return 0


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

World generator

#!/usr/bin/env python3
import json
import math
import subprocess
import sys
import struct
import zlib
from pathlib import Path
from xml.sax.saxutils import escape


def run_gdalinfo(path: Path) -> dict:
    result = subprocess.run(
        ["gdalinfo", "-json", "-stats", str(path)],
        check=True,
        text=True,
        stdout=subprocess.PIPE,
    )
    return json.loads(result.stdout)


def band_min_max(info: dict) -> tuple[float, float]:
    band = info["bands"][0]
    metadata = band.get("metadata", {})
    stats = metadata.get("", {})

    min_value = stats.get("STATISTICS_MINIMUM", band.get("minimum"))
    max_value = stats.get("STATISTICS_MAXIMUM", band.get("maximum"))

    if min_value is None or max_value is None:
        raise RuntimeError("DEM statistics are missing. Run gdalinfo -stats on the DEM.")

    return float(min_value), float(max_value)


def pixel_size(info: dict) -> tuple[float, float]:
    transform = info["geoTransform"]
    return abs(float(transform[1])), abs(float(transform[5]))


def write_png(path: Path, width: int, height: int, rows: list[list[int]]) -> None:
    raw = b"".join(b"\x00" + bytes(row) for row in rows)

    def chunk(kind: bytes, data: bytes) -> bytes:
        crc = zlib.crc32(kind + data) & 0xFFFFFFFF
        return struct.pack(">I", len(data)) + kind + data + struct.pack(">I", crc)

    png = b"\x89PNG\r\n\x1a\n"
    png += chunk(b"IHDR", struct.pack(">IIBBBBB", width, height, 8, 2, 0, 0, 0))
    png += chunk(b"IDAT", zlib.compress(raw, 9))
    png += chunk(b"IEND", b"")
    path.write_bytes(png)


def ensure_texture_assets(output_dir: Path) -> None:
    texture_dir = output_dir / "textures"
    texture_dir.mkdir(parents=True, exist_ok=True)

    diffuse = texture_dir / "dem_diffuse.png"
    if not diffuse.exists():
        rows = []
        for y in range(64):
            row = []
            for x in range(64):
                noise = ((x * 31 + y * 17 + (x * y) % 19) % 33) - 16
                row.extend([
                    max(0, min(255, 105 + noise)),
                    max(0, min(255, 120 + noise)),
                    max(0, min(255, 82 + noise // 2)),
                ])
            rows.append(row)
        write_png(diffuse, 64, 64, rows)

    normal = texture_dir / "flat_normal.png"
    if not normal.exists():
        write_png(normal, 8, 8, [[128, 128, 255] * 8 for _ in range(8)])


def world_sdf(dem_path: Path, width_m: float, depth_m: float, z_range: float, z_pos: float) -> str:
    dem_uri = escape(dem_path.name)
    return f'''<?xml version="1.0" ?>
<sdf version="1.9">
  <world name="dem_terrain_demo">
    <gravity>0 0 -9.8</gravity>

    <physics name="1ms" type="ignored">
      <max_step_size>0.001</max_step_size>
      <real_time_factor>1.0</real_time_factor>
    </physics>

    <plugin
      filename="gz-sim-physics-system"
      name="gz::sim::systems::Physics">
    </plugin>
    <plugin
      filename="gz-sim-user-commands-system"
      name="gz::sim::systems::UserCommands">
    </plugin>
    <plugin
      filename="gz-sim-scene-broadcaster-system"
      name="gz::sim::systems::SceneBroadcaster">
    </plugin>

    <light type="directional" name="sun">
      <cast_shadows>true</cast_shadows>
      <pose>0 0 1200 0 0 0</pose>
      <diffuse>0.9 0.86 0.78 1</diffuse>
      <specular>0.25 0.25 0.25 1</specular>
      <direction>-0.45 0.35 -0.82</direction>
    </light>

    <model name="dem_terrain">
      <static>true</static>
      <link name="terrain_link">
        <collision name="terrain_collision">
          <geometry>
            <heightmap>
              <uri>{dem_uri}</uri>
              <size>{width_m:.3f} {depth_m:.3f} {z_range:.3f}</size>
              <pos>0 0 {z_pos:.3f}</pos>
              <sampling>1</sampling>
            </heightmap>
          </geometry>
        </collision>
        <visual name="terrain_visual">
          <geometry>
            <heightmap>
              <uri>{dem_uri}</uri>
              <size>{width_m:.3f} {depth_m:.3f} {z_range:.3f}</size>
              <pos>0 0 {z_pos:.3f}</pos>
              <sampling>1</sampling>
              <texture>
                <diffuse>textures/dem_diffuse.png</diffuse>
                <normal>textures/flat_normal.png</normal>
                <size>50</size>
              </texture>
            </heightmap>
          </geometry>
        </visual>
      </link>
    </model>

    <model name="origin_marker">
      <static>true</static>
      <pose>0 0 5 0 0 0</pose>
      <link name="link">
        <visual name="visual">
          <geometry>
            <sphere>
              <radius>5</radius>
            </sphere>
          </geometry>
          <material>
            <ambient>1 0.1 0.1 1</ambient>
            <diffuse>1 0.1 0.1 1</diffuse>
          </material>
        </visual>
      </link>
    </model>
  </world>
</sdf>
'''


def main() -> int:
    if len(sys.argv) != 3:
        print(f"Usage: {sys.argv[0]} DEM_TIF OUT_WORLD_SDF", file=sys.stderr)
        return 1

    dem_path = Path(sys.argv[1]).resolve()
    out_path = Path(sys.argv[2]).resolve()

    info = run_gdalinfo(dem_path)
    raster_width, raster_height = info["size"]
    pixel_width, pixel_height = pixel_size(info)
    min_elevation, max_elevation = band_min_max(info)

    width_m = raster_width * pixel_width
    depth_m = raster_height * pixel_height
    z_range = max_elevation - min_elevation
    if not math.isfinite(z_range) or z_range <= 0:
        raise RuntimeError("DEM elevation range must be positive.")

    z_pos = -max_elevation

    ensure_texture_assets(out_path.parent)
    out_path.write_text(
        world_sdf(dem_path, width_m, depth_m, z_range, z_pos),
        encoding="utf-8",
    )

    print(f"Wrote: {out_path}")
    print(f"DEM size: {raster_width} x {raster_height} pixels")
    print(f"World size: {width_m:.3f} x {depth_m:.3f} x {z_range:.3f} m")
    print(f"Elevation: min={min_elevation:.3f} max={max_elevation:.3f}")
    print(f"Terrain pos z: {z_pos:.3f}")
    return 0


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

If 10 m is not available

If OpenTopography returns an error for USGS10m, the selected point is probably outside USGS 3DEP coverage or the account does not have access.

For non-U.S. terrain, use a 30 m global dataset. Example API endpoint:

https://portal.opentopography.org/API/globaldem?demtype=COP30&south=...&north=...&west=...&east=...&outputFormat=GTiff&API_Key=...

For a 0.5 km x 0.5 km world, 30 m DEM gives about 17 x 17 samples. For a 2 km x 2 km world, it gives about 67 x 67 samples. It is still useful for large terrain shape, but it will not show small local features.

Common problems

• If Gazebo cannot load the .tif, run gdalinfo dem_500m_10m_utm.tif first. If GDAL cannot read it, Gazebo will not read it.
• If the terrain is huge or tiny, check that the DEM was reprojected to UTM meters before generating the world.
• If the terrain is visually offset from collision, make sure the same <heightmap> values are used in both <collision> and <visual>.
• If objects do not fall, make sure simulation is running. Use gz sim -r -v4 dem_world.sdf or press the play button in the GUI. Also check that the inserted model does not have <static>true</static>.
• If OGRE crashes with a fragment shader error like detailCol0 undeclared, make sure the visual heightmap has a <texture> block with valid diffuse and normal texture files.
• If simulation is slow, resample to 20 m or 30 m with gdalwarp -tr 20 20.
• If the terrain is hard to find in the GUI, check the generated <pos> z value.
• If dem_mesh_world.sdf cannot find model://local_dem_terrain, set GZ_SIM_RESOURCE_PATH="$PWD/models" from the folder that contains the generated models directory.

References

• Gazebo Sim DEM heightmap tutorial
• Gazebo dem_moon.sdf example
• SDFormat heightmap spec
• Moon 60S 30S Fuel model
• OpenTopography API access to USGS 3DEP rasters
• OpenTopography API developer notes