JavaScriptでDEMの水文地形解析をWeb化する：河道次数・流域・GeoPackage出力まで

はじめに

前回記事のr002では、GeoTIFF形式のDEMをブラウザで読み込み、OpenLayersによる2D表示、Three.jsによる3D表示、標高確認、傾斜量、陰影起伏までをJavaScriptで処理できるようにしました。

r002を作っている段階から、次はこのDEMデータ処理基盤を使って、水の流れる方向や集まり方を求め、河道と流域まで扱えるようにしたいと考えていました。

r003では、その続きを実装します。

まず、

DEM
 ↓
Depression Fill
 ↓
D8 Flow Direction
 ↓
Flow Accumulation
 ↓
河道抽出
 ↓
谷次数
 ↓
流域界作成

までをDEMから計算します。

その後、抽出した河道を単なるラスター表示で終わらせず、河道網へStrahler谷次数を付け、谷次数を基準に複数の流域へ分割します。

さらに、解析結果をGISでそのまま利用できるように、

河道：LineString
谷次数別流域：MultiPolygon

へ変換し、1つのGeoPackageへ保存できるようにします。

今回の処理全体は次の流れです。

バックエンドに Python や GIS サーバーは置いていません。DEM の読込、解析、ベクトル化、GeoPackage の生成までブラウザ側の JavaScript で処理します。

今回作るもの

r002の構成を引き継ぐプロジェクトとして、dem_analysis_web_r003 を作成します。
参考のため、ZIPファイルを置いておきますので、ダウンロードして活用してください。

画面では、これまでの DEM 表示、傾斜量、陰影起伏、3D 表示に加えて、次を操作できます。

水文地形解析
├─ 解析解像度
├─ 河道抽出閾値 [km²]
├─ 流域 最小Strahler谷次数
├─ 水文解析を実行
└─ 河道・流域をGeoPackage保存

2D表示では、

Depression Fill
Flow Direction
Flow Accumulation
河道抽出（ラスター）
Strahler谷次数
河道ベクトル
谷次数別流域

を切り替えられます。

1. 河道を抽出

DEM から Flow Accumulation を計算して河道を抽出すると、通常は複数の支川が現れます。

    ┌─ 支川A
────┤
    └─ 支川B ───── 本川
           ┌─ 支川C
───────────┤
           └─ 支川D

河道網を降雨流出解析や流域特性量の計算に使うなら、河道の合流関係に応じて流域も複数に分けた方が扱いやすくなります。

そこで今回は河道へ Strahler 谷次数を付け、その階層を使って流域を作ります。

2. Strahler谷次数

Strahler 方式では、源頭から始まる河道を1次とします。

同じ次数の河道が合流すると次数が1つ上がります。

1次 + 1次 → 2次
2次 + 2次 → 3次
3次 + 3次 → 4次

異なる次数が合流した場合は、大きい方の次数を引き継ぎます。

1次 + 2次 → 2次
2次 + 3次 → 3次

今回の実装では、D8 の receiver 配列から河道セルの上流・下流接続を作り、そのトポロジー順に谷次数を計算します。

概念的には次のような処理です。

for (const streamCell of topologicalOrder) {
  const incoming = upstreamStreamOrders(streamCell);

  if (incoming.length === 0) {
    order[streamCell] = 1;
    continue;
  }

  const maxOrder = Math.max(...incoming);
  const sameMaxCount = incoming.filter((value) => value === maxOrder).length;

  order[streamCell] = sameMaxCount >= 2
    ? maxOrder + 1
    : maxOrder;
}

実際のコードでは各セルで配列を作りすぎないよう、TypedArray と上流カウントを使っています。

3. 谷次数別流域を作る

画面で最小谷次数を選べるようにしました。

1次以上
2次以上
3次以上
4次以上
5次以上
6次以上

考え方は次のとおりです。

1次以上
  → 小さな支川まで含めて細かく流域分割

2次以上
  → ある程度まとまった支川単位

3次以上
  → より大きな流域単位

指定した次数以上の河道について、

より高い次数へ合流する直前
解析範囲の流出端

を流域の出口候補にします。

その後、DEM の各有効セルから D8 receiver を下流へたどります。

DEM cell
  ↓
receiver
  ↓
receiver
  ↓
流域出口

最初に到達した出口の ID をそのセルの subbasinId とします。

これにより、各セルは原則として1つの流域だけへ所属します。

4. 大きなDEMでは水文解析LODを使う

傾斜量や陰影起伏は周囲の数セルがあれば計算できます。そのため window に分割して処理しやすいです。

Flow Accumulation はそうはいきません。

上流 10
  ↓
    35
  ↓
   120
  ↓
下流

上流側の寄与量が下流まで連続して伝わる必要があります。

そこで、大きな GeoTIFF は水文解析だけ専用 LOD を作ります。

設定値は src/config.js にまとめています。

hydrology: {
  fullMaxCells: 1_500_000,
  maxPixels: 1_200_000,
  maxAxis: 1600,
  resampleMethod: 'nearest',
  defaultStreamThresholdKm2: 0.5,
}

