ドローンを使うことで地形測量は大幅に省力化されましたが、水面下の地形は屈折の影響で実際より浅くなってしまいます。この影響をQGISでできるだけ簡単に補正する方法をお教えします。
川底まで見えるドローンの写真を用意しよう
ドローンの写真から川底の地形を再現するには、当然川底が見える浅い／きれいな河川の写真でなければなりません。本稿で使用した写真を以下に示します。65枚の上空写真と7点のGCPを用いてMetashapeで作成したオルソ画像とDEMです。
拡大してみると、下の写真のように川底の小石までよく見えていることがわかります。川底の小石などが見えないとステレオ視で距離を計算するための対象点が見つかりませんので、補正する以前に川底の地形がそもそも信用できず、本稿の方法は適用できません。
ドローンの写真を使った河川測量精度をあげるためには、水深やにごりに注意する他に、太陽の反射が少ない曇りの日や太陽高度が低い時間帯を選んだり、波が小さくなる風が弱い日を選んだりした方がいいよ。
水面の屈折の影響
水面の屈折により水底が浅く見える現象は、みなさん小学校で習ったのではないでしょうか。どれくらい浅く見えるかは視点や距離によって変化するので正確な議論は難しいのですが、ここではドローンを用いて実験的に求めた補正係数1.42を採用することにします。つまり、ドローンから見える水深より実際の水深は1.42倍深いということです。
水面の位置を求める
水深を使って補正するのですから、まず水面の位置を求める必要があります。そのため、オルソ画像をよく見て、水際が判別しやすい場所にポイントを打っていきます。勾配が急な場所ではポイントの位置が少しずれるだけで標高が大きく変化するので、できるだけなだらかな地形の水際を選んでポイントを打ちます。下図では新しくwatersurfaceというシェープファイルを作成してポイントを打っています。整数フィールドidと実数フィールドwsl (water surface level; 水面高の意味)をシェープファイル作成時に定義してあるのですが、ポイント作成時には何も入力しなくても構いません。下段が完成した様子です。（本稿ではidは最後まで使用しません。）
水面高をDEMから読み取る
作成したwatersurface.shpのポイント位置の標高をDEMから読み取ります。
下図のように属性テーブルを開いて①フィールド計算機を起動し、水面高を保存するフィールド②wslを選択して③④「式」に以下を入力します。
raster_value( '230117_Metashape_DEM_2db2ba43_cfd0_40cf_aee9_4d5725421fe4',1, $geometry )
といっても手で入力する必要はなく、中央の⑤リストボックスから選択してダブルクリックしていけば良いのです。raster_valueは「ラスタ」カテゴリ、Metashape_DEMは「地図レイヤ」カテゴリ、$geometryは「ジオメトリ」カテゴリの中にあります。第二引数の「1」とその前後のコンマだけキーボードから入力します。「OK」をクリックすると、下図右のように点の位置のDEMの数値が読み取られます。
水面高が同じとみなせる横断線を作成する
次に、ラインシェープファイルws_xsections.shpを作成し、watersurface.shpのポイントに交わり流れに垂直な横断線をマウスで描いていきます。ws_xsections.shpにもwslフィールドを作成し、値はwatersurface.shpから手作業でコピーします。フィールド計算機で自動的にコピーする方法もありますが、上流のwslが下流より低いなどの不都合があればこの段階で修正するために敢えて手作業で値を設定しています。
横断線を元に水面高ラスタを作成する
TIN内挿ツールを使います。入力ベクタは「ws_xsections」、内挿対象の属性を「wsl」とし、ws_sectionsのwslを「ブレークライン（改行）」、内挿方法は「線形」を選びました。領域は「ws_sections」に合わせ、出力レイヤで指定したファイルには水面高ラスタが出力され、TIN内挿出力に指定したファイルにはTINの形状を示すシェープファイルが出力されます。
実行結果は以下のようになります。
DEMラスタを再投影する（オプション）
本稿で作成した地形は流況計算プログラムの地形データとして使用することを想定しています。その目的にはDEMの解像度が高すぎますので、再投影してDEMの解像度を小さくしておきます。解像度を小さくする必要がなければこの節は読み飛ばして結構です。
もとのDEMのセルは0.6mm角ですので、ここでは2cm角に変更してみます。GDALの「再投影（warp）」ツールを使い、下図のようにパラメータを設定します。2cm角の設定は、「変更先CRSの単位での解像度[オプション]」で行っています。変換先CRSはEPSG2445（平面直角座標系）ですので、単位はmです。再投影ツールのパラメータ設定はたくさんあるので、下図では本来１つの縦長ダイアログであるものを２つに分けて表示してあります。
DEMラスタをポイントシェープファイルに変換する
「ラスタをベクタ化（pixcels to points）ツール」を使います。前節の再投影を行った場合は再投影されたDEMを、再投影しなかった場合は元のDEMを「ラスタレイヤ」に設定します。「点ベクタ」には出力するポイントシェープファイル名（ここではMetashape.shp）を設定します。
このツールに限らず、出力ファイル名を指定するときには、右端のボタンをクリックして適切なフォルダを選んでね。ファイル名を直接ボックスに入力してしまうと、どこに保存したかわからなくなったり、書き込み禁止のフォルダに出力してエラーになったりするわよ！
出力されたシェープファイルは以下のようになります。デフォルトではポイントの径が大きく、またポイントのアウトライン（縁取り）の黒が目立ってしまい。全体に真っ黒になってしまうので、右下のレイヤプロパティからVALUE（DEMの値＝地盤高）に基づく色分け表示をするように設定し、またシンボルの定義を変更してアウトラインを透明、ポイント径を0.5mm程度に小さくしてあります。中央下はMetashape.shpの属性テーブル、左下は左隅の拡大図です。ここまで拡大すると、ポイントの背景にオルソ画像が透けて見えます。
ポイントシェープファイルから地形データを作成する
さあ、ここからが本番です。作成目標となる地形データは、X, Y, Z座標のcsvファイルです。
この節の作業では、膨大なポイント数のデータをフィールド計算機を使って書き換えたり削除したりするよ。作業の節目節目で属性テーブルを保存してね！保存しないと、作業の取り消しができるように過去のデータを記憶しようとしてQGISがパンクすることがあるよ。
X, Y座標を作成する
Metashapeレイヤの属性テーブルでフィールド計算機を起動し、Xフィールドを新規作成してX座標（$x）を書き込みます。Y座標も同じように操作します。
水面高を水面高ラスタから読み込む
同じくフィールド計算機を使い、水面高ラスタレイヤ（wsl_rst）から該当する位置の水面高を読み込みます。この操作は本稿の「水面高をDEMから読み込む」と同じですね。
水深を求める
水深は水面高から地盤高を減算すれば求まります。負の値は水面より上に地盤があるということですから、０にします。式は以下のようになります。このif文の文法は、Excelと同じく、第一引数が条件、第二引数が条件が真の時の値、第三引数が条件が偽の時の値です。
if ( "WSL" - "VALUE" >=0, ("WSL" - "VALUE") * 1.42, 0)
下図のように設定することで、計算された水深は小数点付き数値の新フィールド「depth」に書き込まれます。
補正された地盤高を求める
水深が正値であれば、水面高ー水深が河床の地盤高、水深が正値でなければ、元のDEMの値が地盤高になります。
if ("depth" > 0, "WSL" - "depth", "VALUE")
最後に、不要なフィールドは消してしまいます。
以上で、水面の屈折を補正した地形データがMetashapeレイヤの属性テーブルとしてできあがりました。
地形データをCSVファイルとしてエクスポートする
Metashapeレイヤを右クリックし、「新規ファイルに地物を保存」でCSVファイルに出力できます。
できあがった地形データを流況計算ソフトにインポートしてみる
最後に、地盤高を流況計算ソフトiRICに読み込んでみました。
オブジェクトブラウザの「地理情報」ー「Elevation」を右クリックして「インポート」－「点群データ」を選ぶと「インポート設定」ダイアログが開きます。下図左のようにCSVファイルを選ぶと、下図右のように地形データが読み込まれました。
工数が多くて大変でしたね。でも１つづつ順番にやっていけば大丈夫。Have fun!
コメント