Mapterhornの地形データをCesiumで表示する

Mapterhornというめちゃくちゃ良いオープンデータプロジェクトがあります。全球の地形タイルを無料で配信してくれていて、しかもZXY形式の他にPMTiles形式もあります。これを使えばMapLibreなんかで3D地形がサクッと表示できます。

以下のようなビューワーも用意されています。

ところが、Cesiumで使おうとすると話が変わってきます。全然サクッといかず、結構ハマります。

この記事では、Mapterhornの地形データをCesium上で動かすまでに踏んだ技術的なステップを書いていきます。データ形式の変換、カスタムデコーダーの差し込み、Web Workerによる並列化まで順を追って説明していきます！

Mapterhornとは

Mapterhornは、オープンデータのDEM（数値標高モデル）を集約して、Webマップ向けの地形タイルとして配信するプロジェクトです。Oliver Wipfli氏がNLnetの支援を受けて開発しています。

弊社MIERUNEもサポートしています！
https://www.mierune.co.jp/news/0p_mmrdem

特徴はこんな感じです。

• Terrarium形式のRGBエンコーディングでWebP画像として配信
• 512×512ピクセルのタイルサイズ
• PMTilesアーカイブ or https://tiles.mapterhorn.com/{z}/{x}/{y}.webp のZXYエンドポイント
• 全球30m解像度、欧州を中心に0.4m〜20mの高解像度データ
• MapLibre GLではネイティブに encoding: "terrarium" を指定するだけで動く

MapLibreユーザーはとても便利に利用できますが、Cesiumユーザーにとってはそうではないです。

なぜCesiumにそのまま読み込めないのか

Cesiumが地形データとして受け付けるのは、基本的にquantized-meshというCesium独自のバイナリ形式です。

quantized-meshは、大規模な地形データをWebでストリーミング配信するために設計されたフォーマットで、TIN（不整三角形網）をタイル単位で効率的にエンコードします。具体的には、各タイルが以下のような構造を持っています。

• ヘッダ: タイルの境界矩形、最小/最大標高などのメタデータ
• 頂点データ: 経度・緯度・高さの3属性を、それぞれ0〜32767の16bit整数に量子化して格納し、差分値をzigzagエンコーディングで圧縮
• インデックス: 三角形メッシュの頂点インデックス配列
• 辺の頂点リスト: 隣接タイルとの接合に使う東西南北の辺上の頂点インデックス
• 拡張データ（オプション）: 法線ベクトル、水面マスク、メタデータJSON等

つまり、ピクセルの集合であるラスター画像とはまるで別物のデータ構造です。Mapterhornが配信するTerrarium形式のWebPタイルを、このquantized-meshに変換する必要があります。

Mapterhornのタイルが直接使えない理由は以下のとおりです。

1. エンコーディングの違い: MapterhornはTerrarium形式です。cesium-martiniが標準で対応しているTerrain-RGB形式とはデコード式が異なります
2. 座標系の違い: タイルはWeb Mercator（EPSG:3857）で配信されますが、CesiumはEPSG:4978（ECEF: 地球中心の直交座標系）で動きます
3. 標高 vs 楕円体高: タイルに入っているのはジオイド基準の標高であって、Cesiumが内部で使う楕円体高ではありません。ただ、今回の実装ではジオイド補正は行わず、標高値をそのまま楕円体高として扱っています
4. フォーマット変換: PNG/WebPのラスターピクセルから、三角形メッシュ（quantized-mesh）への変換が必要です

つまり、やるべきことはこうです。

WebP画像（Terrarium encoding, EPSG:3857タイル）
  → RGBデコードで標高値を抽出
  → MARTINIで適応的三角形メッシュを生成
  → quantized-meshにパッキング
  → タイルのZXYから地理的範囲を算出し、quantized-meshと紐づける
  → Cesiumがレンダリング時にECEF座標へ変換

結構大変ですねー。楽しそうです。

cesium-martiniライブラリを利用する

ここで登場するのがcesium-martiniです。

