とらりもん - 数値標高モデル (DEM) を用いた流路データの作成 Diff

2012/02/06 筑波大学 田中健太郎

!はじめに
本章では、数値標高モデル (Digital Elevation Mode) を用いて、流路データを作成し, 各流路の長さや流路ごとの流域面積を推定する。


!データの取得

これから扱うのは、ASTERのGDEMデータである。
※2011年09月30日よりバージョン2が公開された。

データを取得するためには、下記のサイトに登録する必要がある。

http://www.gdem.aster.ersdac.or.jp/

* 追記: 2014年7月現在、上記のサイトはすでに閉鎖されていた。代わりに[[こちら|http://gdem.ersdac.jspacesystems.or.jp/index.jsp]]のサイトに登録する。
* 再追記:2018年2月現在では、追記で紹介されているページもリンクが切れている模様。(2018/02/12 菊島)

会員登録後 (無料)に"Login"し, "Search" を選択.

すると下記の画面が表示されるので, 赤枠の"Start" をクリックし, データの欲しい範囲をクリックする.

{{attach_view(gdem1.png)}}


今回注文するのは,
*ASTGTM2_N35E138
*ASTGTM2_N35E137
*ASTGTM2_N36E137
*ASTGTM2_N36E138
とする. (Nは画像左下の緯度, Eは経度を示す.)

{{attach_view(gdem2.png)}}



"Next" ボタンをクリックすると, 確認画面が表示されるので, 注文したい範囲と一致していれば "Next" をクリックする.

使用目的を尋ねられるので, 選択し, "Agree" をクリック.


ファイルを個別にダウンロードするか, 一括ダウンロードするかを尋ねられる.


{{attach_view(gdem3.png)}}



個別の場合は, ファイル名右横の"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

{{attach_view(gdem4.png)}}


これでASTGTM2_N35E137,ASTGTM2_N35E138,ASTGTM2_N36E137,ASTGTM2_N36E138の4枚の画像が結合され, 1枚の画像が作成されていることがわかる.
[[以前|http://pen.envr.tsukuba.ac.jp/~RStiger/hiki2/?%A5%E2%A5%B6%A5%A4%A5%AF%BD%E8%CD%FD]]|モザイク処理]]においても, モザイク処理を行った.
ここでは, 以前とは別の方法として, 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

{{attach_view(gdem5.png)}}

ベクター化が完了し, ベクターデータとしての河道網データが作成されていることがわかる.
次に流路ごとの流域マップに, 河道網データを重ねて表示してみる.

d.erase
d.rast basin
d.vect stream

{{attach_view(gdem6.png)}}

!流路の長さおよび流域面積の出力

流路データの解析には, 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 筑波大学 牧野 悠