Cesium encountered an error while loading custom terrain

[wenqin]

I created terrain tiles using this script：

import gzip
import json
import os
import tempfile

import numpy as np
import rasterio as rio
from quantized_mesh_encoder import encode
from rasterio.transform import from_bounds
from rasterio.warp import calculate_default_transform, reproject, Resampling

from models.is_intersect_contains import Rectangle, check_relationship, get_intersection

# ---------- 参数 ----------
DEM_FILE = r'D:\tmp\dem\input\DEM_4326.tif'  # 原始 DEM
OUT_DIR = r'D:\tmp\dem\output'  # 输出根目录
MIN_ZOOM = 0  # 最小层级
MAX_ZOOM = 12  # 最大层级（12 级约 50 m 分辨率）
TARGET_SIZE = 65
# ---------------------------

WGS84 = {'init': 'EPSG:4326'}


def ensure_wgs84(src_path):
    """返回 WGS84 的临时路径（如已是 4326 则原样返回）"""
    with rio.open(src_path) as src:
        if src.crs == 'EPSG:4326':
            return src_path
        tmp = tempfile.NamedTemporaryFile(suffix='.tif', delete=False).name
        transform, width, height = calculate_default_transform(src.crs, WGS84, src.width, src.height, *src.bounds)
        meta = src.meta.copy()
        meta.update({'crs': WGS84, 'transform': transform, 'width': width, 'height': height})
        with rio.open(tmp, 'w', **meta) as dst:
            for i in range(1, src.count + 1):
                reproject(source=rio.band(src, i),
                          destination=rio.band(dst, i),
                          src_transform=src.transform,
                          src_crs=src.crs,
                          dst_transform=transform,
                          dst_crs=WGS84,
                          resampling=Resampling.bilinear)
        return tmp


def tile_bounds(x, y, z):
    """返回某TMS瓦片的经纬度边界 [west, south, east, north]"""
    n_col = 2.0 ** (z + 1)
    n_row = 2.0 ** z
    x_step_col = 360.0 / n_col
    y_step_row = 180.0 / n_row
    x_start_west = -180.0 + x_step_col * x
    x_end_east = x_start_west + x_step_col
    y_start_south = -90.0 + y_step_row * y
    y_end_north = y_start_south + y_step_row
    # print(f'自己计算的值：{x_start_west, y_start_south, x_end_east, y_end_north}')
    return [x_start_west, y_start_south, x_end_east, y_end_north]


