DEM (Digital Elevation Model) を標高以外の任意のデータに準用することを試します。このようなモデルを本投稿では DEM-like と呼ぶことにします。
DEM-like なモデルを使ってラスタ形式のタイル単位でのリクエストを行うことで、標高のように地表に対してべったりと連続して分布した形式のデータを一定のコストで配信できます。
べったりしたデータの例としては、地磁気や重力、日照時間や降水量、気圧などの気象情報、地価などが挙げられます。
一方で、カテゴリルデータのポリゴン（例えば法律で定められた区域など）はベクトルタイルで配信する方が好ましいでしょう。

DEM

DEM は任意の地点における標高を何らかのモデルに変換したものです。DEM への変換は、3次元座標（緯度・経度・標高）の集合（R³）からある DEM のモデルの集合への写像と言い換えられます。

f: R^3 \rightarrow DEM

単に美麗な DEM の PNG タイルを作るだけなら ƒ は好きなように決めて良いですが、データ配信に使う場合は元の標高の値が DEM のモデルに含まれる値から計算できる方が便利です。つまり、ƒ は単射で、かつ逆写像 ƒ^-1 の関数が計算可能である必要があります。

国土地理院の DEM

国土地理院の PNG 標高タイルの元アイデアとなっている文献は以下のようです。

標高を h、標高分解能を 0.01 m、色の各成分の値を R、G、B とすると以下の式になります。h は負でない実数で、R, G, B は 0 以上かつ 255 以下の整数です。

\lceil \frac{h}{0.01} \rceil = 256^2 R + 256 G + B

標高は高々 8,900 m ですから、標高分解能を 0.01 m = 1 cm としても 16,777,216色で十分表現できます。

Mapbox Terrain DEM

Mapbox の DEM は国土地理院のものと計算式が微妙に違っていて、分解能と切片が異なります。これは PNG タイルの色味に影響します。

\lceil \frac{h}{0.1} \rceil = 256^2 R + 256 G + B - 100000

DEM タイルを作ったり比較したりする

以下の画像は同じデータをもとにそれぞれの定義で DEM タイルを作ったものです。左上が標高0m、右下が標高10,000m、左下と右上が標高 5,000m になります。定義上、青成分が繰り返し濃くなったり薄くなったりします。

国土地理院	Mapbox

PNG タイルの生成コマンド

DEM を含むラスタタイルは GDAL を使うと作成できます。

上記の DEM 比較画像の作成方法は以下の通りです。GeoTIFF として得られます。GeoTIFF から PNG タイルにスライスする方法は後述します。

　以下の JavaScirpt で適当に標高の分布データを GeoJSON として作成

import tilebelt from '@mapbox/tilebelt'
const [x, y, z] = [28672, 12800, 15] // とりあえず日本付近のタイル番号で画像を作っているが、場所に特に意味はない
const geojson = {
  type: 'FeatureCollection',
  features: []
}
// 0m - 10000m まで、標高グラデーションを作成
for (let dx = 0; dx < 256; dx++) {
  for (let dy = 0; dy < 256; dy++) {
    const feature = {
      type: 'Feature',
      properties: {
        elevation: Math.ceil(1_000_000 * (((dx + dy) / 2) ** 3) / 256 ** 3) / 100,
        elv_r: 0,
        elv_g: 0,
        elv_b: 0,
      },
      geometry: tilebelt.tileToGeoJSON([x + dx, y + dy, z]),
    }
    geojson.features.push(feature)
  }
}
// 標準出力を `elevation.geojson` として保存
console.log(JSON.stringify(geojson, null, 2))

国土地理院と Mapbox それぞれのモデルに対して標高から RGB をそれぞれ計算

国土地理院と Mapbox のどちらのモデルも標高から RGB の3つの値が一意に定まります。
今回は GDAL の ogrinfo コマンドで GeoJSON プロパティに直接書き込みました。以下のコマンドは国土地理院と Mapbox の場合とで別々に実行してください。

# 国土地理院
$ ogrinfo elevation.geojson -dialect SQLITE -sql "UPDATE elevation SET
  elv_r = (elevation / 0.01) & (POW(256, 3) - 1) >> 16,
  elv_g = (elevation / 0.01) & (POW(256, 2) - 1) >> 8,
  elv_b = (elevation / 0.01) & (POW(256, 1) - 1) >> 0