このライブラリは、内部でMapboxのMARTINIアルゴリズムを利用しています。
MARTINIは、標高値の正方グリッド（Float32Array）からRTIN（Right-Triangulated Irregular Network）による適応的な三角形メッシュを生成するアルゴリズムです。平坦な場所は大きな三角形で、急峻な場所は細かい三角形で表現するため、頂点数を大幅に削減できます。cesium-martiniはこのMARTINIの出力をquantized-mesh形式にパックし、CesiumのTerrainProviderとして差し込めるようにしてくれます。

ただし、デフォルトではMapbox Terrain-RGB形式しか読めません。

Terrain-RGBのデコード式:

height = -10000 + (R × 65536 + G × 256 + B) × 0.1

Terrariumのデコード式:

height = R × 256 + G + B / 256 - 32768

見ての通り、全然違います。エンコード体系が根本的に異なるので、そのまま使うと標高値がめちゃくちゃになります。

カスタムデコーダーを差し込む

ありがたいことに、cesium-martiniはプロバイダー側でカスタムデコーダーを受け取れる設計になっています。TerrainDecoderインターフェースを実装すれば、デコード処理を丸ごと差し替えられます。

Terrarium用のデコーダーはこう書けます。

// terrarium-worker.ts

/** Terrarium encoding: height = r * 256 + g + b / 256 - 32768 */
function terrariumDecode(r: number, g: number, b: number): number {
    return r * 256 + g + b / 256 - 32768;
}

Rチャネルが256m単位、Gが1m単位、Bが1/256m単位（約3.9mm精度）を表現していて、-32768のオフセットで深海の海底地形もカバーします。

この関数をrgbTerrainToGridに渡すと、512×512のRGBAピクセルデータから513×513のFloat32Array地形グリッドが生成されます。
1ピクセル増えるのは、MARTINIアルゴリズムがグリッドサイズ 2^n + 1 を要求するためです。足りない右端の列と下端の行は、隣接する値を複製して埋めています。

import { rgbTerrainToGrid, createQuantizedMeshData } from '@macrostrat/cesium-martini';
import ndarray from 'ndarray';
import Martini from '@mapbox/martini';

const martiniCache: Record<number, Martini> = {};

function buildQuantizedMesh(parameters: TerrainWorkerInput, transferableObjects: Transferable[]) {
    const { imageData, tileSize = 512, errorLevel, maxVertexDistance } = parameters;

    // RGBAピクセルをndarrayで構造化
    const pixels = ndarray(
        new Uint8Array(imageData),
        [tileSize, tileSize, 4],   // 512×512ピクセル、各4チャネル
        [4, 4 * tileSize, 1],      // ストライド: ピクセル間隔、行間隔、チャネル間隔
        0
    );

    // Terrariumデコーダーで標高グリッドに変換
    const terrain = rgbTerrainToGrid(pixels, terrariumDecode);

    // MARTINIインスタンスをキャッシュ（同じtileSizeなら使い回す）
    martiniCache[tileSize] ??= new Martini(tileSize + 1);

    const tile = martiniCache[tileSize].createTile(terrain);
    const mesh = tile.getMesh(errorLevel, Math.min(maxVertexDistance, tileSize));

    // Cesium用のquantized-meshデータを生成
    const res = createQuantizedMeshData(tile, mesh, tileSize, terrain);
    transferableObjects.push(res.indices.buffer);
    transferableObjects.push(res.quantizedVertices.buffer);
    if (res.quantizedHeights) {
        transferableObjects.push(res.quantizedHeights.buffer);
    }
    return res;
}

ここで起きていることを整理します。

1. WebPからピクセルデータを取得: CanvasのgetImageData()でRGBA配列にします
2. ndarrayでストライド指定: メモリ上のRGBAデータを[row, col, channel]の3次元配列として扱います
3. Terrariumデコード: 各ピクセルのRGB値から標高（メートル）を計算します
4. MARTINIメッシュ生成: 標高値のグリッドから適応的な三角形メッシュを生成します
5. quantized-mesh変換: createQuantizedMeshDataがグリッド座標（0〜512）を0〜32768に正規化し、タイル内の相対位置として量子化します

ここまでがWorker内の処理です。Worker内では頂点のタイル内相対位置（0〜32768）を計算するだけで、地球上のどこに位置するかはまだ決まっていません。

