はじめに

GrassGIS,QGISを使って国土地理院の基盤地図情報の数値標高モデルから集水域や流路を作成する方法を紹介します。いくつか同様のwebページはあるのですが、かゆい所に手が届かない感じがあるので、書いてみました。

地理院DEMの変換

　MIERUNE Inc.が公開しているプラグイン「QuickDEM4JP」で、国土地理院の標高データ（XML）をラスタ（tif）に変換します。この際、CRSは地理座標系のWGS84を指定します。

　DEMが複数図枠の場合は、それぞれで変換して、「結合（godal_merge）」で結合します。
　その後、「再投影（warp）」でリサンプリング法とnodata値、変換先のCRS単位での解像度を指定して投影座標に変換します。プロファイルに高圧縮を指定するとラスタのファイルサイズを小さくできます。
　最初に地理座標で変換するのは、元の標高データが地理座標なので、個別に投影変換すると複数図枠の境界が合わなくなること、セルサイズを切りのいい5mにするためです。

nodata値の穴埋め（Fill NoData）

　5mDEMは水域がnodataになっており、水文解析をするためにはgdalのFill NoDataを使って河川を埋めます。そのままだと海まで埋め立てられてしまいますので、埋める範囲を指定するラスタレイヤを作成して、処理範囲を指定します。
入力データ　国土地理院5mDEM

処理範囲のポリゴンの作成

「区分表で再分類」を使ってデータ範囲が1のラスタを作成します。

　「ラスタをベクタ化（polygonize）」でポリゴンに変換します。

　「孔を削除」でポリゴンの中の孔を埋めます。

　これで、海を除いた範囲のラスタ（処理範囲）のポリゴンができました。

nodata値の穴埋め

　Fill NoDataを使って、NoDataを埋めます。内挿値を検索する距離（ピクセル単位）は、埋めたい NoDataの幅/セルサイズ/2 より大きい数字を指定します。/2は両側から埋めるので、半分で良いからです。

　海まで埋まってしまいますので、「マスクレイヤで切り抜く」を使って、処理範囲ポリゴンで切り抜き海をNoDataに戻します。その際、「マスクレイヤの領域に切り抜き範囲を一致させる」のチェックを外します。これは、チェックが入っていると、マスクレイヤが領域になり、ラスタの原点がずれるためです。

　これで、陸地の水域に標高値が入った標高ラスタができました。

凹地の平滑化（r.fill.dir）

　DEMに凹地（周りの土地より低い場所）があると、水の行き先がないので流路が途切れてしまいますので、凹地を埋める（標高値を高くする）処理が必要とされています。grassのr.watershedでは、この処理が組み込まれているので不要ですが、項目だけ挙げます。

流向ラスタの作成

　水は高い方から低い方に流れるのでDEMから各セルの流向を求めることができます。また、各セルに単位雨量が降るものとして、流向に基づいて下流側に累積させることにより、累積流量を求めることができます。
grassではDEMから流向、累積流量ラスタを作成するプロセッシングが２つあります。

r.terraflow
r.watershed
　r.watershedのドキュメントによると、r.terraflowは、ArcGIS（のFlow Direction）と同様、つまり、DEMの傾斜方向で流向を決定しています。一方、r.watershedはDEM データエラーの影響を最小限に抑えるように設計された A^T least-cost search algorithmを使用しているとのことです。grassでは後の処理でr.watershedの出力形式を使用するので、基本的にはこちらを使えばいいと思いますが、r.terraflowの方が処理内容は理解しやすいです。

単一流向 (SFD、D8) と複数流向 (MFD)

　累積流量の計算において、SFDとMFDがあり、MFDは、隣接セルへの傾斜を比例配分の重み付け係数として使用して、標高の低いすべての隣接セルに分配されるのに対し、SFDは１方向のみに分配されます。
　r.terraflow、r.watershedともMFDがデフォルトですが、流向は1方向に決定され、流域は流向から作成するため、流域作成との整合性を考えるとSFDを使用するのが良いと思います。

