JavaScriptでDEM解析Webを作る ー 2D・3D表示と大容量GeoTIFF対応
はじめに
DEMをWebブラウザ上で表示し、そのまま簡単な地形解析まで行える環境を作ってみます。
最初のプロトタイプでは、GeoTIFF形式のDEMを読み込み、OpenLayersによる2D表示とThree.jsによる3D表示を同じ画面で扱えるようにします。さらに、標高確認、傾斜量、陰影起伏までをJavaScript側で処理します。
対象とする機能は、
機能一覧
- 2D地図表示
- 3D地形表示
- クリック地点の標高取得
- 標高段彩
- 傾斜量
- 陰影起伏
- 大容量GeoTIFFのLOD表示
- Web Workerによる分割解析
までを、JavaScript中心でまとめます。
バックエンド側にFastAPIやGDALサーバーは置かず、Docker Composeでビルドした静的WebアプリをNginxから配信する構成です。
このプロトタイプのプロジェクト名は次のとおりです。
dem_analysis_web_r002
参考のために、ZIPファイルを置いておきますので、ダウンロードして活用してください。
このプロトタイプで実装する機能
主な機能は次のとおりです。
| 機能 | 内容 |
|---|---|
| GeoTIFF読込 | ローカルのGeoTIFFをブラウザから選択 |
| COG対応 | Cloud Optimized GeoTIFFも利用可能 |
| 2D表示 | OpenLayers上へDEMを表示 |
| 3D表示 | Three.jsでDEMを地形メッシュ化 |
| DEM全体表示 | 読み込んだGeoTIFFの範囲へ自動フィット |
| 標高段彩 | 標高値に応じた色分け |
| クリック標高 | 地図上をクリックして標高を取得 |
| 傾斜量 | Horn法による傾斜角計算 |
| 陰影起伏 | DEMからHillshadeを生成 |
| LOD表示 | 大きなDEMは縮小表示用データを使用 |
| 全セル分割解析 | GeoTIFFをwindow単位に分割して解析 |
| 高速LOD解析 | 縮小DEM上で傾斜・陰影を計算 |
| 解析進捗 | Web Workerの処理状況を表示 |
| 解析キャンセル | 長時間処理を途中で停止 |
| 3D標高強調 | 地形の高さ倍率を変更 |
| 3Dワイヤーフレーム | 地形メッシュ構造を確認 |
| 座標系対応 | EPSG:3857、EPSG:4326、JGD2011平面直角座標系 |
全体構成
今回の構成は比較的シンプルです。
Docker Composeでは、Node.js + Viteでフロントエンドをビルドし、実行時はNginxから静的ファイルを配信します。
Docker Compose
│
└─ web
├─ build
│ └─ Node.js + Vite
│
└─ runtime
└─ Nginx
バックエンドAPIを置いていないため、構成自体は軽量です。
使用するライブラリ
主に次を使います。
OpenLayers
GeoTIFF.js
Three.js
Proj4js
Vite
Web Worker
Nginx
Docker Compose
役割を分けると次のようになります。
| ライブラリ | 用途 |
|---|---|
| OpenLayers | 2D地図表示 |
| GeoTIFF.js | GeoTIFF読込、window読込、LOD読込 |
| Three.js | 3D地形表示 |
| Proj4js | 座標変換 |
| Web Worker | 傾斜量・陰影起伏のバックグラウンド計算 |
| Vite | フロントエンドビルド |
| Nginx | Web配信 |
フォルダー構成
後から水文解析を追加できるよう、表示、DEM読込、解析、3Dを分離しています。
参考のために、ZIPファイルを置いておきますので、ダウンロードして活用してください。
dem_analysis_web_r002/
├── Dockerfile
├── docker-compose.yml
├── package.json
├── index.html
│
├── nginx/
│ └── default.conf
│
├── public/
│ └── data/
│ └── sample_dem.tif
│
├── scripts/
│ ├── convert_to_cog.sh
│ ├── smoke_test.mjs
│ ├── terrain3d_smoke_test.mjs
│ └── large_raster_smoke_test.mjs
│
├── src/
│ ├── main.js
│ ├── config.js
│ │
│ ├── analysis/
│ │ └── terrainMath.js
│ │
│ ├── dem/
│ │ ├── DemReader.js
│ │ └── sampleElevation.js
│ │
│ ├── map/
│ │ ├── baseLayers.js
│ │ ├── createMap.js
│ │ └── demLayer.js
│ │
│ ├── render/
│ │ ├── elevationPreview.js
│ │ └── rasterOverlay.js
│ │
│ ├── ui/
│ │ └── ui.js
│ │
│ ├── utils/
│ │ ├── projection.js
│ │ ├── projectionRegistry.js
│ │ └── rasterLod.js
│ │
│ ├── view3d/
│ │ ├── Terrain3DViewer.js
│ │ └── terrainMesh.js
│ │
│ ├── workers/
│ │ └── terrain.worker.js
│ │
│ └── styles/
│ └── main.css
│
└── docs/
└── qiita_dem_analysis_web_r002.md
DEM解析の処理をanalysis/、大容量DEM用の読込制御をutils/、3D表示をview3d/へ分けています。
r003ではanalysis/hydrology/を追加する予定なので、この段階から地形解析と水文解析を混ぜない構成にしています。
Docker Composeで起動する
起動はDocker Composeから行います。
cd dem_analysis_web_r002
docker compose up --build
ブラウザから次を開きます。
http://localhost:8080/
ソースを変更したあと、古いビルドが残っている場合は一度作り直します。
docker compose down
docker compose build --no-cache
docker compose up
GeoTIFFを読み込む
画面のGeoTIFFを開くからDEMを選択します。
ローカルファイルはBlobとしてGeoTIFF.jsとOpenLayersへ渡します。
const sourceInfo = {
blob: file,
bands: [1],
};
ファイルを読み込むと、
- GeoTIFFメタ情報取得
- CRS判定
- bbox取得
- 表示用LOD作成
- 標高統計計算
- 2D表示
- 3Dメッシュ更新
- DEM範囲へ地図をフィット
という順で処理します。
GeoTIFFの大きさに合わせて表示する
DEMの場所や大きさはファイルごとに異なります。
固定した中心座標やZoomへ移動するのではなく、GeoTIFF自身のbboxを取得し、OpenLayersのView.fit()へ渡します。
fitMapToCurrentDem({ duration = 0 } = {}) {
if (!this.current) return;
const fitExtent = toMapExtent(
this.current.metadata.bbox,
this.current.metadata.projection,
);
this.map.updateSize();
this.map.getView().fit(fitExtent, {
padding: [32, 32, 32, 32],
maxZoom: 19,
duration,
});
}
これで、縦長、横長、狭い範囲、広い範囲のGeoTIFFでも、読み込み後にDEM全体が地図内へ収まります。
画面にはDEM全体を表示ボタンも用意しています。
JGD2011平面直角座標系へ対応する
国土地理院や業務用DEMでは、緯度経度やWeb Mercatorだけではなく、JGD2011の平面直角座標系を使う場合があります。
r002ではProj4jsを使い、
EPSG:6669 ~ EPSG:6687
のJGD2011平面直角座標系I~XIXを登録します。
処理の流れは次のとおりです。
GeoTIFF CRS
↓
EPSGコード判定
↓
Proj4jsへ登録
↓
OpenLayers projectionへ登録
↓
EPSG:3857へ表示
DEMの解析値は元の格子座標を使い、Web地図表示時だけOpenLayersの座標系へ変換します。
標高段彩を表示する
GeoTIFFから取得した標高値を色へ変換し、OpenLayers上へ表示します。
標高
↓
標高範囲を正規化
↓
色へ変換
↓
RGBA
↓
Canvas
↓
PNG Blob
↓
OpenLayers
表示用の標高配列は、元DEMをそのまま使うのではなくLOD化したデータを使います。
これにより、元DEMが大きくても地図操作を軽くできます。
クリック地点の標高を取得する
クリック標高のためにDEM全体をメモリへ展開する必要はありません。
地図上でクリックした位置から、
地図座標
↓
DEM座標系へ変換
↓
row / col
↓
1 x 1 window読込
↓
標高
とします。
GeoTIFF.jsでは次のように1セルだけ読みます。
async readPixel(row, col) {
const raster = await this.image.readRasters({
window: [col, row, col + 1, row + 1],
samples: [0],
interleave: true,
});
return Number(raster[0]);
}
大きなDEMほど、この方式が効いてきます。
3D地形表示
3D表示にはThree.jsを使います。
3D側ではDEM全セルを頂点へ変換しません。
表示専用にさらに間引き、
表示LOD
↓
最大 約180 x 180頂点
↓
Three.js地形メッシュ
とします。
元DEMの解像度と3D表示用メッシュを分離することで、3D描画の軽さと解析データを切り分けています。
3D画面では、
- 回転
- ズーム
- パン
- 標高強調
- ワイヤーフレーム
- 視点リセット
を操作できます。
NoData部分には三角形を作らず、地形に穴が空いた状態として扱います。
傾斜量を計算する
傾斜量は3×3セルを使うHorn法で計算します。
z1 z2 z3
z4 z5 z6
z7 z8 z9
X方向、Y方向の標高勾配を求め、
const dzdx =
((z3 + 2 * z6 + z9) -
(z1 + 2 * z4 + z7)) /
(8 * dx);
const dzdy =
((z7 + 2 * z8 + z9) -
(z1 + 2 * z2 + z3)) /
(8 * dy);
const slopeRad = Math.atan(
Math.sqrt(dzdx * dzdx + dzdy * dzdy),
);
const slopeDeg = slopeRad * 180 / Math.PI;
として傾斜角を求めます。
NoDataを含む3×3セルは解析対象外にします。
陰影起伏を計算する
傾斜と斜面方位からHillshadeを生成します。
DEM
↓
dz/dx・dz/dy
↓
傾斜
↓
斜面方位
↓
太陽方位・高度
↓
Hillshade
r002では傾斜量と陰影起伏を同じWeb Worker内で計算します。
Web Workerへ解析を分離する
傾斜・陰影計算をメインスレッドで実行すると、処理中にOpenLayersやUI操作が重くなります。
そこで、
main.js
↓
terrain.worker.js
↓
GeoTIFF.js
↓
terrainMath.js
と分けます。
メインスレッドは、
地図
UI
3D表示
進捗
を担当します。
Worker側は、
GeoTIFF読込
window分割
傾斜量
陰影起伏
を担当します。
大容量DEMはLODで表示する
大きなDEMで重要なのは、ファイル全体のセル数とブラウザが常時保持するセル数を分けることです。
r002では初期表示用LODに次の上限を設定しています。
preview: {
maxPixels: 1_500_000,
maxAxis: 1600,
resampleMethod: 'nearest',
}
つまり表示時には、
元DEM
↓
最大 約150万セル
↓
OpenLayers / Three.js
とします。
元DEMが3000万セルでも、最初から3000万個のFloat32値をメインスレッドへ展開するわけではありません。
全セル解析はwindow単位で行う
傾斜・陰影を元解像度で計算する場合でも、DEM全体を一括でWorkerへ読み込みません。
GeoTIFF
↓
Window 1
↓
解析
↓
破棄
↓
Window 2
↓
解析
↓
破棄
↓
...
と順番に処理します。
既定では1windowあたり、
最大 約2,000,000セル
を目安にしています。
3×3近傍計算ではwindow境界で計算が切れないように、外周へ1セルのhaloを付けます。
これによりwindow境界でも同じHorn法を使えます。
Strip TIFFとTile / COGを分ける
GeoTIFF内部の格納方法も処理効率に影響します。
r002では、
Strip TIFF
Tile TIFF
COG
を判定し、解析windowの作り方を変えます。
Strip TIFF
Strip TIFFではX方向へ細かく分割しすぎないよう、
と横長windowを使います。
Tile / COG
Tile / COGでは、
のようにタイル構造に合わせて処理します。
高速LOD解析と全セル分割解析
画面では2つの解析方法を選択できます。
解析方式
高速LOD解析
全セル分割解析
用途は次のように分けています。
| 解析方式 | 用途 |
|---|---|
| 高速LOD解析 | 地形の傾向をすぐ確認する |
| 全セル分割解析 | 元解像度の各セルを計算する |
高速LOD解析では縮小DEM上でHorn法を実行します。
全セル分割解析では元DEMをwindowごとに読み込みます。
解析結果も表示用LODへ集約する
元解像度の全セルを計算しても、結果画像まで全セル分をブラウザメモリへ残す必要はありません。
r002では表示用解析結果を、
analysis: {
outputMaxPixels: 3_000_000,
outputMaxAxis: 2048,
}
へ制限しています。
元解像度DEM
↓
全セル計算
↓
傾斜 / 陰影
↓
表示用LODへ集約
↓
OpenLayers
「全セルを計算すること」と「全セル結果を画面用RGBAとして保持すること」を分けています。
読み込み可能なDEMサイズの目安
r002では、GeoTIFFファイルそのものに固定のファイルサイズ上限は設定していません。
GeoTIFF.jsのwindow読込とLODを利用しているため、
ファイル全体
≠
ブラウザが常時保持する配列
となります。
実装上、ブラウザが一度に扱う主な上限は次のとおりです。
| 処理 | 既定上限 |
|---|---|
| 初期表示LOD | 約1,500,000セル |
| LOD最大辺 | 1600セル |
| 全セル解析の1window | 約2,000,000セル |
| 解析結果表示 | 約3,000,000pixel |
| 解析結果最大辺 | 2048pixel |
| 3D表示メッシュ | 最大約180 × 180頂点 |
したがって、元DEMについては数千万セル級でも、LOD表示やwindow分割を使って扱える設計になっています。
Float32 DEMを非圧縮の標高配列として考えると、おおよその容量は次のとおりです。
| セル数 | Float32標高値のみ |
|---|---|
| 1,000万セル | 約38 MiB |
| 2,000万セル | 約76 MiB |
| 3,000万セル | 約114 MiB |
| 5,000万セル | 約191 MiB |
GeoTIFFファイルの実容量は、圧縮方式、NoData領域、Tile / Strip、Overviewの有無によって変わります。
そのため、
最大ファイル容量 = ○○MB
という固定値ではなく、
元DEMは大きくてもよい
↓
表示LODを制限
↓
解析windowを制限
↓
3Dメッシュを制限
という考え方でメモリ使用量を抑えています。
COGへ変換する
通常のGeoTIFFも読み込めますが、Web表示ではCOGにしておく方が扱いやすくなります。
r002には変換用スクリプトを用意しています。
scripts/convert_to_cog.sh
例です。
./scripts/convert_to_cog.sh \
input_dem.tif \
output_dem_cog.tif
中ではGDALを使います。
gdal_translate \
"$input" \
"$output" \
-of COG \
-co COMPRESS=DEFLATE \
-co BLOCKSIZE=512 \
-co OVERVIEWS=AUTO \
-co OVERVIEW_RESAMPLING=AVERAGE
COGにしておくと、
内部Tile
Overview
Range Request
を利用しやすくなります。
NginxからGeoTIFF / COGを配信する
サーバー上へDEMを置く場合に備え、NginxではRange Requestを使える設定にします。
location ~* \.(tif|tiff)$ {
add_header Accept-Ranges bytes always;
add_header Cache-Control "public, max-age=86400" always;
try_files $uri =404;
}
静的ファイル配信側は、
sendfile on;
tcp_nopush on;
etag on;
としています。
ローカルファイル選択だけで使う場合には必須ではありませんが、COGをWebサーバーへ配置するときにそのまま使えます。
画面構成
r002では左側の操作パネルと、2D / 3D表示領域を分けています。
表示モードは、
2D
3D
2D + 3D
を切り替えられます。
r002でできること
ここまでで、ブラウザ上から次の流れを実行できます。
GeoTIFFを開く
↓
DEM全体へ自動ズーム
↓
標高段彩表示
↓
クリック標高確認
↓
2D / 3D切替
↓
傾斜量計算
↓
陰影起伏計算
↓
解析結果表示
大きなDEMでは、
表示 → LOD
3D → メッシュ間引き
標高取得 → 1セルwindow
解析 → window分割
結果表示 → 出力LOD
としているのがr002のポイントです。
r003へ進める内容
r002では、GeoTIFFの読込、2D・3D表示、標高確認、傾斜量、陰影起伏までを一通り実装しました。
この構成を考えていた段階から、次に追加したかったのがDEMを使った水文地形解析です。
傾斜量や陰影起伏だけでも地形の特徴は確認できますが、河川や流域を扱う用途では、
この場所の水はどちらへ流れるか
どこに水が集まりやすいか
どこから河道として扱うか
河道網をどのように整理するか
河道ごとの集水範囲をどう分けるか
までDEMから求められると、利用できる範囲がかなり広がります。
そこでr003では、r002で作ったGeoTIFF読込、LOD、Web Worker、座標変換、OpenLayers、Three.jsの仕組みをそのまま使いながら、水文地形解析を追加していきます。
基本となる流れは次のように考えています。
最初にDEMの局所的な窪地を補正し、その地形から流向を求めます。
DEM
↓
Depression Fill
↓
Flow Direction
次に、各セルへ流れ込む上流セルを積算してFlow Accumulationを求めます。
Flow Direction
↓
Flow Accumulation
↓
集水面積
一定以上の集水面積を持つセルを河道として抽出します。
Flow Accumulation
↓
河道抽出閾値
↓
河道網
ここまではラスター処理ですが、その後は河道網として扱いやすくするため、Strahler谷次数を付与します。
1次谷 + 1次谷
↓
2次谷
2次谷 + 2次谷
↓
3次谷
谷次数が付けば、河道を単純なセル列として見るだけではなく、
小さな支川
↓
支川
↓
主要河道
という河道ネットワークの階層として整理できます。
さらに、この河道網を基準に複数のサブ流域へ分割します。
河道網
↓
Strahler谷次数
↓
合流点・河道出口
↓
谷次数別サブ流域
r002で作った地形表示と組み合わせると、
DEM
├─ 標高
├─ 傾斜
├─ 陰影起伏
├─ Flow Direction
├─ Flow Accumulation
├─ 河道
└─ サブ流域
を同じWeb画面上で切り替えて確認できるようになります。
ラスターだけで終わらせない
r003では河道抽出結果をラスター表示するだけで終わらせず、D8の接続関係を使って河道をLineStringへ変換することも考えています。
河道セル
↓
D8 receiver
↓
源頭・合流点・流出点を判定
↓
河道区間をトレース
↓
LineString
河道区間には、
segment_id
Strahler谷次数
河道長
上流集水面積
始点種別
終点種別
などを属性として持たせます。
谷次数別サブ流域についてもラスター表示だけではなく、ポリゴンへ変換します。
サブ流域IDラスター
↓
境界追跡
↓
Polygon / MultiPolygon
これにより、解析結果をWeb画面で見るだけではなく、GISデータとして次の処理へ渡せます。
GeoPackageへまとめる
河道とサブ流域は、最終的に1つのGeoPackageへまとめる構成にします。
hydrology.gpkg
├─ streams
│ └─ LINESTRING
│
└─ subbasins
└─ MULTIPOLYGON
streamsには河道ネットワークを、subbasinsには谷次数別のサブ流域を保存します。
こうしておけば、Web上で解析した結果をQGISなどへ持っていき、
河道確認
流域面積集計
土地利用との重ね合わせ
降雨データとの結合
流出解析用データ作成
にも利用できます。
r002の仕組みをそのまま使う
r003は別のシステムを作り直すのではなく、r002をそのまま拡張します。
r002
├─ GeoTIFF.js
├─ LOD
├─ window読込
├─ Web Worker
├─ OpenLayers
├─ Three.js
└─ 座標変換
↓ そのまま利用
r003
├─ Depression Fill
├─ Flow Direction
├─ Flow Accumulation
├─ 河道抽出
├─ Strahler谷次数
├─ 谷次数別サブ流域
├─ 河道LineString
└─ GeoPackage
ただし、傾斜量・陰影起伏とFlow Accumulationでは計算の性質が異なります。
傾斜量や陰影起伏は3×3近傍があれば、haloを付けてwindow単位に計算できます。
一方、Flow Accumulationは上流側の寄与を下流へ引き継ぐ必要があります。
上流
↓
↓
↓
下流
windowを完全に独立させると、境界を越えたところで上流側の集水量が失われます。
そのためr003では、大きなDEMの場合、
元DEM
↓
水文解析用LOD
↓
一枚の連続した解析格子
↓
Depression Fill
↓
D8
↓
Flow Accumulation
という構成を使います。
r002で表示用LODを用意したのと同じ考え方で、水文解析には水文解析専用のLODを用意する形です。
このように、r002で作ったDEM表示・地形解析の基盤を残したまま、水の流れ、河道網、サブ流域へ段階的に広げていく予定です。
おわりに
r002では、GeoTIFFをWebブラウザへ読み込んで表示するところから一歩進めて、
2D
3D
標高
傾斜
陰影
大容量DEM
を一つのWebアプリとして扱えるようにしました。
大きなDEMを扱う場合も、すべてを一度にメモリへ載せるのではなく、
表示LOD
1セル読込
解析window
3D表示メッシュ
を用途ごとに分けています。
次のr003では、このDEM基盤を使って河道と流域を作るところまで進めます。
参考リンク
-
GeoTIFF.js
https://geotiffjs.github.io/geotiff.js/ -
OpenLayers
https://openlayers.org/ -
OpenLayers GeoTIFF Source
https://openlayers.org/en/latest/apidoc/module-ol_source_GeoTIFF.html -
OpenLayers COG example
https://openlayers.org/en/latest/examples/cog.html -
Three.js
https://threejs.org/ -
GDAL COG driver
https://gdal.org/en/stable/drivers/raster/cog.html -
Docker Compose
https://docs.docker.com/compose/