大事なのは、LOD に落としたあとも水文解析格子は一枚につながった状態で計算することです。

5. Depression Fill (窪地処理)

DEM に局所的な窪地があると、D8 の流れがそこで止まることがあります。

r003 では Priority-Flood を使って窪地を処理します。

元DEM
 ↓
外周セルをpriority queueへ登録
 ↓
低標高側から展開
 ↓
周辺セルが現在水位より低い場合は持ち上げる
 ↓
Filled DEM

Fill 深も Float32Array で保持し、2D/3D 表示に使えます。

6. D8 Flow Direction (流向)

各セルから8近傍を調べ、最も急な下り方向を receiver とします。

NW   N   NE
 W   C    E
SW   S   SE

Fill 後に同標高の平坦面が生じた場合は、Priority-Flood で記録した排水親セルをフォールバックとして使います。

そのため、単純に「標高差0なら流れなし」とするより平坦部で河道が途切れにくくしています。

7. Flow Accumulation (累積流量)

receiver から各セルの indegree を求め、上流から下流へトポロジカル順に累積します。

源頭
 ↓ 1
 ↓ 3
 ↓ 8
 ↓ 27
出口

セル数だけでなく、解析格子のセル面積を掛けて km² に換算します。

河道抽出閾値を 0.5 km² にした場合は、

upstream_area >= 0.5 km²

のセルを河道とします。

8. 河道ラスターをLineStringへ変換する

河道のベクトル化では、ラスターの外周をトレースして線にする方法は使っていません。

水文解析ですでに D8 receiver があるため、その接続関係をそのまま利用します。

区間は、

source
confluence
outlet

で分割します。

出力する属性は次のとおりです。

segment_id
strahler
length_m
cell_count
upstream_area_km2
start_type
end_type

length_m は座標値の差をそのまま使わず、水文解析時に求めたメートル単位のセル寸法から D8 経路長を積算します。

9. 流域ラスターをMultiPolygonへ変換する

流域は subbasinId という整数ラスタで得られます。

これを GIS で利用できるようにポリゴン化します。

処理は src/analysis/hydrology/subbasinVector.js へ分離しました。

subbasinId raster
 ↓
隣接セルとIDが異なる辺を抽出
 ↓
境界edgeをリングへ接続
 ↓
外周 / 穴を分類
 ↓
MultiPolygon

なぜPolygonではなくMultiPolygonか

通常は1つの流域が1つの Polygon になります。

ただし、ラスタでは、

1 0 1
0 0 0

のように同じ ID が斜めに接する場合があります。

また NoData 領域が内部にあると穴を持つ場合もあります。

そこで出力 Geometry は常に MultiPolygon とし、

複数成分
内部リング

を失わないようにしています。

属性は次のとおりです。

subbasin_id
stream_order
cell_count
area_km2
upstream_area_km2
outlet_x
outlet_y

GeoPackageへ保存するときは出口座標も EPSG:4326 へ変換し、outlet_lon / outlet_lat として保存します。

10. OpenLayersでベクトル表示する

河道と流域はラスター画像ではなく VectorLayer で表示します。

河道GeoJSON
 ↓
OpenLayers GeoJSON format
 ↓
VectorSource
 ↓
VectorLayer

流域も同様です。

MultiPolygon GeoJSON
 ↓
VectorSource
 ↓
VectorLayer

流域は subbasin_id ごとに色を変えています。

元 DEM の座標系が JGD2011 平面直角座標系でも、OpenLayers の読込時に EPSG:3857 へ変換して地図へ重ねます。

11. GeoPackageをブラウザで生成する

GeoPackage 出力には @ngageoint/geopackage を使います。

r003 では1ファイルに2つの Feature Table を作ります。

<DEM名>_hydrology_strahler.gpkg
│
├─ streams
│   geometry = LINESTRING
│
└─ subbasins
    geometry = MULTIPOLYGON

GeoPackage 生成処理は src/export/geopackageExport.js にまとめました。

概略は次のようになります。

const geoPackage = await GeoPackageAPI.create();

createFeatureTable({
  geoPackage,
  tableName: 'streams',
  geometryType: GeometryType.LINESTRING,
  geometryTypeName: 'LINESTRING',
  boundingBox,
  properties: streamColumns,
});

createFeatureTable({
  geoPackage,
  tableName: 'subbasins',
  geometryType: GeometryType.MULTIPOLYGON,
  geometryTypeName: 'MULTIPOLYGON',
  boundingBox,
  properties: subbasinColumns,
});

await geoPackage.addGeoJSONFeaturesToGeoPackage(
  streamFeatures,
  'streams',
  false,
  500,
);

await geoPackage.addGeoJSONFeaturesToGeoPackage(
  subbasinFeatures,
  'subbasins',
  false,
  500,
);

const bytes = await geoPackage.export();

GeoPackage 内の2レイヤは EPSG:4326 としています。

元 DEM が EPSG:6677 などでも、保存前に OpenLayers で EPSG:4326 へ変換します。

12. sql.js WASMをDocker buildへ組み込む

GeoPackage JS はブラウザで SQLite を扱うために sql.js の WebAssembly を使います。