Workerから結果が返ると、cesium-martiniのメインスレッド側（terrain-provider.tsのprocessTile）で以下の処理が行われます。

1. WebMercatorTilingScheme.tileXYToRectangle(x, y, z)がタイルのZXYから地理的範囲（tileRect、ラジアン）を算出
2. そのtileRectからOrientedBoundingBox（タイルの3D境界ボックス）を計算
3. Workerから返ってきたquantized-meshの頂点データと、このtileRect・境界ボックスを一緒にQuantizedMeshTerrainData（Cesiumのクラス）に渡す

ここでやっているのは「このタイルが地球上のどこにあるか」という位置情報の付与であって、頂点座標そのものの変換ではありません。実際の座標変換はCesiumの内部で行われます。
Cesiumは、量子化された0〜32768の頂点座標をtileRectの範囲にマッピングして地理座標に展開し、レンダリング時にECEF（EPSG:4978）へ変換します。

Web Workerで並列化

さて、デコーダーの差し替えはできました。しかしもう一つ問題があります。

cesium-martiniのビルトインWorkerFarmTerrainDecoderも、重い計算はWorkerに委譲しています。ただしWorkerを1つしか使わないようでした。タスクはキューイングされて逐次処理されるので、地形タイルが大量にリクエスされると明らかにボトルネックになります。

cesium-martiniはTerrainDecoderというインターフェースでデコード処理を抽象化しています（Cesium本体のAPIではなく、cesium-martini独自の設計）。Cesiumがタイルを1枚必要とするたびに、MartiniTerrainProviderは以下の2つのメソッドを順に呼びます。

1. requestTileGeometry: 同時リクエスト数を制御する。上限を超えていればundefinedを返して拒否し、Cesium側にリトライさせる。通過すればタイル画像の取得とピクセル抽出を行う
2. decodeTerrain: ピクセルデータを受け取り、実際の重い計算（デコード → メッシュ生成 → quantized-meshパッキング）を行う

このインターフェースを実装すれば、Workerの数やルーティング方法を自由に決められます。そこでWorkerを複数起動して負荷分散するWorker Poolを自作しました。1枚のタイルが処理される流れはこうなります。

Cesium: タイルが必要
  │
  ▼
MartiniTerrainProvider.requestTileGeometry(x, y, z)
  │
  ▼  decoder.requestTileGeometry({x,y,z}, processTile)
PooledTerrainDecoder (自作 / メインスレッド)
  │  inProgress > maxRequests なら undefined を返して拒否
  │  そうでなければ processTile を実行
  │    → タイル画像をfetch、ピクセルデータを抽出
  │
  ▼  decoder.decodeTerrain(params, pixelData.buffer)
  │  pendingが最も少ないWorkerHandleを選択
  │
  ▼  workerHandle.send({ id: 42, payload }) → postMessage
WorkerHandle (自作 / メインスレッド)
  │  メッセージIDを振ってPromiseを保持
  │
  ▼  ---- スレッド境界 (Transferableでゼロコピー) ----
terrarium-worker.ts (自作 / Workerスレッド)
  │  self.onmessage で受信
  │  buildQuantizedMesh() を実行
  │
  ▼  self.postMessage({ id: 42, payload: result })
WorkerHandle (メインスレッド)
  │  onmessage でID=42のPromiseをresolve
  │
  ▼  QuantizedMeshResult が返る
MartiniTerrainProvider
      tileRectと合わせてCesiumに渡す

ここに登場する自作のクラスは3つです。順に見ていきます。

terrarium-worker.ts

Worker側で動くプログラムです。メインスレッドからpostMessageで送られてきたピクセルデータを受け取り、Terrariumデコード → MARTINIメッシュ生成 → quantized-meshパッキングまでを実行して、結果を返します。

// terrarium-worker.ts（Worker側）

self.onmessage = function (msg) {
    const { id, payload } = msg.data;
    if (id == null) return;
    const objects: Transferable[] = [];
    const res = buildQuantizedMesh(payload, objects);
    // 同じIDを付けて返す → メインスレッド側でPromiseが解決される
    self.postMessage({ id, payload: res }, objects);
};

