MapLibre GL JS+標高タイルで3D地形を表示する方法

# はじめに

:::note info
この記事は [MapLibre Advent Calendar 2023](https://qiita.com/advent-calendar/2023/maplibre) 18日目の記事です。
:::

MapLibre GL JSでは地形データがあれば、[3D地形を表示する](https://maplibre.org/maplibre-gl-js/docs/examples/3d-terrain/)ことができます。地形データ（以下、標高タイル）は[DEM（数値標高モデル）](https://fgd.gsi.go.jp/download/ref_dem.html)とも呼ばれ、MapLibre GL JSではDEMを地形データとして利用しますが、[Terrain-RGB形式（Mapbox）（もしくはTerrarium-RGB形式（Mapzen））というエンコーディング（標高換算式）の標高タイルが利用されます](https://maplibre.org/maplibre-style-spec/sources/#raster-dem)。無償で提供されている標高タイルには以下のようなものがあります（観測範囲）。今回は、無償で利用可能な標高タイルをMapLibre GL JSで3D地形として表示する方法について説明します。

| 標高タイルの種類 | エンコーディング <br>（標高換算式） | 仕様 |
| ---- | ---- | ---- |
| 国土地理院 <br> 標高タイル | 独自仕様 | https://maps.gsi.go.jp/development/demtile.html |
| 産業技術総合研究所 <br> シームレス標高タイル | 独自仕様<br>（地理院とほぼ同じ） | https://tiles.gsj.jp/tiles/elev/tiles.html
| AWS(Mapzen) Terrain Tiles | Terrarium-RGB | https://www.mapzen.com/blog/terrain-tile-service/ |

## 国土地理院 標高タイル
国土地理院（以下、地理院）から[標高タイル（DEM）](https://maps.gsi.go.jp/development/ichiran.html#dem)が配信されていますが、今回は、地理院 標高タイルのうち、DEM10B PNG形式の標高タイルを使用して、MapLibre GL JSで3D地形を表示します。なお、地理院 標高タイルは、10mの地上解像度で配信されており、詳細な地形表現が可能です。**ただし、Terrain-RGB形式ではないエンコーディングを利用しているため、Terrain-RGB形式に変換する必要があります**。Terrain-RGB形式と地理院 標高タイルの考え方の違いについては、下記の記事が詳しいです。

https://qiita.com/Kanahiro/items/e22594b738655a189c1d

## 産業技術総合研究所 シームレス標高タイル
産業技術総合研究所（以下、産総研）から[シームレス標高タイル](https://tiles.gsj.jp/tiles/elev/tiles.html)が配信されていますが、今回は、シームレス標高タイルのうち、[統合DEM](https://tiles.gsj.jp/tiles/elev/tiles.html#h_mixed)を使用して、MapLibre GL JSで3D地形を表示します。なお、産総研 シームレス標高タイル（統合DEM）は、地理院 標高タイルの10mDEM（DEM10B）に加え、5mDEM（DEM5A、DEM5B、DEM5C）や[ASTER GDEM 003](https://asterweb.jpl.nasa.gov/gdem.asp)、[GEBCO Grid](https://www.gebco.net/data_and_products/gridded_bathymetry_data/)が統合された形で配信されており、陸域の地形以外にも海底地形の表現が可能です。**ただし、Terrain-RGB形式ではないエンコーディングを利用しているため、Terrain-RGB形式に変換する必要があります**。

## AWS(Mapzen) Terrain Tiles
AWS(Mapzen) Terrain Tilesでは[全球の地形データがAWS S3からタイル形式で配信されています](https://www.mapzen.com/blog/terrain-tile-service/)。AWS(Mapzen) Terrain TilesのエンコーディングはTerrarium-RGB形式ですので、MapLibre GL JSでそのまま3D地形として表示できます。

# 標高タイルのRGB値のエンコーディング（=標高換算式）
- Terrain-RGB形式と標高タイルのRGB値のエンコーディング（=標高換算式）は下記のとおりになります。
（引用元：[『現場のプロがわかりやすく教える 位置情報エンジニア養成講座』](https://www.shuwasystem.co.jp/book/9784798068923.html)）

## Terrain-RGB形式
- RGBピクセル値を3桁の256進数として捉え、横軸が縦軸を標高とした傾き0.1の単調増加関数になります。
```Terrain-RGB形式:Terrain-RGB形式の標高換算式
標高=-10000+((R値×256×256+G値×256+B値)×0.1)
```
- 地理院 標高タイルをTerrain-RGB形式に変換したものです。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205045/d2a7465e-5f82-b22c-c198-a8643e8fd985.png)
**富士山付近**

## 地理院 標高タイル
- 標高タイルの仕様です。

https://maps.gsi.go.jp/development/demtile.html

- 地理院 標高タイルも、Terrain-RGB形式とベースの考え方は全く同じで、計算式が異なり下記のとおりになります。
```標高タイル:標高タイルの標高換算式
標高=(R値×256×256+G値×256+B値)×0.01
# ただし[R,G,B] = [128,0,0]を無効値とする
# ただし[128,0,1]以上の場合2^24を引き去る
```
- 地理院 標高タイルです。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205045/8ee98d0e-bc3c-eee7-b1ad-01b911f66bd1.png)
**富士山付近**

:::note warn
産総研 シームレス標高タイルの仕様は、地理院 標高タイルの仕様とほぼ同一ですが、無効値が異なります。地理院タイルでは(R, G, B)=(128, 0, 0)が無効値となります。
https://tiles.gsj.jp/tiles/elev/tiles.html
:::

## Terrarium-RGB形式
Terrarium-RGB形式の仕様です。

https://docs.safe.com/fme/2017.1/html/FME_Desktop_Documentation/FME_ReadersWriters/terraintilesaws/terraintilesaws_reader.htm

- Terrarium-RGB形式の標高タイルも、Terrain-RGB形式とベースの考え方は全く同じで、計算式が異なり下記のとおりになります。

```Terrarium-RGB形式:Terrarium-RGB形式の標高換算式
標高=(R値×256+G値+B値÷256)-32768
```

- AWS(Mapzen) Terrain Tilesです。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205045/c6262ce2-3dac-d5c5-68cb-c7d7dd06c48a.png)
**富士山付近**

- エンコーディング（標高換算式）を比較すると下図のようになります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205045/72d2b3df-8171-04c3-7251-c60e4e65e9f0.png)
出典）[3 次元地図データをウェブ地図に表示するための技術的検討（国土地理院、令和2年度）](https://www.gsi.go.jp/common/000235797.pdf)

# MapLibre GL JSと地理院 標高タイルで3D地形を表示する
MapLibre GL JSと地理院 標高タイルで3D地形を表示するには、下記の2通りの方法があるようです（観測範囲）。基本的には、どちらの方法でも3D地形を表示することができます。方法1は、すでにホスティングされている、[地理院 標高タイル](https://maps.gsi.go.jp/development/ichiran.html#dem)を利用できるメリットがあり、一方で、変換モジュールによるTerrain-RGB形式への変換が必要というデメリットもあります。方法2は、事前に地理院 標高タイルをTerrain-RGB形式に変換しているので、方法1のような変換モジュールによるTerrain-RGB形式への変換が不要というメリットがあり、一方で、すでにホスティングされている、[地理院 標高タイル](https://maps.gsi.go.jp/development/ichiran.html#dem)が利用できないというデメリットもあります。

- **方法1：地理院 標高タイルを変換モジュールを用いてTerrain-RGB形式に変換して、3D地形として表示する。**

https://qiita.com/Kanahiro/items/1e9c1a4ad6be76b27f0f

- **方法2：事前にTerrain-RGB形式に変換した、地理院 標高タイルを3D地形として表示する。**

https://qiita.com/hfu/items/915d7fb961d05670cab2

## 方法1：地理院 標高タイルを変換モジュールを用いてTerrain-RGB形式に変換して、3D地形として表示する方法
- 標高タイルをTerrain-RGB形式に変換するモジュールの使い方は下記のとおりです。
- 完全なソースコードは下記のGitHubから入手してください。

https://github.com/shiwaku/gsi-terrain-dem-on-maplibre-gl-js-demo

- 標高タイルをTerrain-RGB形式に変換するモジュールです。
:::note info
変換モジュールは下記の記事のソースコードを参考にさせていただきました。
https://qiita.com/Kanahiro/items/1e9c1a4ad6be76b27f0f
:::
```javascript:maplibre-gl-gsi-terrain-fast-png.js
// maplibre-gl-gsi-terrain
// 【参考】https://qiita.com/Kanahiro/items/1e9c1a4ad6be76b27f0f

// 'fast-png'パッケージから'encode'関数をインポート。これは画像データをPNG形式にエンコードするために使用。
import { encode as fastPngEncode } from 'https://cdn.jsdelivr.net/npm/fast-png@6.1.0/+esm';

// RGB値を元に地形の高さを計算し、その高さに対応する新たなRGB値を返す関数
const gsidem2terrainrgb = (r, g, b) => {
    // まず、RGB値を元に地形の高さを計算
    let height = r * 655.36 + g * 2.56 + b * 0.01;

    // 特定のRGB値(128, 0, 0)は高さ0として扱う
    if (r === 128 && g === 0 && b === 0) {
        height = 0;
    } else if (r >= 128) {
        // Rが128以上の場合は、地形の高さから一定値を引く
        height -= 167772.16;
    }

    // 地形の高さに基準値を加算し、さらにスケーリング
    height += 10000;
    height *= 10;

    // 新たなRGB値を計算
    const tB = (height / 256 - Math.floor(height / 256)) * 256;
    const tG =
        (Math.floor(height / 256) / 256 -
            Math.floor(Math.floor(height / 256) / 256)) *
        256;
    const tR =
        (Math.floor(Math.floor(height / 256) / 256) / 256 -
            Math.floor(Math.floor(Math.floor(height / 256) / 256) / 256)) *
        256;

    // 新たなRGB値を返す
    return [tR, tG, tB];
};

// 地形データを扱うためのプロトコルをmaplibreglに追加
maplibregl.addProtocol('gsidem', (params, callback) => {
    // 新しい画像を作成
    const image = new Image();
    image.crossOrigin = '';

    image.onload = () => {
        // キャンバスを作成し、画像のサイズに合わせる
        const canvas = document.createElement('canvas');
        canvas.width = image.width;
        canvas.height = image.height;

        // 2Dコンテキストを取得し、画像を描画
        const context = canvas.getContext('2d');
        context.drawImage(image, 0, 0);

        // 画像のピクセルデータを取得
        const imageData = context.getImageData(
            0,
            0,
            canvas.width,
            canvas.height,
        );

        // すべてのピクセルについて、RGB値を変換
        for (let i = 0; i < imageData.data.length / 4; i++) {
            const tRGB = gsidem2terrainrgb(
                imageData.data[i * 4],
                imageData.data[i * 4 + 1],
                imageData.data[i * 4 + 2],
            );
            imageData.data[i * 4] = tRGB[0];
            imageData.data[i * 4 + 1] = tRGB[1];
            imageData.data[i * 4 + 2] = tRGB[2];
        }

        // fast-pngのencode関数を使用して画像データをPNG形式にエンコード
        const pngData = fastPngEncode({
            width: canvas.width,
            height: canvas.height,
            data: imageData.data,
        });

        // PNGデータをArrayBufferとしてcallback関数に渡す
        callback(null, pngData.buffer, null, null);

    };

    // 画像のURLを取得し、gsidemプロトコル部分を除去してからimage.srcに設定
    image.src = params.url.replace('gsidem://', '');

    // キャンセル処理を返す(今回は特に何もしない)
    return { cancel: () => { } };
});
```

- 標高タイルをTerrain-RGB形式に変換するモジュール`maplibre-gl-gsi-terrain-fast-png.js`を読み込みます。
```html:html
<head>
    <title>国土地理院 標高タイル（DEM10B）</title>
    <meta charset="utf-8" />
    <meta name="viewport" content="width=device-width, initial-scale=1.0">
    <script src='https://unpkg.com/maplibre-gl@3.4.1/dist/maplibre-gl.js'></script>
    <link href='https://unpkg.com/maplibre-gl@3.4.1/dist/maplibre-gl.css' rel='stylesheet' />
    <script src="https://unpkg.com/pmtiles@2.10.0/dist/index.js"></script>
    <script type="module" src="maplibre-gl-gsi-terrain-fast-png.js"></script>
```

- `map.addSource`の`tiles`を下記のとおりにします。
- `gsidem://https://cyberjapandata.gsi.go.jp/xyz/dem_png/{z}/{x}/{y}.png`とすることで、標高タイルがTerrain-RGB形式に変換されます。
- `map.setTerrain({ 'source': 'gsidem', 'exaggeration': 1 });`とすることで、3D地形が表示されます。
- `exaggeration`は標高を強調する倍率になります。
```javascript:javascript
        map.on('load', () => {
            // 標高タイルソース
            map.addSource("gsidem", {
                type: 'raster-dem',
                tiles: [
                    'gsidem://https://cyberjapandata.gsi.go.jp/xyz/dem_png/{z}/{x}/{y}.png',
                ],
                attribution: '<a href="https://maps.gsi.go.jp/development/ichiran.html#dem" target="_blank">地理院タイル(標高タイル)</a>',
                tileSize: 256
            });

            // 標高タイルセット
            map.setTerrain({ 'source': 'gsidem', 'exaggeration': 1 });

        });
```
- デモサイトです。

https://shiwaku.github.io/gsi-terrain-dem-on-maplibre-gl-js-demo/#12.46/35.34967/138.76069/-36/73

## 方法2：事前にTerrain-RGB形式に変換した、国土地理院 標高タイルを3D地形として表示する方法
- Terrain-RGB形式に変換した、国土地理院 標高タイル（DEM10B PNG形式）を下記のGitHubで公開していますので、これを使用しました。

https://github.com/shiwaku/gsi-dem-10b-terrain-rgb

- 完全なソースコードは下記のGitHubから入手してください。

https://github.com/shiwaku/gsi-terrain-dem-terrain-rgb-on-maplibre-gl-js-demo

- `map.addSource`の`tiles`を下記のとおりにします。
- `map.setTerrain({ 'source': 'gsidem-terrain-rgb', 'exaggeration': 1 });`とすることで、3D地形が表示されます。
- `exaggeration`は標高を強調する倍率になります。
```javascript:javascript
        map.on('load', () => {
            // 標高タイルソース
            map.addSource("gsidem-terrain-rgb", {
                type: 'raster-dem',
                tiles: [
                    'https://xs489works.xsrv.jp/raster-tiles/gsi/gsi-dem-terrain-rgb/{z}/{x}/{y}.png',
                ],
                attribution: '<a href="https://maps.gsi.go.jp/development/ichiran.html#dem" target="_blank">地理院タイル(標高タイル)</a>',
                tileSize: 256
            });

            // 標高タイルセット
            map.setTerrain({ 'source': 'gsidem-terrain-rgb', 'exaggeration': 1 });
        });
```
- デモサイトです。

https://shiwaku.github.io/gsi-terrain-dem-terrain-rgb-on-maplibre-gl-js-demo/#12.46/35.34967/138.76069/-36/73

# MapLibre GL JSと産総研 シームレス標高タイルで3D地形を表示する
MapLibre GL JSと産総研 シームレス標高タイルで3D地形を表示するには、Terrain-RGB形式への変換モジュールを使用します。変換モジュールは、産総研 西岡様よりご提供いただきました。ありがとうございます。

変換モジュールは、産総研の下記のリンクから入手することができます。

https://gsj-seamless.jp/fukui/handson/sample/js/numPngProtocol.js

:::note info
変換モジュールは、現在はライセンスが設定されていませんが、オープンソース（Apache License 2.0）として公開していただけるとのことでしたので、それに準じていただければと思います。
:::

- シームレス標高タイルをTerrain-RGB形式に変換するモジュールの使い方は下記のとおりです。
- 完全なソースコードは下記のGitHubから入手してください。

https://github.com/shiwaku/aist-dem-with-relief-on-maplibre-gl-js

- シームレス標高タイルをTerrain-RGB形式に変換するモジュール（numPngProtocol.js）を読み込みます。

```javascript:javascript
    <script type="module">
        // addProtocolを設定
        maplibregl.addProtocol('numpng', makeNumPngProtocol());

        // Terrain-RGB形式への変換モジュール
        import { makeNumPngProtocol } from './numPngProtocol.js'
```

- `map.addSource`の`tiles`を下記のとおりにします。
- `numpng://tiles.gsj.jp/tiles/elev/mixed/{z}/{y}/{x}.png`とすることで、シームレス標高タイルがTerrain-RGB形式に変換されます。
- `map.setTerrain({ 'source': 'gsidem', 'exaggeration': 7 });`とすることで、3D地形が表示されます。
- `exaggeration`は標高を強調する倍率になります。
```javascript:javascript
        map.on('load', () => {
            // 産総研 シームレス標高タイルソース
            map.addSource("aist-dem", {
                type: 'raster-dem',
                tiles: [
                    'numpng://tiles.gsj.jp/tiles/elev/mixed/{z}/{y}/{x}.png',
                ],
                attribution: '<a href="https://tiles.gsj.jp/tiles/elev/tiles.html" target="_blank">産業技術総合研究所 シームレス標高タイル(統合DEM)</a>',
                tileSize: 256
            });

            // 産総研 シームレス標高タイルセット
            map.setTerrain({ 'source': 'aist-dem', 'exaggeration': 7 });
```

- デモサイトです。

https://shiwaku.github.io/aist-dem-with-relief-on-maplibre-gl-js/#6.37/28.57/141.323/-22.4/69

:::note info
シームレス標高タイル（統合DEM）には[海底の地形データ](https://www.gebco.net/data_and_products/gridded_bathymetry_data/)が含まれるので海底地形が3D表示できます。デモサイトでは、地形の起伏をわかりやすくするために、地形の高さを7倍に強調しています。
:::

# MapLibre GL JSとAWS(Mapzen) Terrain Tilesで3D地形を表示する
- AWS(Mapzen) Terrain Tilesは、エンコーディングが[Terrarium-RGB形式](https://docs.safe.com/fme/2017.1/html/FME_Desktop_Documentation/FME_ReadersWriters/terraintilesaws/terraintilesaws_reader.htm)なので、そのままMapLibre GL JSで3D地形として表示することができます。

- 完全なソースコードは下記のGitHubから入手してください。

https://github.com/shiwaku/mapzen-terrain-tiles-on-maplibre-gl-js

- `map.addSource`の`tiles`を下記のとおりにします。
- エンコーディングは、`encoding: "terrarium"`とします。
- `map.setTerrain({ 'source': 'tilezen-dem', 'exaggeration': 1 });`とすることで、3D地形が表示されます。
- `exaggeration`は標高を強調する倍率になります。
```javascript:javascript
        map.on('load', () => {
            // 標高タイルソース
            map.addSource("tilezen-dem", {
                type: 'raster-dem',
                tiles: ['https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png'],
                attribution: '<a href="https://github.com/tilezen/joerd/blob/master/docs/attribution.md">Tilezen Joerd: Attribution</a>',
                encoding: "terrarium"
            });

            // 標高タイルセット
            map.setTerrain({ 'source': 'tilezen-dem', 'exaggeration': 1 });
```
- デモサイトです。

https://shiwaku.github.io/mapzen-terrain-tiles-on-maplibre-gl-js/#12.46/35.34967/138.76069/-36/73

:::note info
AWS(Mapzen) Terrain Tilesのデータソースは下記のGitHubから確認することができます。
https://github.com/tilezen/joerd/blob/master/docs/data-sources.md
:::

# おわりに
無償で利用可能な標高タイルをMapLibre GL JSで3D地形として表示する方法についてご紹介しました。MapLibre GL JSでの3D地形の表示は、もとの標高タイルのエンコーディングがTerrain-RGB形式（もしくはTerrarium-RGB形式）ではない場合は変換モジュールが必要ですが、エンコーディングがTerrain-RGB形式（もしくはTerrarium-RGB形式）の場合は変換モジュールが不要で利用することができます。最後に、MapLibre GL JS+標高タイル（3D地形）を用いたWebマップの事例をいくつかご紹介します。MapLibre GL JSで3D地形が表示できると、様々な地形表現が可能になります。

##### 産総研 シームレス標高タイル（統合DEM）と赤色立体地図・陰陽図（東京都）を表示するデモサイト
https://shiwaku.github.io/aist-dem-with-sekishoku-and-inyou-on-maplibre-gl-js/

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205045/5617ed0e-00bd-bd84-9d77-1aad8ae26a61.png)

:::note info
東京都の赤色立体地図および陰陽図（ともにラスタータイル）は、[open-hinata](https://kenzkenz.xsrv.jp/open-hinata)の開発者である[@kenzkenz](https://twitter.com/kenzkenz)様からご提供いただきました。ありがとうございます。
https://kenzkenz.xsrv.jp/open-hinata/open-hinata.html
:::

##### 産総研 シームレス標高タイル（統合DEM）と地理院タイル（色別標高図、赤色立体地図、陰影起伏図、傾斜量図、空中写真）をMapLibre GL JSで表示するデモサイト

https://shiwaku.github.io/aist-dem-with-gsi-chiriin-tile-on-maplibre-gl-js/

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205045/b067a40b-6b22-ac58-8b71-d1b479cfaeb3.png)

##### 産総研 シームレス標高タイル（統合DEM）とCS立体図・傾斜区分図・レーザ林相図・樹種ポリゴン・森林資源量集計メッシュをMapLibre GL JSで表示するデモサイト

https://shiwaku.github.io/forest-resource-map-on-maplibre-gl-js/

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205045/6eba7172-e82a-2983-ac2d-7c74f9195fb9.png)

:::note info
CS立体図・傾斜区分図・レーザ林相図・樹種ポリゴン・森林資源量集計メッシュは林野庁がG空間情報センターを通じてオープンデータとして公開している森林資源情報を使用しています。
https://front.geospatial.jp/news/2023/10/3575/
:::

##### 地理院 標高タイル、法務省地図XML、3D都市モデルPLATEAUをMapLibre GL JSで表示するデモサイト

https://shiwaku.github.io/2.5d-map-on-mapLibre-gl-js/

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205045/ea4b3543-ebaa-3a7e-c5f8-24e3e79ecd54.png)