"

# Mapbox
$ ogrinfo elevation.geojson -dialect SQLITE -sql "UPDATE elevation SET
  elv_r = ((elevation + 10000) / 0.1) & (POW(256, 3) - 1) >> 16,
  elv_g = ((elevation + 10000) / 0.1) & (POW(256, 2) - 1) >> 8,
  elv_b = ((elevation + 10000) / 0.1) & (POW(256, 1) - 1) >> 0
"

GeoTIFF 生成

以下は GeoTIFF を生成するコマンドです。

# 焼き込みの元データを用意。解像度はとりあえず 1000 x 1000
$ gdal_rasterize -co "alpha=yes" -burn 0 -burn 0 -burn 0 -burn 0 -ot Byte -ts 1000 1000 elevation.geojson tmp.tif
# アルファチャンネルは255に
$ gdal_rasterize -b 4 -burn 255 elevation.geojson tmp.tif
# rgb の3バンドに焼き込み
$ gdal_rasterize -b 1 -a elv_r elevation.geojson tmp.tif
$ gdal_rasterize -b 2 -a elv_g elevation.geojson tmp.tif
$ gdal_rasterize -b 3 -a elv_b elevation.geojson tmp.tif

DEM-like / DEM-like タイルを定義する

ある DEM への変換写像が単射であるように、 DEM-like でも単射である必要があります。逆にこの条件が成立するならばどんな値を DEM-like に変換して OK なはずです。DEM-like タイルへの変換写像の定義域は任意のベクタ値の集合で、また、使うかどうかは別としてアルファチャンネルもありますので、終域は4バイトのベクタ値の集合です。

DEM-like タイルのサンプル

複数のデータを詰め込んで DEM-like タイルのサンプルを作ってみます。
国土数値情報から年間降水量、年平均気温、年間日照時間のデータを取得し、これを同じ PNG タイルに詰めます。

データを全てダウンロードして展開し、

#! /bin/bash
mkdir -p tmp
MESH_CODES=(3622 3623 3624 3724 3725 3823 3824 3831 3926 3927 3928 4027 4028 4128 4129 4229 4230 4328 4329 4429 4529 4530 4531 4629 4630 4631 4728 4729 4730 4731 4828 4829 4830 4831 4839 4928 4929 4930 4931 4932 4933 4934 4939 5029 5030 5031 5032 5033 5034 5035 5036 5039 5129 5130 5131 5132 5133 5134 5135 5136 5137 5138 5139 5229 5231 5232 5233 5234 5235 5236 5237 5238 5239 5240 5332 5333 5334 5335 5336 5337 5338 5339 5340 5432 5433 5435 5436 5437 5438 5439 5440 5531 5536 5537 5538 5539 5540 5541 5636 5637 5638 5639 5640 5641 5738 5739 5740 5741 5839 5840 5841 5939 5940 5941 5942 6039 6040 6041 6139 6140 6141 6239 6240 6241 6243 6339 6340 6341 6342 6343 6439 6440 6441 6442 6443 6444 6445 6540 6541 6542 6543 6544 6545 6641 6642 6643 6644 6645 6741 6742 6840 6841 6842)
for MESH_CODE in "${MESH_CODES[@]}"; do
  URL="https://nlftp.mlit.go.jp/ksj/gml/data/G02/G02-12/G02-12_$MESH_CODE-jgd_GML.zip"
  DEST="./tmp/$MESH_CODE.zip"
  curl -sL $URL > $DEST
  unzip $DEST -d ./tmp
  sleep 1
done

シェープファイルをマージ

#! /bin/bash
MERGING=no
for SHAPEFILE in `ls ./tmp/*/*.shp`; do
  if [ $MERGING = no ] ; then
    ogr2ogr ./tmp/data.shp $SHAPEFILE
    MERGING=yes
  else
    ogr2ogr -update -append ./tmp/data.shp $SHAPEFILE -nln data
  fi
done

シェープファイルを更新して各色成分の値を格納

まずカラムを作成。

