標高PNGタイルと WebGL による地形表現

これは [FOSS4G Advent Calendar 2017](https://qiita.com/advent-calendar/2017/foss4g) の記事です。

# はじめに

昨年の [FOSS4G Advent Calendar 2016](http://qiita.com/advent-calendar/2016/foss4g) では [地理院標高タイルと Leaflet でつくるCS立体図](https://qiita.com/frogcat/items/b9aaa3df866a89d46ef5) という記事を書きました。一年たっていろいろと変化があったのでまとめてみます。


# ショーケース

## CS立体図日本版

![image.png](https://qiita-image-store.s3.amazonaws.com/0/82494/be4bd828-6256-5240-529e-f6460cb56593.png)

- <https://frogcat.github.io/csmap-gl/>

[地理院標高PNGタイル](http://maps.gsi.go.jp/development/demtile.html) を使った現行最新版です。PNGタイルをWebGLで直接処理することで、処理性能面の不満がかなり解消されました。いまでは画面上のスライダーに追従して即座に CS 立体図を再合成するような処理ができるようになりました。

## CS立体図世界版 (mapbox)

![image.png](https://qiita-image-store.s3.amazonaws.com/0/82494/9a5f05b2-2c39-66ed-03e3-bf06f564266c.png)

- [https://frogcat.github.io/csmap-gl/mapbox.html](https://frogcat.github.io/csmap-gl/mapbox.html#10/20.8062/203.5918)

CS立体図日本版の標高PNGタイルの取得先を [Mapbox Global elevation data](https://blog.mapbox.com/global-elevation-data-6689f1d0ba65) に変更して少し手を入れたところ、全世界の CS 立体図を表示できるようになりました。


## CS立体図世界版 (mapzen)

![image.png](https://qiita-image-store.s3.amazonaws.com/0/82494/2d0caf8f-e40e-405e-639b-28e36626932a.png)

- [https://frogcat.github.io/csmap-gl/mapzen.html](https://frogcat.github.io/csmap-gl/mapzen.html#8/35.214/139.499)

おなじく [mapzen terrain tile service](https://mapzen.com/documentation/terrain-tiles/) に切り替えて少し手を入れてみたところ、こちらは海底地形まで表示できていますね。


## 陰影段彩等高線図

![image.png](https://qiita-image-store.s3.amazonaws.com/0/82494/d1989711-c998-3335-2a14-a5a9a3242435.png)

- <https://frogcat.github.io/leaflet-tilelayer-glue/>

こちらは CS立体図の試行を応用して、地理院標高PNGタイルだけを持ってきて陰影段彩等高線図を作れるかどうかやってみたものです。できてますね。


# 解説

2016年時点では CSV 形式の標高タイルを Javascript + Canvas でがんばって画像化していたのですが、今は PNG 形式の標高タイルを Web GL で処理することで高速な処理ができるようになりました。

## 標高PNGタイル

もともと国土地理院のCSV形式の標高タイルを使っていたのですが、その後 PNG 形式がリリースされました。Javascript 上でPNGから数値データを抜き取って計算して画像合成するような使いかただと、正直転送コスト削減以上のメリットがないなーと思っていたのですが Web GL の入力として考えるととてもよい形式です。

標高PNGタイルは上述のようにいろいろとプロバイダーがいるのですが微妙にポリシーが違うので注意が必要です。上で触れた「ちょっとした調整」は基本的には RGB→標高変換の式の修正です。

```js
//地理院 (別途条件あり)
var h  = 0.01 * (R * 256 * 256 + G * 256 + B) ;

// mapbox
var h = -10000 + (R * 256 * 256 + G * 256 + B) * 0.1;

// mapzen
var h = (R * 256 + G + B / 256) - 32768 ;

```

## WebGL

WebGL については詳しくなかったのですが、一般的な画像処理（モノクロ化とか輪郭抽出とか陰影図生成）に応用できるということでイチから勉強してみました。初心者なのであまり語れることはないのですが、CS立体図や陰影図を作るのに必要な知識は [WebGL Image Processing](https://webglfundamentals.org/webgl/lessons/webgl-image-processing.html) くらいの内容の理解を目標にすればいいのではないかと思います。

参考までに現時点のCS立体図（地理院標高タイル用）の FragmentShader の定義はこれくらいのボリュームのコードです。


```glsl

precision mediump float;
uniform sampler2D image;
uniform vec2 unit;
uniform vec4 argv;
uniform float zoom;
const vec4 rgb2alt = vec4(256 * 256, 256 , 1, 0) * 256.0 * 0.01;
const mat3 conv_c = mat3(vec3(0,-1, 0),vec3(-1, 4,-1), vec3(0,-1, 0));
const mat3 conv_sx = mat3(vec3(-1, 0, 1),vec3(-2, 0, 2),vec3(-1, 0, 1));
const mat3 conv_sy = mat3(vec3(-1,-2,-1),vec3(0, 0, 0),vec3( 1, 2, 1));
const vec3 color_convex  = vec3(1.0,0.5,0.0);
const vec3 color_concave = vec3(0.0,0.0,0.5);
const vec3 color_flat    = vec3(0.0,0.0,0.0);
float conv(mat3 a, mat3 b){
  return dot(a[0],b[0]) + dot(a[1],b[1]) + dot(a[2],b[2]);
}
float alt(sampler2D i,vec2 p){
  return dot(texture2D(i, p), rgb2alt);
}
void main() {
  vec2 p = vec2(gl_FragCoord.x,1.0 / unit.y - gl_FragCoord.y);
  mat3 h;
  h[0][0] = alt(image, (p + vec2(-1,-1)) * unit);
  h[0][1] = alt(image, (p + vec2( 0,-1)) * unit);
  h[0][2] = alt(image, (p + vec2( 1,-1)) * unit);
  h[1][0] = alt(image, (p + vec2(-1, 0)) * unit);
  h[1][1] = alt(image, (p + vec2( 0, 0)) * unit);
  h[1][2] = alt(image, (p + vec2( 1, 0)) * unit);
  h[2][0] = alt(image, (p + vec2(-1, 1)) * unit);
  h[2][1] = alt(image, (p + vec2( 0, 1)) * unit);
  h[2][2] = alt(image, (p + vec2( 1, 1)) * unit);
  float z = 10.0 * exp2(14.0 - zoom);
  vec2 cs = h[1][1] > 4000.0 ? vec2(0) : clamp(vec2(
    conv(h,conv_c),
    length(vec2(conv(h , conv_sx),conv(h , conv_sy)))
  ) * vec2(argv[0] / z,argv[1] / z), -1.0 ,1.0);
  gl_FragColor = vec4(
    cs[0] > 0.0 ? mix(color_flat,color_convex,cs[0]) : mix(color_flat,color_concave,-cs[0]) ,
    cs[1]
  );
}
```


## Leaflet

上記のデモ群は Leaflet 上で動いています。
Leaflet で PNG タイルを Web GL で加工して出力するプラグインとして [Leaflet.TileLayer.GL](https://gitlab.com/IvanSanchez/Leaflet.TileLayer.GL) があるのですが、こちらはひとつのPNGタイルごとにひとつの Canvas を作って描画する仕組みなので、段彩図だけならばまだしも、陰影図のような処理ではタイル境界部分に難があります。これに対処するために、複数の PNG をまとめたものをひとつの Canvas に描画するプラグイン [leaflet-tilelayer-glue](https://github.com/frogcat/leaflet-tilelayer-glue) というものを作っています。

今後の展開を考えると（たとえばベクトルタイルの背景図に使ったり）、応用先は Leaflet 一辺倒ということもないでしょうから、Leaflet 依存性がない形でコードを維持していくのが重要ですね。ただ、処理のコアは「画像をつなげたものにWebGLでエフェクトかけて出力」というシンプルなものなので案外移植はしやすいのではと思っています。

## カラーリング

CS立体図の着色、最初期はこんなかんじで塗っていたのですが

![image.png](https://qiita-image-store.s3.amazonaws.com/0/82494/92893343-5938-b6a3-aaa1-a746941cd025.png)


一年を経てこんなかんじの塗り方に落ち着きました。


![image.png](https://qiita-image-store.s3.amazonaws.com/0/82494/f9e7bccd-70b9-4829-70fd-091944ecb190.png)

北海道・倶多楽湖周辺 : <https://frogcat.github.io/csmap-gl/#13/42.4934/141.1930>


どうも RGB だけのベタ塗りだとどこかグロテスクというか禍々しい印象の絵になりがちで、ほかの地図とブレンドして使用しようとしても調和させにくいなー、と思っていました。このへんは CS立体図 の関係者各位試行錯誤されているところかと思いますが、最終的にはこんなカラーリングポリシーに落ち着きました。

|         |最初期カラーリング        |現在のカラーリング        |
|---------|----------|----------|
|傾斜     |平坦(白) <br/>緩斜面(ピンク)<br/>急斜面(赤) |平坦(100%透過)<br/>緩斜面(50%透過) <br/>急斜面(0%透過)|
|曲率     |凸面(白) <br/> 平坦(水色) <br/> 凹面(青)|凸面(オレンジ) <br/> 平坦(黒) <br/> 凹面(ネイビー)        |
|透過     |レイヤー全体の透過度を変更        | 特に変更しない（傾斜によってピクセルごとに透過）        |


- RGBA でカラーリングを考えることで、背景地図を選択的・段階的に透過
- オレンジ、ネイビーを使ったカラーユニバーサルデザイン的にやさしい配色


# まとめ

紆余曲折あって、この一年で標高PNGタイルと WebGL を組み合わせることで、ブラウザ上で実用的な速度で地形表現を合成することができるようになりました。
去年も締めで書いたのですが、自治体が高精度の標高タイルをオープンデータとして出す、みたいな事例ができるといいなーと相変わらず思っています。特に PNG 形式はブラウザやモバイル環境での利活用につながるのではないでしょうか。