buildQuantizedMeshの中身は前節で説明したとおりです。ポイントは、メインスレッドから受け取ったidをそのまま返していることです。メインスレッド側はこのIDでどのリクエストへの応答かを特定します。

WorkerHandle

メインスレッド側で、個々のWorkerインスタンスを管理するクラスです。

// pooled-decoder.ts（メインスレッド側）

class WorkerHandle {
    worker: Worker;
    pending: number = 0;  // このWorkerの処理中タスク数
    private resolves: Record<number, (value: unknown) => void> = {};
    private rejects: Record<number, (reason: unknown) => void> = {};

    constructor(worker: Worker) {
        this.worker = worker;
        worker.onmessage = (e) => {
            const { id, payload } = e.data;
            this.pending--;
            this.resolves[id]?.(payload);
            delete this.resolves[id];
            delete this.rejects[id];
        };
    }

    send(payload: unknown, transferables: Transferable[]): Promise<unknown> {
        const id = globalMsgId++;
        this.pending++;
        return new Promise((resolve, reject) => {
            this.resolves[id] = resolve;
            this.rejects[id] = reject;
            this.worker.postMessage({ id, payload }, transferables);
        });
    }
}

やっていることはシンプルです。sendでメッセージIDを振ってPromiseを作り、Workerのonmessageで同じIDのPromiseを解決する。これでpostMessage/onmessageという低レベルなAPIを、await workerHandle.send(data)というPromiseベースの呼び出しに変換しています。pendingは現在処理中のタスク数で、次に説明するPool側の負荷分散に使います。

PooledTerrainDecoder

複数のWorkerHandleを束ねて、cesium-martiniのTerrainDecoderインターフェースを実装するクラスです。デフォルトのWorkerFarmTerrainDecoderとの違いは、Workerの数と振り分け方です。

export class PooledTerrainDecoder implements TerrainDecoder {
    private workers: WorkerHandle[];
    private inProgress = 0;
    private maxRequests: number;

    constructor(workerFactory: () => Worker, poolSize: number, maxRequests?: number) {
        this.workers = Array.from({ length: poolSize }, () =>
            new WorkerHandle(workerFactory())
        );
        this.maxRequests = maxRequests ?? poolSize * 2;
    }

    async decodeTerrain(
        params: TerrainWorkerInput,
        data: ArrayBufferLike
    ): Promise<QuantizedMeshResult> {
        // 最もpendingが少ないWorkerを選択
        const worker = this.workers.reduce((a, b) =>
            (a.pending <= b.pending ? a : b)
        );
        return worker.send(params, [data]) as Promise<QuantizedMeshResult>;
    }
}

コンストラクタの第1引数はWorkerのファクトリ関数（() => Worker）です。Pool自体はWorkerの中身に関知せず、この関数をpoolSize回呼んでインスタンスを生成するだけです。Cesiumがタイルを要求するたびにdecodeTerrainが呼ばれ、pendingが最も少ないWorkerを選んで仕事を振ります。

その他の設計判断は以下のとおり。

• プールサイズ: navigator.hardwareConcurrencyを参照して、CPUコア数に応じて1〜8個のWorkerを起動します
• バックプレッシャー: maxRequestsを超えるリクエストはundefinedを返して拒否します。Cesium側がリトライしてくれるようです
• Transferable: ArrayBufferの所有権をWorkerに移転するので、メインスレッドとWorker間でデータのコピーが発生しません

全体をつなげる

最終的な組み立てはシンプルです。

// create-terrain-provider.ts
import { MartiniTerrainProvider, DefaultHeightmapResource } from '@macrostrat/cesium-martini';
import { PooledTerrainDecoder } from './pooled-decoder';
import TerrariumWorker from './terrarium-worker?worker';

const WORKER_POOL_SIZE = Math.min(navigator.hardwareConcurrency ?? 4, 8);