流向ラスタなどの作成（r.watershed）

　入力は、標高ラスタ、外部流域の最小サイズ、単一フローを使う　を指定します。外部流域の最小サイズは、河川が発生する地点（源頭）の集水面積（セル単位）で、小さな値を指定するほど、河川が上流に延びます。
　出力は、排水先セル数（＝集水面積、累積流量）、排水方向（＝流向）、河川セグメントを指定します。

流域の作成

　河川上の任意地点（例えば流量観測所）の流域（集水域）は、その地点から流向ラスタを上流にたどることにより求められます。
流域を作成するr.water.outletは任意地点の座標と流向ラスタを入力とし、流域ラスタを出力します。さらに、r.to.vectを使ってベクタ化することにより流域ベクタ（ポリゴン）が得られます。

流域ラスタの作成（r.water.outlet）

　排水方向ラスターと出口の座標を指定します。出口の座標は河川セグメントラスタを参照して河川のセルの座標を指定します。

流域ベクタ（ポリゴン）の作成（r.to.vect）

　入力ラスタ（流域ラスタ）、地物のタイプ（area,デフォルト）を指定します。面の角を滑らかにするをチェックするとセルの角がとれます。

流路ベクタ（ライン）の作成

水が地形によって集まって流れるところが流路で、つまり、累積流量が任意の閾値以上のセルをつないだものです。
　grassgisでは、r.stream.extractを使用して流路ベクタ（ライン）が得られますが、DEMから独自に計算するので、r.watershd等の流向と整合しないのが難点です。
　また、河川セグメントラスタからr.to.vectでライン化すると流路が並走するところで線が団子状になってしまいます。細線化（r.thin）する方法もありますが、流向・累積流量と関係なく流路を変更することになります。
　流路のベクタ化の考え方は流向に基づいてラインを引くというシンプルなものなので、QGISの「式によるジオメトリ」で作成する方法を考えました。

河川のポイント化、ラスタ値付加

河川セグメントラスタを「ラスタをベクタ化（pixels to points）」でポイントに変換します。河川（累積流量が閾値以上）のセルがポイントになり、閾値以下のセルはnodataなので出力されません。
次に、ポイントにラスタ値付加をします。「ベクタレイヤにラスタ値付加」を使って、流向ラスタ(fd1)、累積流量ラスタ（fa1)を付加します。

式によるジオメトリで流路ラインを作成

　流向ラスタの方向は、右上を1として反時計回りで定義されています。
https://grass.osgeo.org/grass-stable/manuals/r.watershed.html

　流向ラスタ（fd1）を使って、ポイントから1ピクセル分の流路ラインを作成します。

/*
１ピクセル分の流路ラインを作成、セルサイズ5
*/

make_line( @geometry,
case 
when   "fd1" = 1 then   translate( @geometry,5,5)
when   "fd1" = 2 then   translate( @geometry,0,5)
when   "fd1" = 3 then   translate( @geometry,-5,5)
when   "fd1" = 4 then   translate( @geometry,-5,0)
when   "fd1" = 5 then   translate( @geometry,-5,-5)
when   "fd1" = 6 then   translate( @geometry,0,-5)
when   "fd1" = 7 then   translate( @geometry,5,-5)
when   "fd1" = 8 then   translate( @geometry,5,0)
else
 null
end
)

累積流量を流路の太さで表現する

　累積流量（fa1フィールド）を使って、累積流量（集水面積）をシンボルで表現することができます。ここでは、連続値による定義で、累積流量が大きいほど太くしました。fa1はベクタの範囲外からの流入があると値が負になるので、値を絶対値としています。

おわりに

　最初は、「式によるジオメトリで流路ラインを作成」の所だけ書くつもりだったのですが、長々と書いてしまいました。読んでいただいた方ありがとうございます。