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/