export function createTerrainProvider() {
    const resource = new DefaultHeightmapResource({
        url: 'https://tiles.mapterhorn.com/{z}/{x}/{y}.webp',
        maxZoom: 17,
        tileSize: 512
    });

    const decoder = new PooledTerrainDecoder(
        () => new TerrariumWorker(),    // Viteの?workerインポートでバンドル
        WORKER_POOL_SIZE,
        WORKER_POOL_SIZE * 3
    );

    return new MartiniTerrainProvider({
        resource,
        decoder,
        detailScalar: 1.0,         // メッシュの細かさ
        minimumErrorLevel: 0.1,    // 最小誤差0.1m
        fillPoles: true            // 極地を補完
    });
}

MartiniTerrainProviderに対して、標準のDefaultHeightmapResource（タイルの取得を担当）と、自作のPooledTerrainDecoder（Terrariumデコード + quantized-mesh変換を担当）を渡すだけです。

Viteの?workerサフィックスのおかげで、Worker用のTypeScriptファイルが自動的にバンドル・インライン化されます。ビルド設定で特別なことをする必要はありません。

そしてCesiumのViewerに接続します。

// create-viewer.ts
const viewer = new Viewer(container, {
    baseLayer: new ImageryLayer(
        new UrlTemplateImageryProvider({
            url: 'https://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/{z}/{x}/{y}.jpg',
            credit: '国土地理院'
        })
    ),
    requestRenderMode: true,
    // ...その他UI設定
});

viewer.terrainProvider = terrainProvider;
viewer.scene.globe.depthTestAgainstTerrain = true;

// 富士山を見る
viewer.camera.setView({
    destination: Cartesian3.fromDegrees(138.7274, 35.3606, 15000),
    orientation: {
        heading: CMath.toRadians(0),
        pitch: CMath.toRadians(-30),
        roll: 0
    }
});

地理院の航空写真タイルをベースレイヤーにして、Mapterhornの地形データを重ねます。富士山が3Dで眼前にどーんと現れます。いいですねー。

データ変換パイプラインの全体像

改めて、タイルがブラウザに届いてから3Dメッシュになるまでの全工程を整理します。

[Mapterhornサーバー]
  WebPタイル (512×512, Terrarium encoding, EPSG:3857)
        │
        ▼
[メインスレッド / processTile]
  fetch → ImageBitmap → Canvas.drawImage → ImageData (Uint8ClampedArray)
        │
        ▼ Transferable (ゼロコピー)
[Web Worker (×1〜8)]
  (1) ndarray化: RGBA配列を [row, col, channel] の3次元構造に
  (2) Terrariumデコード: R×256 + G + B/256 - 32768 → Float32 標高値
  (3) 地形グリッド生成: 512×512px → 513×513の標高グリッド
  (4) MARTINIメッシュ生成: 適応的三角形分割
  (5) quantized-mesh生成: 頂点をUint16に量子化、インデックス配列を構築
        │
        ▼ Transferable (ゼロコピー)
[メインスレッド / processTile]
  (6) tileXYToRectangle: タイルのZXYから地理的範囲 (tileRect) を算出
  (7) tileRectからOrientedBoundingBox（3D境界ボックス）を計算
  (8) QuantizedMeshResult + tileRect + OBB → QuantizedMeshTerrainData を生成
        │
        ▼
[Cesium]
  量子化された頂点座標をtileRect内の地理座標にマッピング → ECEF変換 → WebGLレンダリング

重い処理（ステップ1〜5）はすべてWorker上で実行されます。メインスレッドが担うのはタイルの取得（fetch）と、Workerから返ってきた結果への座標情報の付与（ステップ6〜8）です。

まとめ

Mapterhornは本当にありがたいプロジェクトです。高品質な地形タイルがオープンデータとして誰でも使えます。ただ、MapLibre GLでは簡単に3D地形が出ますが、Cesiumではそうはいきません。それなりの変換パイプラインを自前で組む必要があります。

ただ、cesium-martiniのカスタムデコーダー機構が思いのほか柔軟で、Terrarium形式のデコーダーをきれいに差し込めました。さらに、デフォルトの単一Worker構成をPooled Workerに置き換えたことで、パフォーマンスもいい感じになりました。

皆さんもCesiumでMapterhornの地形タイルを動かしてみてください！