この記事は、MIERUNE Advent Calendar 2021の19日目の記事です。
https://kanahiro.github.io/webgl-gis-example/
はじめに
WebGLをGISに応用…といっても、Mapbox GL JS(MapLibre GL JS)やCesium、Deck.glなど、既に多数のWebGISライブラリでWebGLは標準的に活用されています。しかしGPUの能力を活かした並列処理を可能とするWebGLは、GISの世界でまだまだ活用の可能性があると考えています。
ということで、WebGLのGIS応用について色々書いていきたいと思います。今回はDEM編です。次回のネタ・時期は未定です。
DEMとは
いわゆる標高データです。GISでDEMといえば、ピクセル単位に標高値が焼き込まれた画像ファイルのことを言います。標高データがあれば傾斜を計算したり、水域を算出したり、いろいろとGIS的にうれしいことができます。
WebGLの世界にDEMを取り込む
DEMが画像なら、WebGLにはテクスチャとして取り込むのがベストでしょう。GISの世界ではDEMはTIFFファイルであることが多く、ピクセル値には標高値がそのまま保存されています。しかしWebGLはウェブの世界なので、ウェブで使える軽量なフォーマットにする必要があります。今回はPNGとします。
標高値をRGB値にエンコード
TIFF-DEMは見た目だけで考えるといわゆるグレースケール画像です。なのでこれをそのままPNGに変換…としてしまいそうですが、これはナンセンスです、というかもったいない。TIFFのピクセル値は実数値ですがPNGはRGBAそれぞれ8bitの整数値です。もしTIFFをグレースケールPNGに変換したら、PNGが表現できる標高の分解能は256段階になってしまいます（=8bit）。しかしPNGはRGBの3つで8bit、つまり24bitの情報を持てます（16777216段階、アルファもあるので32bitにもできる）。
ここで、GISの世界にはTerrain-RGBという標高値とRGB値の変換ルールが存在します。以下の式で標高値をRGB値に変換してPNGを作成します。
height = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
標高値とRGBの変換については以下の記事が詳しいです。 https://qiita.com/Kanahiro/items/e22594b738655a189c1d#rgb%E5%80%A4%E3%81%AE%E6%A8%99%E9%AB%98%E6%8F%9B%E7%AE%97
ということで手順は割愛しますが、基盤地図情報DEMからTerrain-RGB準拠のPNGファイルを作成しました、富士山周辺です。
エキゾチックな色合いです、なんとなく、これが富士山っぽいなというのは読み取れますが…。このように、TerrainRGBは人間が見て理解できるようには出来ていません。
テクスチャを表示してみる
とりあえずPNGをテクスチャとして取り込んで表示してみましょう。私のお気に入りのluma.glでやります。
npm install @luma.gl/engine @luma.gl/webgl
全コード。
import { AnimationLoop, Model } from '@luma.gl/engine';
import { Buffer, clear, Texture2D } from '@luma.gl/webgl';
const loop = new AnimationLoop({
//@ts-ignore
onInitialize({ gl }) {
const texture = new Texture2D(gl, {
data: 'fujirgb.png',
});
const model = new Model(gl, {
vs: `
attribute vec2 position;
varying vec2 uv;
uniform sampler2D texture;
vec2 getuv(vec2 xy) {
return vec2((xy.x + 1.0) * 0.5, 1.0 - 0.5 * (xy.y + 1.0));
}
void main() {
uv = getuv(position);
gl_Position = vec4(position.xy, 0.0, 1.0);
}
`,
fs: `
uniform sampler2D texture;
varying vec2 uv;
void main() {
vec4 texValue = texture2D(texture, uv);
gl_FragColor = texValue;
}
`,
uniforms: {
texture,
},
attributes: {
position: new Buffer(
gl,
new Float32Array([
-1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, -1,
]),
),
},
vertexCount: 6,
});
return { model };
},
//@ts-ignore
onRender({ gl, model}) {
clear(gl, { color: [0, 0, 0, 1] });
model.draw();
},
});
loop.start({ canvas: 'canvas' });
これでオリジナルのPNG画像と同じようにキャンバスに描画されます。
グレースケール表示してみる
TerrainRGBのカラフル画像だと、標高値の理解が直感的ではありません。富士山っぽいところはなんとなくわかりますが…。ということで、RGB値から標高値を復元し、グレースケール表示してみます。Fragment Shaderのみいじります。
fs: `
uniform sampler2D texture;
varying vec2 uv;
float rgb2height(vec3 rgb) {
return -10000.0 + (rgb.r * 6553.6 + rgb.g * 25.6 + rgb.b * 0.1) * 255.0;
}
void main() {
vec4 texValue = texture2D(texture, uv);
float height = rgb2height(texValue.rgb);
gl_FragColor = vec4(vec3(height / 4000.0), 1.0); 0-4000mでグレースケール化
}
`,
よくあるDEMになりました。TerrainRGBより直感的ですね。
傾斜量図を作ってみる
DEMから傾斜量を求めるアルゴリズムは色々あるみたいですが、今回は雑に、x方向とy方向の傾斜量の平均をとってみたいと思います。またFragment Shaderをいじります。
fs: `
uniform sampler2D texture;
varying vec2 uv;
float rgb2height(vec3 rgb) {
return -10000.0 + (rgb.r * 6553.6 + rgb.g * 25.6 + rgb.b * 0.1) * 255.0;
}
float slope() {
float canvas_px_x = 1.0 / 1024.0;
float canvas_px_y = 1.0 / 1024.0;
//123
//456
//789
//vec4 p1 = texture2D(texture, uv + vec2(-canvas_px_x, -canvas_px_y));
vec4 p2 = texture2D(texture, uv + vec2(0, -canvas_px_y));
//vec4 p3 = texture2D(texture, uv + vec2(canvas_px_x, -canvas_px_y));
vec4 p4 = texture2D(texture, uv + vec2(-canvas_px_x, 0));
//vec4 p5 = texture2D(texture, uv + vec2(0, 0));
vec4 p6 = texture2D(texture, uv + vec2(canvas_px_x, 0));
//vec4 p7 = texture2D(texture, uv + vec2(-canvas_px_x, canvas_px_y));
vec4 p8 = texture2D(texture, uv + vec2(0, canvas_px_y));
//vec4 p9 = texture2D(texture, uv + vec2(canvas_px_x, canvas_px_y));
float h2 = rgb2height(p2.rgb);
float h4 = rgb2height(p4.rgb);
float h6 = rgb2height(p6.rgb);
float h8 = rgb2height(p8.rgb);
float distPerPixel = 38.0;
return ((abs(h6-h4) + abs(h8-h2)) * 0.5) / distPerPixel;
}
void main() {
float slope = slope();
gl_FragColor = vec4(vec3(slope), 1.0);
}
`,
するとこう。傾斜が急なエリアが白くなり、より視認性が増しました。 コード的には、周囲のピクセル値を取得してきて集計するのがポイント。
陰影図
fs: `
uniform sampler2D texture;
varying vec2 uv;
float rgb2height(vec3 rgb) {
return -10000.0 + (rgb.r * 6553.6 + rgb.g * 25.6 + rgb.b * 0.1) * 255.0;
}
float shade() {
float canvas_px = 1.0 / 1024.0;
//123456789
vec4 p1 = texture2D(texture, uv);
vec4 p2 = texture2D(texture, uv + vec2(canvas_px * 1.0));
vec4 p3 = texture2D(texture, uv + vec2(canvas_px * 2.0));
vec4 p4 = texture2D(texture, uv + vec2(canvas_px * 3.0));
vec4 p5 = texture2D(texture, uv + vec2(canvas_px * 4.0));
vec4 p6 = texture2D(texture, uv + vec2(canvas_px * 5.0));
vec4 p7 = texture2D(texture, uv + vec2(canvas_px * 6.0));
vec4 p8 = texture2D(texture, uv + vec2(canvas_px * 7.0));
vec4 p9 = texture2D(texture, uv + vec2(canvas_px * 8.0));
float h1 = rgb2height(p1.rgb);
float h2 = rgb2height(p2.rgb);
float h3 = rgb2height(p3.rgb);
float h4 = rgb2height(p4.rgb);
float h5 = rgb2height(p5.rgb);
float h6 = rgb2height(p6.rgb);
float h7 = rgb2height(p7.rgb);
float h8 = rgb2height(p8.rgb);
float h9 = rgb2height(p9.rgb);
float intensity = 0.0;
if (h1 > h2) {
intensity += 0.125;
};
if (h1 > h3) {
intensity += 0.125;
};
if (h1 > h4) {
intensity += 0.125;
};
if (h1 > h5) {
intensity += 0.125;
};
if (h1 > h6) {
intensity += 0.125;
};
if (h1 > h7) {
intensity += 0.125;
};
if (h1 > h8) {
intensity += 0.125;
};
if (h1 > h9) {
intensity += 0.125;
};
return intensity;
}
void main() {
float shade = shade();
gl_FragColor = vec4(vec3(shade), 1.0);
}
`,
アルゴリズムが雑なので月面みたいな雑な見た目ですが、まあそれなりに陰影図らしくはなっています。もうちょっと頭を使うとすれば、今は右下45度からの光だけを考慮していますが、右方向・下方向の光も加算してやれば、もう少し良い見た目になると思います。
合算してみる
最後にこれまで描画してみた3つの画像を合算してみます。
void main() {
vec4 texValue = texture2D(texture, uv);
float height = rgb2height(texValue.rgb);
float slope = slope();
float shade = shade();
gl_FragColor = vec4(vec3(0.4 * height / 4000.0 + slope * 0.3 + shade * 0.3), 1.0);
}
あんなカラフルな画像からこんな面白い画像が作れました。
今回の全コード
https://github.com/Kanahiro/webgl-gis-example
おわりに
GISの世界では、DEMがあれば色々楽しいことが出来ます。WebGLでも標高値テクスチャがあれば楽しいことがあるとわかりました。そんなわけでWebGLをGISに絡めて開拓中です。お仕事お待ちしております。