チュートリアル: 井戸データとDEMから地下水位マップを作成
| Site: | OpenCourseWare for GIS |
| Course: | 水文地質学のためのGISトレーニング |
| Book: | チュートリアル: 井戸データとDEMから地下水位マップを作成 |
| Printed by: | Guest user |
| Date: | Thursday, 25 April 2024, 2:04 PM |
Description
1. 概要
このチュートリアルを終えると、以下のことができるようになります。
- 井戸（Borehole/Wells）のデータセットをSDIからQGISに読み込む
- ポリゴン内の井戸の密度を計算する
- SRTM 1-Arc Second DEMのダウンロード
- DEMをクリップして再投影する
- DEMのスタイル設定
- 井戸の位置の標高をサンプリングし、標高属性と比較する。
- 属性テーブルで補正・計算を行う
- 井戸の地下水位をラスタに補間する。
2. GeoNode SDIから井戸のデータセットを読み込む
まずはGeoNodeのSDIからQGISに井戸のデータセットを読み込むところから始めます。
このチュートリアルではGeoNodeのSDIであるOrange-Senqu River Basin GIS Serverから取得するBorehole Database in the Stampriet Transboundary Aquifer（Stampriet越境帯水層の井戸データベース）を使用します。
1. Borehole Database in the Stampriet Transboundary Aquiferのメタデータをチェックします。 詳細情報タブを見て、レイヤーの属性を確認します。
2. QGISで新規プロジェクトを立ち上げます。
3. ツールバーのデータソースマネージャーを開く ボタン 
をクリックします。
4. GeoNode タブを選択します。
5. 新規 ボタンをクリックし新しいサービスの接続情報を作成します。
6. 新しいGeoNode接続を作成ダイアログで名前にORASECOMと入力し、http://gis.orasecom.orgをURLにコピーします。
7. 接続テストをクリックします。
テストが成功すると、このポップアップが表示されます。
接続に失敗する場合は、インターネット接続環境やURLを確認してください。
8. OKをクリックしてポップアップを閉じます。
9. ダイアログのOKをクリックして閉じると、新しい接続が追加されます。
10. 接続ボタンをクリックします。
11. ポップアップでエラーが表示されますが、無視してください。エラーメッセージを削除するには、OKをクリックします。
すると、GeoNode上のレイヤーが一覧表示されます。
Stampriet帯水層のデータを使用します。
12. フィルタにstamprietと入力します。
13. stasboreholes WFS レイヤをクリックします (上の図で示されています).
WMSデータのような画像だけではなく、データ自体にアクセスしたいので、WFS Webサービスを使用します。
14. 追加をクリックすると、地図キャンバスにレイヤーが追加されます。
15. ダイアログを閉じます。
インターネットからの読み込みには時間がかかります。
読み込みが完了するまでお待ちください。すると、このように表示されます。
次のセクションでは、このWFSレイヤーをGeoPackageに保存し、ローカルで使用できるようにします。
3. WFSデータをGeoPackageにエクスポート
井戸のデータセットは、まだWFSフォーマットで、緯度・経度の座標を持つ地理座標系（GCS）になっています。標高を補間するためには、このファイルをローカルに保存し、UTM Zone 34S/WGS-84に再投影する必要があります。
1. レイヤパネルでgeonode:stasboreholesレイヤを右クリックします。
2. エクスポート | 地物の保存...を選択します。
3. ベクタレイヤを名前をつけて保存...ダイアログで形式にGeoPackageを選択します。
ボタンをクリックしハードディスクでGeoPackageを保存したいフォルダ(e.g. Z:\Stampriet)を参照します。Stampriet_Data.gpkgというファイル名にします。レイヤ名はBoreholesと入力します。
4. 
ボタンをクリックし、出力するレイヤの投影法を変更します。
5. 座標参照系の選択ダイアログでフィルタに32734（EPSGコード）と入力し、WGS 84 / UTM Zone 34S投影法を検索します。投影法の名前をクリックしOKを押下しダイアログのウィンドウに戻ります。
6. OKをクリックすると、エクスポートが実行されます。
7. After the Stampriet_Data Boreholes レイヤが地図キャンバスに追加されたら、geonode:stasboreholesレイヤを右クリックして、レイヤの削除...を選択し、確認画面でOKを押してレイヤを削除します。
これで、Boreholesレイヤーのローカルコピーができました。
8. プロジェクトの OTF投影を変更します。レイヤパネルで Stampriet_Data Boreholes を右クリックし、レイヤのCRS｜レイヤのCRSをプロジェクトのCRSに設定を選択します。
9. プロジェクトを保存します。
次のセクションでは、特定のエリアの井戸を選択します。
4. 特定のエリアにある井戸を選択
このチュートリアルでは、Stampriet周辺の井戸にのみ関心があります。
想定する境界線のポリゴンは、チュートリアルのデータに含まれています。
1. Stampriet_boundary シェープファイルをプロジェクトに追加します。
2.レイヤパネルでStampriet_boundaryを右クリックしレイヤの領域にズームを選択します。
3. 境界線を赤い枠線でスタイリングします。
ここで、境界内にある井戸を切り抜いて、GeoPackageに保存します。
4. メインメニューから ベクタ | 空間演算ツール | 切り抜く（clip）...を選択します。
5. 切り抜く（clip）ダイアログで入力レイヤとしてStampriet_Data Boreholesを、オーバーレイレイヤとしてStampriet_boundaryを選択します。切り抜いた結果をStampriet_Data GeoPackageにというboreholes_clippedレイヤ名で保存します。
6. 実行をクリックします。処理後にダイアログを閉じます。
7. レイヤパネルからStampriet_Data Boreholesレイヤを削除します。
次のセクションでは、このエリアの井戸の密度を計算します。
5. 井戸の密度の算出
これで、調査地域（帯水層の場合もあります）内の井戸が揃いましたので、井戸の密度を計算することができます。
1. プロセッシングツールボックスを開く：メインメニューのプロセッシング｜ツールボックスを選択します。
2. プロセッシングツールボックスでベクタ解析 | ポリゴン内の点の数を選択します。
3. ポリゴン内の点の数ダイアログでStampriet_boundaryをポリゴンとして、boreholes_clippedをポイントとして選択します。カウント属性名をNo of boreholesに変更します。カウントの出力をStampriet_Data GeoPackageにStampriet study areaというレイヤ名で保存します。
この出力の名前を使っているのは、このツールがStampriet_boundaryポリゴンのコピーを作成し、No of boreholesフィールドを追加しているからです。GeoPackageにはまだ調査地域の境界がなかったので、このように追加します。
4. 実行をクリックします。処理後に閉じるをクリックします。
Note that if your input has multiple polygons (e.g. different aquifers) it will produce the count per polygon. You can apply the next steps then to calculate the borehole density for all polygons.
5. Stampriet_boundaryのスタイルを Stampriet study areaにコピーします。 レイヤパネルでStampriet_boundaryを右クリックし、スタイル | スタイルのコピー | 全スタイルカテゴリを選択します。さらにStampriet study areaを右クリックしスタイル | スタイルの貼り付け | 全スタイルカテゴリを選択します。
6. レイヤパネルからStampriet_boundaryを削除します。
それでは、Stampriet study areaの属性テーブルを見てみましょう。
7. レイヤパネルで、Stampriet study areaを右クリックし、属性テーブルを開くを選択します。
調査対象地域には64の井戸があることがわかります。
次に、調査対象地域のポリゴンの面積を計算する必要があります。
8. 属性テーブルで 
ボタンをクリックし編集モードに切り替えます。
9. 
をクリックしてフィールド計算機を開きます。
10. フィールド計算機ダイアログで、出力する属性（フィールド）の名前にArea (km2) 、フィールド型に小数点付き数値（real）を選択し新しいフィールドを作成します。
11. 式の下で中央のパネルからジオメトリ | $areaを選択し、$areaをダブルクリックして式に追加します。AREA は、多角形の面積を地図の単位で計算します。ここではm2とします。次のようにして式を完成させます。
$area / 1000000
ボタンを使って割り算を追加します。
下の図のようになります。
次に、井戸の密度を計算します。
13. をクリックしてフィールド計算機を開きます。
14. フィールド計算機ダイアログで、出力する属性（フィールド）の名前にBorehole density per km2 、フィールド型に小数点付き数値（real）を選択し新しいフィールドを作成します。以下の式を使用します。
"No of boreholes" / "Area (km2)"
フィールドを追加するには、中央パネルのフィールドと値にあるフィールド名をダブルクリックしてください。 ボタンを使って割り算を追加します。
16. をクリックして編集モードをオフにしポップアップで保存をクリックして結果を保存します。
計算では、調査対象地域の井戸の密度は0.2/km2になりました。
次のセクションでは、プロジェクトにDEMを追加します。
6. SRTM 1-Arc Second DEMのダウンロード
このセクションでは、プロジェクトにDEMを追加していきます。ここでは、Shuttle Radar Topography Mission (SRTM) 1 Arc Secondプロダクトを使用します。これはフリーのDEMで、南緯56度から北緯60度をカバーし、赤道上では約30mの空間分解能を持ちます。
SRTMのタイルは、USGS Earth Explorerのウェブサイトからダウンロードするか、プラグインを使ってダウンロードすることができます。ここでは、QGISのSRTM Downloaderプラグインを使用します。
1.メインメニューでプラグイン | プラグインの管理とインストールを選択します。
2. SRTM-Downloader プラグインを検索し、プラグインをインストールしたらプロジェクトに戻ります。
3. ツールバーでこのボタン 
をクリックしてSRTM-Downloaderダイアログを開きます。
4. SRTM-Downloaderのダイアログで、Set canvas extentをクリックします（調査地域の境界にズームしていることを確認してください）。Output-pathには、このプロジェクトのデータが保存されているフォルダを参照します。
5. Downloadをクリックします。
6. 初回は、ログイン情報を入力するポップアップが表示されます。SRTMタイルをダウンロードするには、 https://urs.earthdata.nasa.gov/users/new でアカウントを作成する必要があります。
7. ユーザー名とパスワードを入力します。Save Credentialsをクリックし、OKをクリックします。
これで、調査対象地域をカバーするSRTMタイルのダウンロードが始まります。
8. ポップアップでOKをクリックします。
9. ダイアログを閉じます。
7. DEMを切り取って再投影
SRTMタイルが地図キャンバスに追加されています。S25E018.hgtレイヤーにズームすると、より広い範囲をカバーしていることがわかります。また、投影は地理座標系（GCS）で、緯度・経度座標は度単位です（EPSG:4326）。
このため、DEMを切り取って再投影し、GeoPackageに追加する必要があります。
1. レイヤパネルで S25E018.hgt を右クリックしエクスポート | 名前をつけて保存...を選択します。
8. DEMのスタイリング
をクリックしてレイヤスタイルパネルを開きます。
これで、プロジェクトは下の図のようになります。
次のセクションでは、DEMの標高を井戸の属性テーブルに追加していきます。
9. 井戸にDEMの標高を付与する
これから、井戸の位置でDEMの標高を付与していきます。
そのためには、Point sampling toolプラグインをインストールする必要があります。なお、プラグインのインストールにはインターネット接続が必要です。
1. メインメニューでプラグイン | プラグインの管理とインストールを開きます。
2. プラグインダイアログでPoint sampling tool を検索しインストールします。
3. ツールバーで 
ボタンをクリックしPoint Sampling Toolダイアログを開きます。
4. In the Point Sampling Tool ダイアログのGeneralタブでboreholes_clippedをLayer containing sampling points.として選択します。boreholes_clippedとDEMレイヤからの全てのフィールドを選択します。なお、このリストに表示するには、レイヤパネルでこれらのレイヤーをチェックする必要があります。このツールは既存のGeoPackageに追加することができないので、Output point vector layerをESRI Shapefileに保存します。名前は Boreholes_Z.shp とします。
5. 次にFieldsタブをクリックします。そこでは、必要に応じて出力フィールドの名前を変更することができます。ここでは、DEMフィールドをZにだけ変更します。
6. OKをクリックします。処理後にダイアログを閉じます。
7. Boreholes_Z レイヤの属性テーブルをチェックします。
Z フィールドが追加されていることがわかります。
次のセクションでは、標高とDEMからのZ値の散布図を作成してみます。
10. 報告された標高とDEMの値を比較
Boreholes_Z属性テーブルには、報告された標高とDEMからの値があります。この2つのフィールドを散布図にして比較してみましょう。
elevationフィールドには、多くのデータなしの値が含まれています。散布図ではこれらを考慮しないことにしましょう。
1. Boreholes_Zの属性テーブルを開きます。
2. 式による地物選択 ボタン 
をクリックします。
3. 式を書きます。
"elevation" < 0
4. 地物を選択をクリックし、ダイアログを閉じます。
10個の地物が選択されました。編集モードに切り替えて、ゴミ箱をクリックすれば、属性テーブルから削除することができます。しかし、散布図では、選択を反転させて、選択されたフィーチャーだけを使用するようにできます。
5. 
をクリックして選択を反転させます。
散布図を作成するには、Data Plotlyプラグインをインストールする必要があります。
6. プラグインマネージャーからData Plotlyプラグインをインストールします。
7. ツールバーで 
アイコンをクリックしてData Plotly パネルを開きます。
8. Data Plotly パネルでScatter PlotをPlot typeとしてBoreholes_Z.shpをLayerとして選択します。Use only selected featuresのチェックを入れます。elevation をX fieldに、ZをY fieldにします。その他はデフォルトのままにしておきます。
9. 
をクリックして凡例設定タブにいきます。
10. Show legend のチェックを外します。Plot title を Borehole vs DEM elevationに変更します。 Legend titleを削除します。X label にBorehole elevation (m)を、Y labelにDEM elevation (m) を入力します。残りはデフォルトのままにします。
11. パネルの下のCreate Plotをクリックします。
プロットを確認します。
次のセクションでは、標高の欠損値をDEMの標高に置き換えてみます。
11. nodataの値をDEMの標高に置き換える
散布図が妥当に見えるので、Boreholes_Z.shpレイヤーのelevationフィールドのnodata値をDEMからの値に置き換えてみます。
1. Boreholes_Z.shp レイヤから属性テーブルを開きます。
2. 
をクリックして選択を全て解除します。
3. 
をクリックして編集モードに切り替えます。
4. 出力をelevationに変更し、
ボタンをクリックすると式が作成されます。
5. 式ダイアログで次のような式を作ります。
if( "elevation" < 0, "Z" , "elevation" )
これは、標高フィールドの値が0より小さい場合はZフィールドの値を使用し、そうでない場合はelevationフィールドの値を使用するという意味になります。
6. OKをクリックします。
これで、式は属性テーブルで次のように表示されるはずです。
7. 全て更新をクリックします。
これで、すべての欠損値がDEMからの標高に置き換えられました。
8. をクリックし編集モードをオフにし保存をクリックして結果を保存します。
次のセクションでは、井戸調査で得られた地下水の絶対量を計算します。
12. 井戸の地下水位の算出
井戸の属性テーブルには、地下水位までの深さのデータがあります。
このセクションでは、井戸の水位の絶対標高を計算します。
1. Boreholes_Z.shpの属性テーブルを開きます。
2. 
をクリックして編集モードに切り替えます。
まず、depth_mフィールドの深さのデータが欠損しているフィーチャーを削除します。
3. 
をクリックしてQGIS式による選択ダイアログを開きます。
4. 次のような式を作ります。
"depth_m" < 0
5. 地物の選択をクリックします。閉じるをクリックしてダイアログを閉じます。
6. 
をクリックして選択された地物を削除します。
7. 
をクリックして今までの編集内容を保存します。
8. 
をクリックしてフィールド計算機ダイアログを開きます。
9. フィールド計算機ダイアログで出力する属性（フィールド）の名前をGW_level_mにし、フィールド型を小数点付き数値（real）にして新しいフィールドを作成します。その他の設定はデフォルトのままで、以下の式を入力します。
"elevation" - "depth_m"
10. OKをクリックします。
11. 結果を確認し、 をクリックして編集モードをオフにします。保存をクリックして結果を保存します。
13. 井戸の地下水位をラスタに補間する
井戸の地下水位がわかったところで、そのポイントをラスタに補間します。
ここでは2つの方法を使います。ここでは、IDW（Inverse Distance Weighting）とThiessenポリゴンの2つの方法を使います。
まず、IDWによる補間を行ってみましょう。
1. プロセッシングツールボックスで内挿 | IDW内挿（逆距離荷重法）を選択します。
2. IDW内挿（逆距離荷重法）ダイアログでBoreholes_Z.shpを入力ベクタとして選択します。GW_level_mを内挿対象の属性（フィールド）として使用します。
をクリックしてこれらの設定をリストに追加します。距離係数（p）を2のままにしておきます。Stampriet study area レイヤーから領域を選択します。ピクセルサイズを30mに変更し、補完したラスタをGW_level_IDW.tifとしてプロジェクトフォルダに保存します。
3. 実行をクリックします。処理後に閉じるをクリックします。
その結果をカラーランプでスタイリングしてみましょう。
4. レイヤパネルでGW_level_IDWレイヤを選択し 
をクリックしてレイヤスタイルパネルを開きます。
5. レイヤスタイルパネルでレンダラーを単バンド擬似カラーに変更し、全カラーランプ | RdBuを選択します。
レイヤースタイルパネルは以下のようになります。
そして、補間されたレイヤーは下の図のようになります。視覚化を向上させるために、さまざまなランプやクラスを試すことができます。IDW補間は滑らかな出力をしていますが、ボーリングの数が少ないため、個々のボーリングや外れ値に多くのバイアスがかかっていることがわかります。
それでは、Thiessenポリゴン（別名：nearest neighbourまたはVoronoi tesselation）法を適用してみましょう。
6. メインメニューでラスタ | 解析 | グリッド (最近傍法)を選択します。
7. グリッド (最近傍法)ダイアログでBoreholes_Z.shpレイヤを入力レイヤ（点）として選択します。詳細パラメータの下で GW_level_m を 内挿するZ値の属性として選択します。内挿（最近傍法）した結果をGW_level_Thiessen.tifに保存します。その他はデフォルトのままにします。
8. 実行をクリックします。閉じるをクリックしてダイアログを閉じます。
9. 結果をスタイルします。
このアルゴリズムは研究領域全体をカバーしておらず、入力ポイントレイヤのいわゆる凸包に限定されています。このアルゴリズムは、IDW内挿法よりも滑らかではない結果をもたらします。しかし、データが少ない場合には、最も近い値が最良の近似値となる可能性があります。