prebuild で WASM を公開ディレクトリへコピーするようにしました。

package.json は次のようにしています。

{
  "scripts": {
    "prebuild": "node scripts/copy_geopackage_wasm.mjs",
    "build": "vite build"
  }
}

コピー先は、

public/vendor/geopackage/sql-wasm.wasm

です。

アプリ側では、

setSqljsWasmLocateFile(
  (file) => `/vendor/geopackage/${file}`
);

として参照します。

これなら Docker の runtime 側は従来どおり Nginx だけで構成できます。

13. GeoPackageの属性

streams

カラム	型	内容
segment_id	INTEGER	河道区間ID
strahler	INTEGER	谷次数
length_m	REAL	区間長[m]
cell_count	INTEGER	セル数
upstream_area_km2	REAL	末端での上流面積[km²]
start_type	TEXT	source / confluence
end_type	TEXT	confluence / outlet

subbasins

カラム	型	内容
subbasin_id	INTEGER	流域ID
stream_order	INTEGER	対応谷次数
cell_count	INTEGER	セル数
area_km2	REAL	面積[km²]
upstream_area_km2	REAL	出口上流面積[km²]
outlet_lon	REAL	出口経度
outlet_lat	REAL	出口緯度

14. 画面操作

基本操作は次の順です。

1. GeoTIFFを開く
2. 河道抽出閾値を入力
3. 最小Strahler谷次数を選択
4. 水文解析を実行
5. Strahler谷次数を確認
6. 河道ベクトルを確認
7. 谷次数別流域を確認
8. 河道・流域をGeoPackage保存

GeoPackage 出力ボタンは、水文解析が完了して河道 LineString と流域 Polygon の両方が生成された場合だけ有効になります。

15. フォルダー構成

dem_analysis_web_r003/
├── Dockerfile
├── docker-compose.yml
├── package.json
├── index.html
├── README.md
├── VALIDATION.md
│
├── nginx/
│   └── default.conf
│
├── public/
│   └── data/
│       └── sample_dem.tif
│
├── scripts/
│   ├── convert_to_cog.sh
│   ├── copy_geopackage_wasm.mjs
│   ├── smoke_test.mjs
│   ├── terrain3d_smoke_test.mjs
│   ├── large_raster_smoke_test.mjs
│   ├── hydrology_smoke_test.mjs
│   ├── stream_network_smoke_test.mjs
│   ├── subbasin_polygon_smoke_test.mjs
│   ├── geopackage_export_smoke_test.mjs
│   └── ui_wiring_test.mjs
│
├── src/
│   ├── main.js
│   ├── config.js
│   │
│   ├── analysis/
│   │   ├── terrainMath.js
│   │   └── hydrology/
│   │       ├── hydrologyMath.js
│   │       ├── hydrologyRender.js
│   │       ├── streamNetwork.js
│   │       └── subbasinVector.js
│   │
│   ├── dem/
│   │   ├── DemReader.js
│   │   └── sampleElevation.js
│   │
│   ├── export/
│   │   └── geopackageExport.js
│   │
│   ├── map/
│   │   ├── baseLayers.js
│   │   ├── createMap.js
│   │   ├── demLayer.js
│   │   ├── streamVectorLayer.js
│   │   └── subbasinVectorLayer.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
│   │   └── hydrology.worker.js
│   │
│   └── styles/
│       └── main.css
│
└── docs/
    └── qiita_dem_analysis_web_r003.md

解析ロジック、表示、GIS出力を分離しているので、後から処理を差し替えやすくしています。

16. 起動

Docker Compose で起動します。

cd dem_analysis_web_r003
docker compose up --build

ブラウザで開きます。

http://localhost:8080/

以前のビルドが残っている場合は、

docker compose down
docker compose build --no-cache
docker compose up

として、ブラウザ側も Ctrl + F5 で再読み込みします。

17. テスト

npm test

今回追加した主なテストは2つです。

流域ポリゴン化

subbasin polygon smoke test: OK
MultiPolygon / diagonal components / holes: OK

通常のポリゴンだけでなく、

斜め接触
複数成分
穴

を確認しています。

GeoPackage出力

geopackage export smoke test: OK
streams=LINESTRING / subbasins=MULTIPOLYGON: OK

Feature Table の Geometry Type、出力ボタン、ビルド時 WASM 配置などをチェックします。

また、HTML と JavaScript の DOM ID が食い違うと実行時エラーになりやすいため、ui_wiring_test.mjs で ID の欠落と重複も検査します。

18. 今後やりたいこと

GeoPackage まで出せるようになると、DEM 解析結果を QGIS 等へ渡しやすくなります。

次に追加したいのは、流域ごとの地形特性量です。

流域
├─ 面積
├─ 平均標高
├─ 最大標高
├─ 平均傾斜
├─ 主河道長
├─ 河道密度
├─ 流域形状係数
└─ 出口標高

ここまで揃えば、降雨流出モデルへ渡す前処理データとしてかなり使いやすくなります。

また、地図上で任意の河道点をクリックし、その点を流域出口として上流域を再抽出する機能も追加したいところです。

参考