def raster_range_in_TMS(global_bounds, z):
    min_lon, min_lat, max_lon, max_lat = global_bounds
    n_col = 2 ** (z + 1)
    n_row = 2 ** z
    x_step_col = 360 / n_col
    y_step_row = 180 / n_row
    x_min = int((180 + min_lon) // x_step_col)
    y_min = int((90 + min_lat) // y_step_row)
    x_max = int((180 + max_lon) // x_step_col)
    y_max = int((90 + max_lat) // y_step_row)
    return x_min, y_min, x_max, y_max


# 裁剪与重采样
def clip_and_resample(dataset, tile_bounds, data_bounds, target_size):
    t_rect = Rectangle(tile_bounds[0], tile_bounds[1], tile_bounds[2], tile_bounds[3])
    r_rect = Rectangle(data_bounds[0], data_bounds[1], data_bounds[2], data_bounds[3])
    check_result = check_relationship(t_rect, r_rect)
    # 切片边界包含数据边界或二者相同
    if check_result == 1 or check_result == 2:
        window = dataset.window(data_bounds[0], data_bounds[1], data_bounds[2], data_bounds[3])
        dst_array = read_tif(dataset, window, tile_bounds, target_size)
    # 数据边界包含切片边界或与二者相交
    elif check_result == 3 or check_result == 5:
        new_bounds = get_intersection(tile_bounds, data_bounds)
        window = dataset.window(new_bounds[0], new_bounds[1], new_bounds[2], new_bounds[3])
        dst_array = read_tif(dataset, window, tile_bounds, target_size)
    else:
        return None
    return dst_array


def read_tif(dataset, window, tile_bounds, target_size):
    # print(f'数据的nodata值为：{dataset.nodata}')
    data = dataset.read(1, window=window, boundless=True, fill_value=dataset.nodata)
    if dataset.nodata is not None:
        data = np.where(np.isclose(data, dataset.nodata), 0, data)
    data = np.where(data < -20000, 0, data)
    if data.size == 0 or np.all((data == 0)):
        return None
    # 目标transform（覆盖整个瓦片）
    dst_transform = from_bounds(tile_bounds[0],
                                tile_bounds[1],
                                tile_bounds[2],
                                tile_bounds[3],
                                target_size,
                                target_size)
    # 创建一个纯零数组，行列为target_size，数据类型为32位浮点型
    dst_array = np.zeros((target_size, target_size), dtype=np.float32)
    reproject(
        source=data,
        destination=dst_array,
        src_transform=dataset.window_transform(window),
        src_crs=dataset.crs,
        dst_transform=dst_transform,
        dst_crs=dataset.crs,
        resampling=Resampling.bilinear
    )
    dst_array = np.nan_to_num(dst_array, nan=0.0)
    return dst_array


def build_layer_json(out_dir, minz, maxz, bounds, available_list):
    """生成 Cesium 必需的 layer.json"""
    layer = {
        "tilejson": "1.0",
        "name": "terrain",
        "version": "1.0.0",
        "scheme": "tms",
        "tiles": ["{z}/{x}/{y}.terrain"],
        "minzoom": minz,
        "maxzoom": maxz,
        "bounds": bounds,
        # "valid_bounds": bounds,
        "projection": "EPSG:4326",
        "format": "quantized-mesh-1.0",
        # "metadataAvailability": 10,
        "available": [available_list]
    }
    with open(os.path.join(out_dir, 'layer.json'), 'w') as f:
        json.dump(layer, f, indent=2)


def numpy_grid_to_terrain(height_map, bounds):
    rows, cols = height_map.shape
    west, south, east, north = bounds

    # --- 1. Positions (不变) ---
    lons = np.linspace(west, east, cols)
    lats = np.linspace(north, south, rows)
    xx, yy = np.meshgrid(lons, lats)

    positions = np.column_stack((
        xx.flatten(),
        yy.flatten(),
        height_map.flatten()
    ))

    # --- 2. Indices (改进版：生成 N行3列) ---

    # 这里的 rows/cols 是顶点数，格子数要减 1
    h_cells = rows - 1
    w_cells = cols - 1

    # 生成左上角顶点的索引网格
    # 这里的 i, j 对应格子的行和列
    i = np.arange(h_cells)
    j = np.arange(w_cells)
    jj, ii = np.meshgrid(j, i)

    # 计算每个格子左上角的全局索引 k
    # k 是一个 (h_cells, w_cells) 的矩阵
    k = ii * cols + jj

    # 定义每个格子的四个顶点索引 (利用广播机制)
    # top_left = k
    # top_right = k + 1
    # bottom_left = k + cols
    # bottom_right = k + cols + 1

    # 构建两个三角形 (逆时针 CCW)
    # Triangle 1: TopLeft -> BottomLeft -> TopRight
    # 结果 shape: (h_cells, w_cells, 3)
    t1 = np.stack([k, k + cols, k + 1], axis=-1)

    # Triangle 2: TopRight -> BottomLeft -> BottomRight
    # 结果 shape: (h_cells, w_cells, 3)
    t2 = np.stack([k + 1, k + cols, k + cols + 1], axis=-1)

    # --- 关键步骤：交替合并 ---
    # 我们希望结果是：[Cell0_T1, Cell0_T2, Cell1_T1, Cell1_T2, ...]
    # 使用 dstack 将它们堆叠在一起，然后 reshape

    # combined shape: (h_cells, w_cells, 2, 3) -> 每个格子有2个三角形，每个三角形3个点
    combined = np.stack([t1, t2], axis=2)

    # Reshape 成 (N, 3)
    # N = 总格子数 * 2
    indices = combined.reshape(-1, 3)

    # 确保是整数类型 (uint32 或 uint16)
    indices = indices.astype(np.uint32)

    # --- 3. 编码 ---
    print(f"Indices Shape: {indices.shape}")  # 现在这里会显示 (8192, 3) 之类的
    print(f"第一行三角形: {indices[0]}")  # 方便调试查看

    with open('output_v2.terrain', 'wb') as f:
        # ⚠️ 注意：如果 encode 函数报错，说明它强制要求 1D 数组。
        # 此时只需在传参时加 .flatten() 即可，不影响上面的逻辑构建。
        encode(f, positions, indices.flatten())

    return positions, indices


def main():
    os.makedirs(OUT_DIR, exist_ok=True)
    wgs84_path = ensure_wgs84(DEM_FILE)
    with rio.open(wgs84_path) as ds:
        global_bounds = [ds.bounds.left, ds.bounds.bottom, ds.bounds.right, ds.bounds.top]
        print(f'栅格范围是：{global_bounds}')
        available_list = []
        for z in range(MIN_ZOOM, MAX_ZOOM + 1):
            print(f'=== zoom {z} ===')
            # 计算当前层级瓦片范围
            x_min, y_min, x_max, y_max = raster_range_in_TMS(global_bounds, z)
            print(f'X范围：{x_min,x_max},Y范围{y_min,y_max}')
            available_dict = {"startX": x_min, "startY": y_min, "endX": x_max, "endY": y_max}
            available_list.append(available_dict)
            tile_count = 0
            for x in range(x_min, x_max + 1):
                for y in range(y_min, y_max + 1):
                    # 根据当前层级的x,y,x获取当前的切片的实际范围。
                    bounds = tile_bounds(x, y, z)
                    # 判断切片范围与数据范围的位置关系，获取正确的“原始数据数组”
                    grid_data = clip_and_resample(ds, bounds, global_bounds, TARGET_SIZE)
                    if grid_data is None:
                        continue
                    # print(z, x, y, arr.shape, np.unique(arr))
                    # 构造三角形
                    # 1. 将 NaN 替换为 0 (防止编码报错)
                    grid_data = np.nan_to_num(grid_data, nan=0.0)
                    positions, indices = numpy_grid_to_terrain(grid_data, bounds)

                    out_path_dir = os.path.join(OUT_DIR, str(z), str(x))
                    os.makedirs(out_path_dir, exist_ok=True)
                    file_path = os.path.join(out_path_dir, f'{y}.terrain')

                    with gzip.open(file_path, 'wb') as f:
                        encode(f, positions, indices=indices)

                    tile_count += 1
                    if tile_count % 10 == 0:
                        print(f'  已生成 {tile_count} 个瓦片...', end='\r')
            print(f'Zoom {z} 完成，共生成 {tile_count} 个瓦片')
    # 如果进行了栅格坐标系转换，则移除新建栅格
    if wgs84_path != DEM_FILE:
        os.remove(wgs84_path)
    build_layer_json(OUT_DIR, MIN_ZOOM, MAX_ZOOM, global_bounds, available_list)
    print('all done!')


if __name__ == '__main__':
    main()

Then deploy the generated tiles to Nginx for access by Cesium. The Nginx configuration is as follows:

server {listen 80;server_name localhost;

location /terrain/ {
    alias /usr/share/nginx/html/;
    autoindex on;
    
    # CORS设置
    add_header Access-Control-Allow-Origin * always;
    add_header Access-Control-Allow-Methods "GET, POST, PUT, DELETE, OPTIONS" always;
    
    # layer.json
    location ~* \.json$ {
        add_header Content-Type application/json;
        add_header Access-Control-Allow-Origin * always;
    }
    
    # terrain文件 - 关键修复
    location ~* \.terrain$ {
        add_header Content-Type "application/vnd.quantized-mesh;extensions=octvertexnormals" always;
        add_header Access-Control-Allow-Origin * always;
		add_header Content-Encoding gzip always;
		# gzip on;
    }
}
}

Load this terrain in Cesium. The code is as follows:

// main.js

import './style.css'

import * as Cesium from 'cesium'

import './widgets.css'




window.CESIUM_BASE_URL = '/Cesium/'





/* 2. 创建 Viewer，把地形挂进去 */

const viewer = new Cesium.Viewer('cesiumContainer', {

  homeButton: false,

  sceneModePicker: false,

  baseLayerPicker: false,

  geocoder: false,

  navigationHelpButton: false,

  vrButton: false,

  fullscreenButton: false,

  timeline: false,

  animation: false,

  infoBox: true

})

const terrainProvider = await Cesium.CesiumTerrainProvider.fromUrl('http://localhost:8080/terrain');

viewer.terrainProvider = terrainProvider




// 隐藏 logo

viewer.cesiumWidget.creditContainer.style.display = 'none'




/* 3. 飞到你的数据区域（示例：西安附近） */

viewer.camera.setView({

  destination: Cesium.Cartesian3.fromDegrees(108.86, 35.65, 20000),

  orientation: {

    heading: Cesium.Math.toRadians(0),

    pitch: Cesium.Math.toRadians(-45),

    roll: 0

  }

})

At this point, the browser can load the base map, layer.json, and 0/1/0.terrain, but no matter how much you zoom in, it cannot load the next level of terrain tiles. Why is this?

layer.json

{
  "tilejson": "1.0",
  "name": "terrain",
  "version": "1.0.0",
  "scheme": "tms",
  "tiles": [
    "{z}/{x}/{y}.terrain"
  ],
  "minzoom": 0,
  "maxzoom": 15,
  "bounds": [
    108.75213517708231,
    35.50852907780295,
    108.97661924097244,
    35.801902431140846
  ],
  "projection": "EPSG:4326",
  "format": "quantized-mesh-1.0",
  "available": [
    [
      {
        "startX": 1,
        "startY": 0,
        "endX": 1,
        "endY": 0
      },
      {
        "startX": 3,
        "startY": 1,
        "endX": 3,
        "endY": 1
      },
      {
        "startX": 6,
        "startY": 2,
        "endX": 6,
        "endY": 2
      },
      {
        "startX": 12,
        "startY": 5,
        "endX": 12,
        "endY": 5
      },
      {
        "startX": 25,
        "startY": 11,
        "endX": 25,
        "endY": 11
      },
      {
        "startX": 51,
        "startY": 22,
        "endX": 51,
        "endY": 22
      },
      {
        "startX": 102,
        "startY": 44,
        "endX": 102,
        "endY": 44
      },
      {
        "startX": 205,
        "startY": 89,
        "endX": 205,
        "endY": 89
      },
      {
        "startX": 410,
        "startY": 178,
        "endX": 410,
        "endY": 178
      },
      {
        "startX": 821,
        "startY": 357,
        "endX": 821,
        "endY": 357
      },
      {
        "startX": 1642,
        "startY": 714,
        "endX": 1643,
        "endY": 715
      },
      {
        "startX": 3285,
        "startY": 1428,
        "endX": 3287,
        "endY": 1431
      },
      {
        "startX": 6570,
        "startY": 2856,
        "endX": 6575,
        "endY": 2862
      },
      {
        "startX": 13141,
        "startY": 5712,
        "endX": 13151,
        "endY": 5725
      },
      {
        "startX": 26282,
        "startY": 11424,
        "endX": 26303,
        "endY": 11450
      },
      {
        "startX": 52565,
        "startY": 22848,
        "endX": 52606,
        "endY": 22901
      }
    ]
  ]
}

[Luke_McKinstry]

Hi @wenqin,

Thanks for your post and welcome to the Cesium community. I do not see any glaring mistakes in the layer.json file or CesiumJS code.

A good first step would be to narrow down where the issue is occurring. Are you able to determine if the problem is in the terrain generation, or in how you are serving it to CesiumJS? For testing serving and loading, perhaps you could try load a known working sample of quantized mesh terrain data using your nginx setup? And for testing your terrain generation process, perhaps you could try serving it to CesiumJS using a quanitized mesh terrain server which is known to work? GitHub - heremaps/quantized-mesh-viewer: Render custom quantized mesh tiles in Cesium.js and debug individual tiles using THREE.js renderer. could be a tool to help you do this.

Please let us know if you are able to narrow the problem down further and how we can help next.
Best,
Luke

[wenqin]

Hello, I refined the issue as you suggested and found that the problem still lies with my slice creation. I manually coded the terrain tiles and resolved the issues I encountered. The current result is as follows:

There are two problems：
（1）The basemap only appears in areas where I have terrain. The western hemisphere is blank, as shown.

（2）The basemap disappears when zoomed in to a certain level, specifically around zoom level 8.

（3）When I zoom the map, terrain tiles beyond level 0 are still not being requested. I feel like the configuration in layer.json isn’t taking effect.

I’m not sure how to troubleshoot this next. What do you think the issue might be?