はじめに

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

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

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

ラスタタイル関係の新機能

まずは、ラスタタイルの色を変更する関係のスタイルの新設定について、GitHub を覗いたり、使ってみたりして調べた結果をメモします。（公式ドキュメントが出たら、そちらを参照ください……。）

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分割というのが今回の曲者でした。

地理院の標高タイルを表示してみる

まず、標高タイルの仕様を確認してみると、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() を用いてスタイルを追加する例となります。

// 分解能
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 への追加は以下の通り。

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 で黄色......となるはずなのですが、真っ白のままです。

※背景は国土地理院の最適化ベクトルタイル

いろいろと試してみると、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 の陰影表現（hillshade レイヤ）や Mapbox GL JS v2 での地形の3次元表現（terrain 機能）の際も、標高タイルへの対応は課題となっておりました。

この課題については、国土地理院も調査研究年報で報告しています。

3次元地図データをウェブ地図に表示するための技術的検討（川井他、令和2年度）
※ 図-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 において、Mapbox GL JS v1 に直接手を入れることで陰影表示を標高タイルに対応させたようです（GitHub を見ると、本家にはない "gsi" というエンコード方式を設定できるようになっています）。

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

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

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

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

0 m 未満を無視した設定

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

// 分解能
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
  }
});

一応、それっぽくなったのでは？

※背景は国土地理院の最適化ベクトルタイル

感想

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

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