"raster-dem" source doesn't work with custom dem tiles.

[starwarcat]

mapbox-gl-js version: v2.8.2

browser: Any

Steps to Trigger Behavior

1. I am trying to use my own 1 meter high resolution RGB tiles to setTerrain, rather than using mapbox://mapbox.mapbox-terrain-dem-v1, for a small area.
2. The code is:
   map.addSource("rwanda-dem", { "type": "raster-dem", "tiles": ["/public/images/raster-dem/{z}/{x}/{y}.webp"], "tileSize": 514, "maxzoom": 18, "minzoom": 1, });
3. The area's terrain was generated correctly with elevation, excep the whole area is floating from other part of the map where no elevation covered.

[fr-an-k]

Haven't tried myself yet, but supposedly it needs this pixel encoding:
https://github.com/tilezen/joerd/blob/master/docs/formats.md#terrarium
Possibly limited to PNG

[abrman]

I just fixed something like this an hour ago.

What is happening is that the entire map is at sea level as it has no height data, your data probably accurate hovers in the 150-500m range? What you need to do is adjust your entire heightmap so that ideally your lowest point is at 0. I'm assuming you're following the elevation = -10000 + (({R} * 256 * 256 + {G} * 256 + {B}) * 0.1) formula.

So basically you get the elevation of your lowest point and subtract it from your tile-set.

[mttreed]

Facing the same issue.

@abrman thanks for the tip, however, I tried this offset trick and while it may have helped a bit, the issue persists. See screenshots and python code used for RGB encoding with the offset. For context I am trying to use the Swiss government's high resolution DEMs inside Mapbox GL JS and Mapbox iOS SDK.

Mapbox GL JS:

iOS:
  
  

Any ideas?

#!/usr/bin/env python3
import numpy as np
from osgeo import gdal
import os
import sys

# Input and output files
input_file = "test_small_area_3857.tif"
output_file = "test_small_area_terrain_rgb_offset.tif"  # New filename to avoid overwriting

print(f"Processing {input_file} to Mapbox Terrain-RGB format with elevation offset...")

# Open the input DEM
ds = gdal.Open(input_file)
if ds is None:
    print(f"Error: Could not open {input_file}")
    sys.exit(1)

# Read elevation data
band = ds.GetRasterBand(1)
nodata = band.GetNoDataValue()
elevation = band.ReadAsArray().astype(np.float32)

# Get statistics on elevation data (ignoring NoData values)
valid_mask = elevation != nodata if nodata is not None else np.ones_like(elevation, dtype=bool)
if np.any(valid_mask):
    min_elev = np.min(elevation[valid_mask])
    max_elev = np.max(elevation[valid_mask])
    mean_elev = np.mean(elevation[valid_mask])
    print(f"Original elevation statistics:")
    print(f"  Minimum: {min_elev:.2f}m")
    print(f"  Maximum: {max_elev:.2f}m")
    print(f"  Mean: {mean_elev:.2f}m")
    print(f"  Range: {max_elev - min_elev:.2f}m")
    print(f"  NoData value: {nodata}")
else:
    print("Warning: No valid elevation data found!")
    sys.exit(1)

# Replace NoData values with a suitable elevation
if nodata is not None:
    nodata_count = np.sum(~valid_mask)
    if nodata_count > 0:
        print(f"Replacing {nodata_count} NoData pixels with 0m")
        elevation[~valid_mask] = 0

# ===== NEW STEP: Apply elevation offset =====
# Store original elevation for comparison
elevation_original = elevation.copy()

# Offset elevations by subtracting the minimum elevation
# This makes the lowest point effectively at 0m elevation
print(f"\n==== APPLYING ELEVATION OFFSET ====")
print(f"Subtracting minimum elevation ({min_elev:.2f}m) from all values")
elevation_offset = elevation - min_elev

# Calculate new statistics after offset
min_offset = np.min(elevation_offset[valid_mask])
max_offset = np.max(elevation_offset[valid_mask])
print(f"Offset elevation statistics:")
print(f"  New minimum: {min_offset:.2f}m")
print(f"  New maximum: {max_offset:.2f}m")
print(f"  New range: {max_offset - min_offset:.2f}m")
print(f"  Relative relief preserved: {(max_offset - min_offset) == (max_elev - min_elev)}")