$ ogrinfo ./tmp/data.shp -sql "ALTER TABLE data ADD r INTEGER"
$ ogrinfo ./tmp/data.shp -sql "ALTER TABLE data ADD g INTEGER"
$ ogrinfo ./tmp/data.shp -sql "ALTER TABLE data ADD b INTEGER"

今回は以下のルールと国土数値情報の仕様に基づき色を決めます。

	赤	緑	青
データ	年間日照時間	年間降水量	年平均気温
シェープファイルのフィールド名	`G02_071`	`G02_014`	`G02_053`
格納されたデータの単位	0.1時間	0.1mm	0.1°C
値の範囲。最大値は想定	1,500時間から2500時間	500mm から 4,000mm	30°C から 10°C（大きいほど減らす）
分解能	5時間	20mm	0.1°C

SQL で RGB の各値を計算します。

$ ogrinfo ./tmp/data.shp -dialect SQLITE -sql "UPDATE data SET
    r = ((G02_071 / 10) - 1500) / 5,
    g = ((G02_014 / 10) - 500) / 20,
    b = (30 - (G02_053 / 10)) / 0.1
"

GeoTIFF 作成

# 焼き込みの元データを用意。解像度はとりあえず 1000 x 1000　
$ gdal_rasterize -co "alpha=yes" -burn 0 -burn 0 -burn 0 -burn 0 -ot Byte -ts 1000 1000 ./tmp/data.shp ./tmp/data.tif
# アルファチャンネルは255に
$ gdal_rasterize -b 4 -burn 255 ./tmp/data.shp ./tmp/data.tif
# rgb の3バンドに焼き込み
$ gdal_rasterize -b 1 -a r ./tmp/data.shp ./tmp/data.tif
$ gdal_rasterize -b 2 -a g ./tmp/data.shp ./tmp/data.tif
$ gdal_rasterize -b 3 -a b ./tmp/data.shp ./tmp/data.tif

これで3バンドに3つのデータを焼き込めました。

それぞれの色成分ごとの値は以下のようになっています。

赤 (年間日照時間)	緑 (年間降水量)	青 (年平均気温)

測地系と座標系追加

# 元のシェープファイルに測地系と座標系の設定がなかったので書き込む
$ gdalwarp ./tmp/data.tif ./tmp/final.tif -t_srs "+proj=longlat +ellps=WGS84"

ラスタタイルに変換
```
# とりあえず ズーム10まで
$ gdal2tiles.py ./tmp/final.tif ./tmp/tiles -z0-10 --xyz
```
これで以下のように地図に上にオーバーレイして表示できます。

値の復元

ラスタタイルから値を取得して、データを取り出してみます。

# Imagemagick を使って [x, y, z] = [898, 409, 10] のタイルの (100, 100) ピクセルの色情報を取り出す
$ magick ./tmp/tiles/10/898/409.png -format "%[hex:p{100,100}]" info:
3687C8FF # rgb(54, 135, 200)

年間日照時間: 54 * 5 + 1500 = 1770 [時間]

年間降水量: 135 * 20 + 500 = 3700 [mm]

年平均気温: 30 - 200 * 0.1 = 10 [°C]

それらしい値がパースできました。

考察

モデルの設計

今回は1つのバンドに対して1つのデータを焼き込みましたが、その辺りの設計は自由で、複数の値を分解してそれぞれのバンドに焼き込んだりすることができるはずです。ただし、データの次数と分解能はトレードオフです。対数スケールを使うのもありです。

4次元のデータを3バンドに対応させて詰め込む例

こういった複雑なデータ変換を想定しているのは可視化の選択肢の数を最大にするためです。データの入れ物なのと共に画像である DEM や DEM-like のタイルには、データ配信と可視化の機能を持たせられます。ただし、可視化と両立した美しいモデル（=良い変換関数）の設計は職人芸的な世界である気がします。

モデルの配信

分解能や切片などのパラメータを metadata.json のプロパティとして配信できると良いかもしれません。クライアントサイドにこれらの係数をハードコードする必要がなくなります。単純な多項式のモデルなら簡単ですが、そうでない場合は何らかのフォーマットで数式を構造化して記述する必要があるかもしれません。

モデルを JavaScript の関数で記述し、toString() してプロパティに含めて配信することもできなくはないですが、セキュリティ的に難があると思います。