DEMから流域界と流路を検出するニーズは多くあると思います。できるだけ簡単に流域界と流路を検出し、利用しやすいようにベクタデータにする方法を教えます。QGIS3.40LTR対応。
QGISでの流域界の検出は、これまでSAGAツールを使用することが多かったけど、QGIS3.30以降SAGAツールがデフォルトでは使えなくなったんだ。ここでは、GRASSツールを使う方法を解説するよ。SAGA GIS自体はQGISと一緒にインストールされているので、SAGA GISで流域界を作成することは今でもできるよ。
QGISはフリーでオープンソースの地理情報システムです。
DEMの窪地を埋める
流域界と流路を検出する原理は、DEMの勾配の方向を求めることで水が流れる方向を決めていくことです。しかし、窪地があるとそこに流れが滞留してしまいます。（現実世界のように水が溜まってあふれ出すような複雑な計算はしていません。）ですので、最初に窪地を埋めてしまう必要があります。
例として使用するDEMは国土地理院からダウンロードした与論島です。座標系は投影座標系のUTMか平面直角座標系にします。地理座標系ではうまくいきません！
窪地を埋めるために使うツールはGRASSのr.fill.dirを使います。下図のように「標高」にDEMレイヤを、「方位角の形式」にgrassを、窪地なしDEMの出力に適切なファイル名を設定して「実行」します。他にもいくつかの出力がありますが、使用しませんのでここでは「一時ファイルに保存」にしてあります。
r.fill.dirを実行すると、ログにERROR6が赤字で表示されることがあるけど、このエラーは無視して良いです。実行結果は以下のようになりました。
このままではオリジナルとの違いが分かりにくいので、ラスタ計算機でオリジナルとの差を作ってみましょう。流域界の作成とは関係ない説明のための作業なので、実務ではスルーしていただいて結構です。
上図右が減算結果で、グレー～白色で表示されている部分がオリジナルのDEMより高くなった場所（＝埋められた窪地）です。
流域界、流路を生成する
GRASSのr.watershedツールを使います。入力する「標高」は、前節で作成した窪地を埋めたDEMです。このツールは１回の実行でいろいろな出力を行いますが、ここでは「排水先セルの数」、「排水方向」、「流域ラベル」、「河川セグメント」をそれぞれaccumulation.tif, flowdir.tif, basin.tif, segment.tifというファイル名に出力させます。
パラメータとして「外部流域(exterior basin)の最小サイズ」を設定する必要があります。サイズの単位はDEMのセル数です。外部流域とは、解析範囲の端（海岸線など）から外に水が流れ出す地域のことで、この値を小さくすると、小さな流域は無視されます。ここでは1000にしてみました。
このパラメータは河川上流の分岐を判定する集水セル数の閾値にも使われていて、大きな値にすると流域が大括りに、小さくすると細かく分割されます。上流部の流域界の分割にも強く影響するので、「外部流域」というネーミングはふさわしくないですね。
ふさわしくないと言えば、出力の「排水先セルの数」はそのセルに水が流れ込むセルの累積数だし、「流域ラベル」は「検出された流域」のことだし、ネーミングが直感的とは言えないね。使って慣れてね！
出力されたラスタを以下に示します。
セル数が1000以上となる流域に分割され、それぞれの流域にはそれぞれ1本の河川セグメントが生成されていることがわかります。
なお、「排水方向(flowdir)」は本解析では使用しませんが、任意の点から上流の流域を検出する際に使用するので出力しました。上図を見ると海岸線のflowdirは負値になっていることがわかります。flowdirの絶対値は流向（北東を1として反時計回りに東が8）を表し、正負で外部に流出しているかどうかを表しています。
流域界・流路をベクタ化する
前節までで流域界・流路は検出できましたが、ラスタのままよりベクタにした方が使いやすいこともあると思いますので、やっておきましょう。
流路のベクタ化
ツールボックスのGRASSにあるr.stream.extractを使います。厳密にいうと、前項で作成した河川セグメントラスタをベクタ化するのではなく、改めて窪地を埋めたDEMから流路ベクタを作成します。
入力は「標高ラスタ」に窪地を埋めたDEM、「累積流量ラスタ」に前項で作成したaccumulationを設定します。「河川の最低累積流量」は前項の「外部流域(exterior basin)の最小サイズ」と同じ1000を設定します。出力できるものはいろいろありますが、ここでは「河川ID のベクタ」だけに出力ファイル名を設定しましょう。詳細パラメータを開き、「v.out.ogr出力タイプ」をlineにすることを忘れないで！
出力ファイルの形式としてはかならず.gpkgを選んでね。シェープファイルも選べるけど、何も出力されないわ。.gpkg形式が嫌なら、あとでシェープファイルにエクスポートすればいいのよ。
流域界のベクタ化
QGISのメニューバーから「ラスタ」ー「変換」ー「ラスタをベクタ化(polygonize)」を使います。入力レイヤにはbasinを、出力ファイルとして「ベクタ化」枠にファイル名を設定します。
このツールではシェープファイルも出力できるけど、流路にあわせて.gpkg形式にしてみたわ。
以上です。意外に簡単でしょう？Have fun!
コメント