Mapbox GL JS v3 で地理院の標高タイルから色別標高図を作ろうとして上手くいかなかった話

## はじめに

前回の記事で、Mapbox GL JS の新機能について、簡単に紹介してみましたが、今回の Mapbox GL JS v3 では、ラスタタイルの色を変更して表示する機能が実装されています。

https://qiita.com/mg_kudo/items/92283d017f7528f06c5e

公式ドキュメント等はまだ出ていませんが、これを使うことで、国土地理院の標高タイルを使って、標高別に色を塗り分けた地図を作れるのではないかと試行錯誤しましたので、その結果を残します。結論から言うと「**標高タイルの仕様は Mapbox GL JS の思想に合わないのかな**」と思いました。

:::note warn
Mapbox GL JS v3 はまだベータ版であり、公式ドキュメントも出ていないため、この記事をもって結論づけるには時期尚早だと思いますので、その点はご了承ください。
:::


https://github.com/mapbox/mapbox-gl-js

## ラスタタイル関係の新機能
まずは、ラスタタイルの色を変更する関係のスタイルの新設定について、[GitHub を覗いたり](https://github.com/mapbox/mapbox-gl-js/blob/f501fa4498b22f9193a6fe3e8852dc1aeafac812/src/style-spec/reference/v8.json)、使ってみたりして調べた結果をメモします。（公式ドキュメントが出たら、そちらを参照ください……。）
* `raster-color`：この設定に変更したい色を渡します。
* `raster-value`：`raster-color` 内でのみ利用できる設定で、後述の `raster-color-mix` の設定により得られた値を取得できます。
* `raster-color-mix`：RGB 値を一つの値（上述の `raster-value` で取得できるもの）に変換するパラメータを指定します。`[mix.r, mix.g, mix.b, mix.a]` の４パラメータを指定することで、`mix.r * src.r + mix.g * src.g + mix.b * src.b + mix.a` という変換を行うようです。（最後のパラメータは透過度ではなく、固定値のオフセットとのこと。）
* `raster-color-range`：上述の `raster-color-mix` で得られる値の範囲を指定します。この範囲が1024に等分割されて、`raster-color` で利用されるようです。なお、`raster-color` の設定内で加工することで、`raster-color-range` を超えて利用もできるようです。

最後の `raster-color-range` の**1024分割というのが今回の曲者でした**。

## 地理院の標高タイルを表示してみる
まず、[標高タイルの仕様](https://maps.gsi.go.jp/development/demtile.html)を確認してみると、RGB 値から以下のように標高へ変換できると分かります。

* $x = 2^{16}R + 2^8G + B$
* $x < 2^{23}$ の場合　$h = xu$
* $x = 2^{23}$ の場合　$h = NA$
* $x > 2^{23}$ の場合　$h = (x-2^{24})u$

こちらを素直にスタイルに落とし込んでみます。`raster-color` 内では、`let` を用いて、`raster-value` で取り出した値をさらに加工できるようです。以下は、`addLayer()` を用いてスタイルを追加する例となります。 

```JavaScript
// 分解能
const u = 0.01;
// スタイルの追加
map.addLayer({
  "id":"freeRelief",
  "type":"raster",
  "source":"dem",
  "paint":{
    // 範囲設定
    "raster-color-range": [
      0,
      Math.pow(2, 24)
    ],
    // x = 2^16 * R + 2^8 * G + B を計算。
    // この x が "raster-value" で取り出される。
    "raster-color-mix": [
      Math.pow(2, 16) * Math.pow(2, 8),
      Math.pow(2, 8) * Math.pow(2, 8),
      Math.pow(2, 8),
      0
    ],
    // x < 2^23の場合 h = x * u
    // x = 2^23の場合 h = NA (便宜的に 0 とする)
    // x > 2^23の場合 h = (x-2^24) * u
    "raster-color": [
       // mix の値により条件分岐させた結果を
       // 標高値へと計算して "altitude" へ格納。
       "let", "altitude", [
          "case",
          // 2^23 = 8388608 より小さい場合
          ["<", ["raster-value"], Math.pow(2, 23)],
          ["*", ["raster-value"], u],
          // 2^23 = 8388608 より大きい場合
          [">", ["raster-value"], Math.pow(2, 23)],
          ["-", ["*", ["raster-value"], u], Math.pow(2, 24)],
          // NA 値は便宜的に 0 とする
          0
       ],
       [
         "interpolate",
         ["linear"],
         // 上記の "let" で格納した "altitude" を取り出す。
         ["var", "altitude"], 
         // 標高に応じて、グラデーションで色を設定。
         - 10  , "rgb(0,0,255)",
         - 5   , "rgb(0,255,255)",
         + 0   , "rgb(255,255,255)",
         + 5   , "rgb(0,255,0)",
         + 10  , "rgb(255,255,0)",
         + 100 , "rgb(255,155,0)",
         + 500 , "rgb(155,155,0)",
         + 1000, "rgb(155,0,0)",
         + 3777, "rgb(55,0,0)"
       ]
    ],
    "raster-opacity": 0.8
  }
});
```
なお、標高タイルの `source` への追加は以下の通り。
```JavaScript
map.addSource('dem', {
  "type": "raster",
  "minzoom":1,
  "maxzoom":14,
  "tiles":["https://cyberjapandata.gsi.go.jp/xyz/dem_png/{z}/{x}/{y}.png"],
  "tileSize": 256,
  "attribution":"<a href='https://maps.gsi.go.jp/development/ichiran.html#dem' target='_blank'>地理院タイル</a>"
});
```

これでうまくいけばよかったのですが、残念ながらそうはいきませんでした。以下の画像のように、ゼロメートル地帯や標高の低い場所で、グラデーションが思った通りにかかりません。グラデーション設定上は、-10 m から青、0 m で白となり、5 m で緑、10 m で黄色......となるはずなのですが、真っ白のままです。

![試行錯誤その１.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/3e6c3c03-b4b0-9a91-9d88-66930ffb0da0.png)
※背景は国土地理院の最適化ベクトルタイル

いろいろと試してみると、`raster-color-range` の設定が曲者のようです。上述のように、ここで設定された値が 1024 分割されて、`raster-color` に渡されます（つまり、地図へ表示されます）。今回、`raster-color-range` には、0 ～ $2^{24}$ という大きな値がセットされています。$2^{24}$は16,777,216ですから、1024分割すれば、その1 unit は、16,777,216 / 1024 = 16384 となります。つまり、分解能 u = 0.01 を考慮したとしても、色分けの最小単位が約164 m となってしまっているようです。これでは、平野部やゼロメートル地帯等の塗分けはできません。

## 標高タイルとの相性の悪さ

それでは、工夫して対応ができないものでしょうか？

ここで頭を悩ますのが、標高タイルにおける 0 m 前後の計算方法の違いです。標高タイルは、0 m 前後で値が大きく断絶しています。今一度、標高タイルのデコード方法を見てみますが、標高 1 m の時の x は100ですが、標高-1 m の時は、16,777,116となります。つまり、ゼロメートル地帯をカバーする場合、わずか-1 m であっても、非常に大きな数値を `raster-color-range` へ登録する必要があります。

* $x = 2^{16}R + 2^8G + B$
* $x < 2^{23}$ の場合　$h = xu$
  → h=1 の時、$x = h/u = 100$
* $x = 2^{23}$ の場合　$h = NA$
* $x > 2^{23}$ の場合　$h = (x-2^{24})u$
  → h=-1 の時、$x = h/u+2^{24} = 16,777,116$

もともと、Mapbox GL JS v1 の[陰影表現](https://docs.mapbox.com/mapbox-gl-js/style-spec/layers/#hillshade)（`hillshade` レイヤ）や Mapbox GL JS v2 での[地形の3次元表現](https://docs.mapbox.com/mapbox-gl-js/style-spec/terrain/)（`terrain` 機能）の際も、標高タイルへの対応は課題となっておりました。

この課題については、国土地理院も調査研究年報で報告しています。
* [3次元地図データをウェブ地図に表示するための技術的検討（川井他、令和2年度）](https://www.gsi.go.jp/common/000235797.pdf)
※ 図-1に標高タイルと Mapbox GL JS で想定されているデコード方式との違いがまとまっています。

Mapbox GL JS では、`mix.r * src.r + mix.g * src.g + mix.b * src.b + mix.a` のようにデコードできる PNG タイルデータを念頭に開発をしているようで、上記の `hillshade` レイヤも `terrain` 機能も、このようなデコード方式のみがサポートされていました。

妄想ですが、今回の `raster-color` 関連機能も、このデコード方式を前提としつつ、パラメータをスタイル上で設定できるようにしたという流れなのではないかと想像しています。今回は、`raster-color` の設定上、`case` を用いた条件分岐が可能そうでしたので、期待を持っていたのですが、結局 `raster-color-range` の設定に躓いてしまったというところです。

### （参考）先人による対応方法
参考に、国土地理院の標高タイルを Mapbox GL JS 系統へ対応させた事例を列挙してみます。

* 国土地理院は、[地理院地図Vector](https://maps.gsi.go.jp/vector/) において、Mapbox GL JS v1 に[直接手を入れる](https://github.com/gsi-cyberjapan/gsimaps-vector-experiment/blob/master/mapbox/1.13.0/mapbox-gl-gsi.js)ことで陰影表示を標高タイルに対応させたようです（[GitHub](https://github.com/gsi-cyberjapan/gsimaps-vector-experiment/blob/f6d0ffff22670f60cadd5afa85f4bf1d4f6f6866/js/src/map/layer-hillshade2.js#L83C24-L83C24) を見ると、本家にはない `"gsi"` というエンコード方式を設定できるようになっています）。

https://github.com/gsi-cyberjapan/gsimaps-vector-experiment

* 直接手を入れる方法は、以下の記事にも紹介されています。

https://qiita.com/tattii88/items/3d88907c7116d1aeb028

* Mapbox GL JS v1 から派生した MapLibre GL JS の `addProtocol` という機能を用いて、標高タイルをデコードして、色を変換してから表示させるという方法もあるようです。

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

* ServiceWorker を使うという方法もあるようです。以下の記事では、標高タイルへのアクセスを監視し、標高タイルの場合は、色の変換を行うという方法を採用しているようです。

https://qiita.com/frogcat/items/d12bed4e930b83eb3544

* そもそも、Mapbox GL JS に対応したエンコード方式へ変更してしまうのも手です。以下の記事は、webp 形式へ変更した事例です。ただし、標高タイルの 10 m メッシュは、ZL 14 までとはいえ、さすがにデータ量が大きいようです。

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

## 0 m 未満を無視した設定

標高タイルの問題は、0 m を境界に値が大きく異なることが原因のため、0 m 未満（正確には、0 m 未満と NA 値）を無視して、0 m 以上だけで塗り分けることは可能だと考えられます。以下、スタイルに落とし込んでみます。NA 値と0 m 未満は、`raster-color-range` によって上限値となるため、`raster-color` のグラデーション設定の最後に NA 値と0 m 未満用の設定も加えておきます。

```javaScript
// 分解能
const u = 0.01;
map.addLayer({
  "id":"freeRelief",
  "type":"raster",
  "source":"dem",
  "paint":{
    // 範囲を標高（m 単位）で設定。
    "raster-color-range": [
      0,
      4000
    ],
    // x' = 2^16 * R * u + 2^8 * G * u + B * u を計算。
    // この x' が "raster-value" で取り出される。
    "raster-color-mix": [
      Math.pow(2, 16) * Math.pow(2, 8) * u,
      Math.pow(2, 8) * Math.pow(2, 8) * u,
      Math.pow(2, 8) * u,
      0
    ],
    // h = x'（x' < 2^23 * u の場合）のみを考える。
    "raster-color": [
       "interpolate",
       ["linear"],
       ["raster-value"],
       // 標高に応じて、グラデーションで色を設定。
       + 0   , "rgb(255,255,255)",
       + 10   , "rgb(0,255,0)",
       + 50  , "rgb(255,255,0)",
       + 100 , "rgb(255,155,0)",
       + 500 , "rgb(155,155,0)",
       + 1000, "rgb(155,0,0)",
       + 3799, "rgb(55,0,0)",
       // 0 m 以下又は NA 値の色
       + 3800, "rgba(200,200,255,0.5)" 
    ],
    "raster-opacity": 0.8
  }
});
```
一応、それっぽくなったのでは？
![試行錯誤その２.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/d92f0d2e-8a1a-7e13-a27d-668b3d7368d8.png)
※背景は国土地理院の最適化ベクトルタイル

## 感想

今回、Mapbox GL JS v3 で注目していた機能を、やや前のめりに試してみました。今回の記事ではうまくいきませんでしたが、まだ公式ドキュメントがきちんと出ていませんので、他の設定でうまくいく可能性もありますし、`raster-color-range` や `raster-color-mix` 等のパラメータの調整を工夫したりすることでうまくいく可能性もあります。

個人的には、mix のような形ではなく、RGB 値を個別に直接取得・演算して渡せるようになると便利だと思いました。