土砂崩れや浸食などで流出した土砂量をドローンなどで作成したDEMを用いて求めるようなケースを想定しています。ラスタ計算機, ラスタのポリゴンによる切り抜き、ラスタのサーフェス体積ツールを使います。
QGISはフリーでオープンソースの地理情報システムです。
かわのさん
かわのさん：昨今、災害が多いから、使うことが多そうなテクニックね。
すだくん：災害だけでなく、植物の伸長速度の計測とか、工事の土工量の計算とか、応用できる場面はいろいろあるよ。
DEMの準備
下の２枚のDEMは、河川内にひっかかった人工物が、洪水時に高水敷を壊しながら移動した前後をUAVで撮影したものです。DEMに影響する河岸植生は可能な範囲で刈り取ってあります。対象地点周辺にGCP（標定点）を多数配置して撮影し、できるだけ正しくジオリファレンスされるよう、DEM作成時には注意して作業する必要があります。（これはドローン写真の合成技術の問題で、本稿とは関係ありませんが・・）また、体積を求めるためにはDEMのCRSは投影座標系である必要があります。この例ではJGD2000 ／ 平面直角座標系になっています。
マスクとなるポリゴンを作成
ドローンで作成したDEMは辺縁部に歪みがあることが多いですし、河岸の急勾配部分もどうしても位置精度が悪くなりますので、そうした部分は排除し、注目したい範囲だけを囲んだポリゴンを作成してマスクを作ります。
直立護岸とか崖とか高低差が大きい部分は、ほんの少し水平位置がずれただけで土砂崩れなどが起きていなくても２枚のDEMの標高差が大きく検出されてしまうんだ。だから最初からそんな部分は計算に入れないようにした方が良いし、もっと言うと、GCPもそんなところに置かない方がいいよ。GCP自体のマッチングの誤差が大きいと、データが信用できなくなるからね。
メニューから「レイヤ」ー「レイヤを作成」ー「ベクタレイヤを作成」を選び、プロジェクトと同じCRSのポリゴンを作成します。フィールドは重要ではありませんのでデフォルトのままとします。
デジタイイングツールバーの「編集モード切替」ツール（ツールバー下段の鉛筆アイコン）をクリックして編集モードに入り、「ポリゴン地物を追加」ツールをクリックして護岸部分が含まれないように注意しながら計算したい範囲を囲んでいきます。囲み終え、右ボタンをクリックしたら、最後にフィールド値を聞かれますので、とりあえず１としておきます。（なんでも構いません。）このあと、「レイヤ編集内容の保存」（編集ツールバーのディスクアイコン）を実行しておくことを忘れないでください。保存しなくても続けて作業はできますが、ポリゴンでラスタをクリップする時にエラーが発生して原因がわからず苦労することがあります。
高水敷崩壊後のDEMとも比べてみて、問題ないことを確かめておきます。
DEMの減算
ラスタ計算機を使い、土砂流出前ー土砂流出後を計算します。
下図のような結果が得られます。
シンボロジで単バンド疑似カラーを選び、手動で±0.05, 0.5, 1, 1.5, 2の範囲で色を付けてみました。色も手動で設定しています。緑の部分はほとんど変化していないところ、黄色からオレンジの範囲は0.5から1m浸食を受けた部分です。水色から青は堆積した部分です。
シンボロジでうまく色をつけると変化が読み取りやすくなるので、ぜひ工夫してほしいところね。デフォルトのままでは変化に気が付かないこともよくあるものね。
マスクで減算したDEMを切り抜く
「ラスタ」ー「抽出」ー「マスクレイヤによる切り抜き」で、注目範囲を切り抜きます。下図ではnodata値としてゼロを与えていますが、よく使われる-999999999.0の方が良かったと思います。どちらの値を設定したとしてもnodata値として判別され、次のステップの計算からは除外されますが、ゼロは元のラスタにデータとして含まれている可能性があるので、体積には影響しなくても面積には影響する恐れがあるからです。ちなみに、もしマスクレイヤーが保存操作されていないと、シェープファイルの中身が空っぽということになり、このステップでエラーが発生します。
下図の白黒画像部分が切り抜かれた標高差のラスタです。
土砂量を計算する
プロセシングツールボックスから「ラスタのサーフェス体積」を起動します。このツールは、基準値として設定した値とDEMの値との差をとって積分し、体積を求めます。方法として以下の４つが選べます。
- 「基準値を超える場合だけ」
- 「基準値を下回る場合だけ」
- 「基準値を下回る場合は減算」
- 「基準値を下回る場合は加算」
１.の体積は正の値、２.の体積は負の値になります。３.は、基準値を下回った負の体積を１.の正の体積から減算するので１.の体積より大きくなりそうですが、実際には小さくなります。４.は逆に体積が１.より大きくなります。本稿のように流出前後のDEMを減算したラスタを使い、基準値をゼロとした場合、1.は浸食域の体積、2.は堆積域の体積（負値）、3.は浸食域の体積から堆積域の体積（正値）を減算したもの、4.は浸食域の体積と堆積域の体積（正値）を加算したものになります。どの方法をとるかは目的によって異なりますが、本稿では1.「基準値を超える場合だけ」を選びました。
これはちょっとわかりにくいわね。大事な解析なら、上の記事をよく読んでそれでも自信がないようだったら、簡単なラスタデータを準備して実際にやって動作を確認しておいた方がよいかもね。
なお、出力されるのは体積だけではなく、方法で選んだ条件を満たすピクセル数と面積も出力されます。計算結果はhtml形式もしくは表として出力されますが、どちらも同じ値が記録されていますので、一方だけで構いません。
計算結果は下図のようになりました。切り抜いたラスタも切り抜き前と同じ 単バンド疑似カラーにしてあります。左のダイアログは計測ツールで計算対象の長径、短径を概則したもの、右側のhtmlおよびテーブルは計算結果です。先述のようにhtmlとテーブルは同じ値になっています。長径約8.5m×短径約2.5mの矩形内の2/3（≒14m2）が注目範囲で、さらに水色～青色の範囲は方法1.ですから堆積域で計算から除外されます。浸食深の平均は色分け図の目算で0.3m前後と思われますので、領域面積10.6m2, 体積2.76m3は妥当な結果だと思われます。
どんな計算でも結果の妥当性をこのような概算で確認しておくことは大切だね。「ラスタのサーフェス体積」のように動作の説明がわかりにくい機能を使うならなおさらだよ。
ちょっと機能がわかりにくいところもあったけど、落ち着いて手順を踏めば大丈夫でしょう？ Have fun!
コメント