# Use the offset elevation for further processing
elevation = elevation_offset
# ===== END OF NEW STEP =====

# Clip elevation values to Terrain-RGB valid range
original_min = np.min(elevation)
original_max = np.max(elevation)
elevation = np.clip(elevation, -10000, 25500)
clipped_min = np.min(elevation)
clipped_max = np.max(elevation)

if original_min != clipped_min or original_max != clipped_max:
    print(f"Warning: Elevation values were clipped to fit Terrain-RGB range:")
    print(f"  Original range: {original_min:.2f}m to {original_max:.2f}m")
    print(f"  Clipped range: {clipped_min:.2f}m to {clipped_max:.2f}m")
else:
    print("All elevation values within valid Terrain-RGB range.")

# Convert elevation values to Mapbox's Terrain-RGB encoding
# Step 1: Add 10000 to shift the range (to handle negative values)
# Step 2: Multiply by 10 for 0.1m precision
# Step 3: Calculate the RGB components
print("Converting offset elevation values to RGB channels...")
terrain_values = np.floor((elevation + 10000) / 0.1).astype(np.uint32)

# Extract individual RGB channels
r = np.floor(terrain_values / (256 * 256)).astype(np.uint8)
g = np.floor((terrain_values - r * 256 * 256) / 256).astype(np.uint8)
b = (terrain_values - r * 256 * 256 - g * 256).astype(np.uint8)

# Verify conversion with sample points
sample_points = [
    (elevation.shape[0] // 4, elevation.shape[1] // 4),
    (elevation.shape[0] // 2, elevation.shape[1] // 2),
    (3 * elevation.shape[0] // 4, 3 * elevation.shape[1] // 4)
]

print("\nVerification of encoding at sample points:")
for i, (y, x) in enumerate(sample_points):
    # Get both original and offset elevations for comparison
    orig_elev = elevation_original[y, x]
    offset_elev = elevation[y, x]
    rgb = (r[y, x], g[y, x], b[y, x])
    
    # Decode to verify
    decoded_elev = -10000 + ((rgb[0] * 256 * 256 + rgb[1] * 256 + rgb[2]) * 0.1)
    
    print(f"Sample {i+1} at position ({x}, {y}):")
    print(f"  Original absolute elevation: {orig_elev:.2f}m")
    print(f"  Offset elevation: {offset_elev:.2f}m")
    print(f"  RGB values: R={rgb[0]}, G={rgb[1]}, B={rgb[2]}")
    print(f"  Decoded elevation: {decoded_elev:.2f}m")
    print(f"  Difference: {offset_elev - decoded_elev:.6f}m")

# Create the RGB output file
print(f"\nCreating output file: {output_file}")
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(
    output_file,
    ds.RasterXSize,
    ds.RasterYSize,
    3,  # 3 bands for RGB
    gdal.GDT_Byte,
    options=[
        'COMPRESS=DEFLATE',  # Good compression with no quality loss
        'PREDICTOR=2',       # Horizontal predictor for better compression
        'TILED=YES',         # Optimize for tiled reading
        'BLOCKXSIZE=256',    # Standard tile size for web mapping
        'BLOCKYSIZE=256'
    ]
)

# Copy the geospatial reference information
out_ds.SetGeoTransform(ds.GetGeoTransform())
out_ds.SetProjection(ds.GetProjection())

# Write the RGB bands
out_ds.GetRasterBand(1).WriteArray(r)  # Red channel
out_ds.GetRasterBand(2).WriteArray(g)  # Green channel
out_ds.GetRasterBand(3).WriteArray(b)  # Blue channel

# Clean up
out_ds.FlushCache()
out_ds = None
ds = None

print(f"Successfully created offset Terrain-RGB file: {output_file}")
print("This file is now ready for upload to Mapbox.")
print("\nIMPORTANT: This terrain data has been offset so that the lowest point")
print("is at 0m elevation. This eliminates the dramatic cliff effect when")
print("transitioning to areas without custom terrain data.")