数値標高モデル (DEM) を用いた流路データの作成
2012/02/06 筑波大学 田中健太郎
データの取得
これから扱うのは、ASTERのGDEMデータである。 ※2011年09月30日よりバージョン2が公開された。
データを取得するためには、下記のサイトに登録する必要がある。
http://www.gdem.aster.ersdac.or.jp/
- 追記: 2014年7月現在、上記のサイトはすでに閉鎖されていた。代わりにこちらのサイトに登録する。
- 再追記:2018年2月現在では、追記で紹介されているページもリンクが切れている模様。(2018/02/12 菊島)
会員登録後 (無料)に"Login"し, "Search" を選択.
すると下記の画面が表示されるので, 赤枠の"Start" をクリックし, データの欲しい範囲をクリックする.
今回注文するのは,
- ASTGTM2_N35E138
- ASTGTM2_N35E137
- ASTGTM2_N36E137
- ASTGTM2_N36E138
とする. (Nは画像左下の緯度, Eは経度を示す.)
"Next" ボタンをクリックすると, 確認画面が表示されるので, 注文したい範囲と一致していれば "Next" をクリックする.
使用目的を尋ねられるので, 選択し, "Agree" をクリック.
ファイルを個別にダウンロードするか, 一括ダウンロードするかを尋ねられる.
個別の場合は, ファイル名右横の"Download" をクリック.
一括ダウンロード場合は, ファイル名左横にチェックを入れ, "Download" をクリック.
後は待っていれば自動でダウンロードが始まる.
データの解凍および読み込み
unzip ASTGTM2_N35E138.zip unzip ASTGTM2_N35E137.zip unzip ASTGTM2_N36E137.zip unzip ASTGTM2_N36E138.zip cd ASTGTM2_N35E138 r.in.gdal in=ASTGTM2_N35E138_dem.tif out=ASTGTM2_N35E138 --o cd ../ASTGTM2_N35E137 r.in.gdal in=ASTGTM2_N35E137_dem.tif out=ASTGTM2_N35E137 --o cd ../ASTGTM2_N36E137 r.in.gdal in=ASTGTM2_N36E137_dem.tif out=ASTGTM2_N36E137 --o cd ../ASTGTM2_N36E138 r.in.gdal in=ASTGTM2_N36E138_dem.tif out=ASTGTM2_N36E138 --o
課題 読み込んだ画像を, ディスプレイに表示し, 正しく読み込まれているか確認せよ.
つぎに読み込んだ画像を結合し (モザイク処理), 1枚の画像にする.
g.region rast=ASTGTM2_N35E137 g.region n=37 s=35 w=137 e=139 r.series input=ASTGTM2_N35E137,ASTGTM2_N35E138,ASTGTM2_N36E137,ASTGTM2_N36E138 output=ASTGTM2_Tyubu method=maximum d.mon x0 r.colors ASTGTM2_Tyubu color=elevation #elevationカラーにセットする. d.rast ASTGTM2_Tyubu
これでASTGTM2_N35E137,ASTGTM2_N35E138,ASTGTM2_N36E137,ASTGTM2_N36E138の4枚の画像が結合され, 1枚の画像が作成されていることがわかる. 以前においても, モザイク処理を行った. ここでは, 以前とは別の方法として, r.series を使った.
流路および流路の抽出
作成したモザイクDEM画像から, r.watershedコマンドを用いて, 流路および各流路ごとの流域を抽出する.
r.watershed elevation=ASTGTM2_Tyubu basin=basin stream=stream drainage=drainage threshold=10000 ※thresholdは, 外側の分水域の最小セル数であり, 値が小さいほど, 細い河川も抽出対象となる.ここでは, 10000と設定した. ※※r.watershed で指定可能なオプションは, 以下の通り. elevation 入力:標高を示すファイル depression 入力:地盤沈下量を示すファイル flow 入力:1セルの流量を示すファイル disturbed.land 入力:開発地域の割合を示すファイル blocking 入力:曲面流でブロック化した地形を示すファイル accumulation 出力:流れ出るセル数を示すファイル drainage 出力:排出方向を示すファイル basin 出力:分水域のラベルファイル stream 出力:河川を示すファイル half.basin 出力:半盆地のラベルファイル visual 出力:結果を可視化したファイル length.slope 出力:傾斜長と傾斜勾配(LSファクター)を示すファイル slope.steepness 出力:RUSLE用の傾斜勾配(Sファクター)を示すファイル threshold 外側の分水域の最小セル数
r.watershedから作成された流路データ (stream) は, ラスター形式である. しかし一般的に流路データは, ベクターデータとして整備されている. そこで次のように, ラスターtoベクター変換を行う.
r.thin in=stream out=stream_thin -o ※r.thin は, 細線化を行うコマンドである. ライン形式のベクターデータを作成する前には, 必ず細線化処理を行う必要がある. r.to.vect input=stream_thin out=stream_thin d.vect stream
ベクター化が完了し, ベクターデータとしての河道網データが作成されていることがわかる. 次に流路ごとの流域マップに, 河道網データを重ねて表示してみる.
d.erase d.rast basin d.vect stream
流路の長さおよび流域面積の出力
流路データの解析には, r.stream.stats という非常に便利なコマンドが存在している.
※ただしr.stream.statsは, GRASSの拡張コマンドであり, 通常は使用できない. g.extensionコマンドを使用し, モジュールを追加することで, 使用可能となる.
r.stream.statsにより, 流域ごとに平均流路長や平均勾配といった情報を得ることができる. そのためstreamデータは, 対象にしたい流域 (複数可) でマスクをかける必要がある.
r.mapcalc "stream_mask=if(basin==4644,stream,null())" r.stream.stats stream=stream_mask dir=drainage dem=ASTGTM2_Tyubu
例)
Max order | Tot.N.str. | Tot.str.len. | Tot.area. | Dr.dens. | Str.freq. (num) | (num) | (km) | (km2) | (km/km2) | (num/km2) 4644 | 1 | 10.5492 | 29.2529 | 0.3606 | 0.0342 Stream ratios with standard deviations: Bif.rt. | Len.rt. | Area.rt. | Slo.rt. | Grd.rt. Order | Avg.len | Avg.ar | Avg.sl | Avg.grad. | Avg.el.dif num | (km) | (km2) | (m/m) | (m/m) | (m) 4644 | 10.5492 | 29.2529 | 0.0647 | 0.0248 | 262.0000 Order | Std.len | Std.ar | Std.sl | Std.grad. | Std.el.dif num | (km) | (km2) | (m/m) | (m/m) | (m) Order | N.streams | Tot.len (km) | Tot.area (km2) 4644 | 1 | 10.5492 | 29.2529 Order | Bif.rt. | Len.rt. | Area.rt. | Slo.rt. | Grd.rt. | d.dens. | str.freq. 4644 | 0.0000 | 0.0000 | -nan | 0.0000 | 0.0000 | 0.3606 | 0.0342
※実際はたくさんの行が出力されるが, 上記は必要な部分を抜粋した.
更新履歴
2014/07/03 筑波大学 牧野 悠
Keyword(s):
References:[GIS入門]