GRASSのインストール、標高データを用いた地滑り危険度マップの作成

ここでは、GRASSというGISソフトウェアを使った簡単な演習を説明しています。GRASSではリモートセンシング画像も表示、処理できます。このページは、GRASSのGUIでの操作を説明していきます。

[GRASSのインストール、標高データを用いた地滑り危険度マップの作成]
[植生指数 (NDVI) の計算、表示]
[標高データ (SRTM)の表示、植生指数 (NDVI) の3次元表示]
[反射率、輝度温度、標高データを用いた土地被覆分類]
[QGIS, Rを用いた公示地価データの空間統計分析：空間的自己回帰モデル]
[QGIS, Rを用いた公示地価データの空間統計分析：静的な時空間モデリング]

2016/04/21追記：演習1の参考となるWebページ

GRASS GIS 7 クィックスタート

GRASSのインストール

QGISとセットになったGRASSのインストーラーをダウンロードします。

1. インストーラーのダウンロード：download.osgeo.org/osgeo4w/osgeo4w-setup.exeより
2. osgeo4w-setup.exeの起動：インストールに必要なファイルのダウンロード
   osgeo4w-setup.exeをダブルクリック　→　「Download Without Installing」　→　保存先の指定　→　ミラーサイトの指定
3. （ダウンロード終了後）osgeo4w-setup.exeの起動：インストール
   osgeo4w-setup.exeをダブルクリック　→　「Install from Local Directory」　→　インストールディレクトリの指定　→　ダウンロードしたファイルのディレクトリの指定　→　インストールするパッケージの指定（分からなければ「Install」を選択して、全てをインストールする）

GRASSの起動と設定

GRASSのアイコンをダブルクリックすれば起動します。

1. GRASS GIS database directoryの指定：

   データを保存するディレトクリ（フォルダ）を指定。任意の場所で構わないが、スペースなしの英語名のディレトクリを指定する必要がある。一例を下記に示す。

   C:\GRASS_data
   

2. データベースとロケーションの指定：

   GISデータディレクトリ:   C:\GRASS_data
   プロジェクトロケーション:DEM
   Location Title:          PERMANENT
   

3. 新しいロケーションの作成方法を選択
   • EPSG (European Petroleum Survey Group) コードを指定
   • EPSGコード：32653（WGS84楕円体、UTM投影法、ゾーン番号53）

「ロケーション<DEM>が作成されました　デフォルト領域と解像度をセットしますか？」と尋ねられるが、「いいえ」ボタンを押す。その後、「マップセットを作成するか？」とも尋ねられるが、こちらもキャンセルする。すると、GRASSの画面が表示される。

標高データの入手

1. 国土地理院 基盤地図情報ダウンロードサービス
   http://fgd.gsi.go.jp/download/からログイン（登録が必要）
2. 「基盤地図情報数値標高モデル：JPGIS(GML)形式」→「10mメッシュ」を選択。下記の番号のいずれかを選択し、データ（zipファイル）をダウンロードする（班によって異なる）。
   • 523534
   • 523535
   • 523536
   • 523524
   • 523525
   • 523526

メッシュ番号を選択する画面

参考：

• 標準地域メッシュシステム
  • 1次メッシュ：40'×1°
  • 2次メッシュ：5'×7'30"
• 数値標高モデル10 mメッシュ
  • 1/25,000地形図（2次メッシュ）を経度方向に1,125等分、緯度方向に750等分して得られる各格子（0.4"四方）の中心の標高
  • 地理情報標準プロファイル (Japan Profile for Geographic Information Standards: JPGIS) に準拠
  • Geography Markup Language (GML) で記述

QGISを使ったGeotiff画像への変換

上記の数値標高モデル10 mメッシュデータをGRASSで使えるようにする方法の一つが、GeoTIFF画像に変換する方法です。

1. fgddemImporter (QGISプラグイン)のダウンロード
   https://dl.dropboxusercontent.com/u/21526091/qgis-plugins/plugins.xml
2. QGIS上での変換：詳細は 「fgddemImporter プラグイン」 というWebページに記載されています。
   1. fgddemImporterフォルダをQGISのプラグインフォルダ(ユーザフォルダ/.qgis2/python/plugins)にコピー。
   2. QGIS Desktop（QGIS Browserではない）を起動し、「プラグイン」→「プラグインの管理とインストール」→「fgddemImporter」にチェックを入れる。
   3. 「プラグイン」→「基盤情報DEMインポータ」でダウンロードしたZIPファイルを指定し、インポートする。この結果、GeoTIFF画像が生成される。

QGIS上で標高データをインポートして、GeoTIFF画像が生成された様子。

3. GeoTIFF画像の移動
   GRASSのデータ用のフォルダへ移動します。上記の例では、下記のフォルダへ移動します。

   C:\GRASS_data\DEM\PERMANENT
   

注意：523535のメッシュを選択すると、この時点で「FG-GML-5235-35-DEM10B.tif」というGeoTIFF画像が生成されています。ハイフンがついた名前だとGRASSで処理できませんので、「dem.tif」等、ハイフンがつかない分かりやすい名前に変更して下さい。以下では、「dem.tif」という名前に変更されたという前提で話を進めていきます。

標高画像ファイルのインポートと表示

1. 画像ファイルのインポート：
   • ファイル　→　ラスターデータのインポート　→　Common formats import [r.in.gdal]
   • マップ名：dem

   注意：既に「dem」という名前でインポートしている場合には、「既存のファイルに上書きする」をチェックしましょう。

2. 範囲の設定：画像ファイルの範囲をデフォルト領域に設定する
   • 設定　→　領域　→　領域設定 [g.region]
   [入力項目]
   

   • 「既存領域」→ Set region to match raster map(s): dem@PERMANENT
   • 「出力」→ 「現在の領域を表示」をチェック

   注意：この範囲設定を正しく行わないと、領域の範囲が行方向（column）に1画素、列方向（row）に1画素しか存在しない設定になっているため、以下のエッジ抽出等が実行できません。

3. カラーテーブルの設定
   • ラスター　→　カラー調整　→　カラーテーブル [r.colors]
   [入力項目]
   

   • 「Map」→ Name of raster map(s): dem@PERMANENT
   • 「Define」→ Name of color table: elevation

この処理が終わると、この図が表示されます。

地形解析

1. エッジ抽出（水平方向に続くエッジ）：
   • edgeH.txtという名前で、以下の内容を記述したテキストファイルを作成

   TITLE 3x3 edgeH
   MATRIX 3
   1 2 1
   0 0 0
   -1 -2 -1
   DIVISOR 1
   TYPE P
   

   補足：本来はDIVISOR=8でなければならないが、エッジのコントラストを出すために1にしている。上記の演算の意味は講義中に説明する。

   • 画像　→　画像フィルター　→　行列/畳み込みフィルター [r.mfilter]
   [入力項目]
   
   • 「必須」→ 入力するラスターマップ名: dem@PERMANENT
   • 「必須」→ 入力するラスターマップ名: dem_edgeH
   • 「必須」→ フォントファイルへのパス: C:\GRASS_data\DEM\PERMANENT\edgeH.txt
   • カラーテーブルの設定：
     ラスター　→　カラー調整　→　カラーテーブル [r.colors]
     「terrain」のカラーテーブルを選択。

処理後の画像

2. エッジ抽出（鉛直方向に続くエッジ）：
   • edgeV.txtという名前で、以下の内容を記述したテキストファイルを作成

   TITLE 3x3 edgeV
   MATRIX 3
   -1 0 1
   -2 0 2
   -1 0 1
   DIVISOR 1
   TYPE P
   

   • 画像　→　画像フィルター　→　行列/畳み込みフィルター [r.mfilter]
   [入力項目]
   
   • 「必須」→ 入力するラスターマップ名: dem@PERMANENT
   • 「必須」→ 入力するラスターマップ名: dem_edgeV
   • 「必須」→ フォントファイルへのパス: C:\GRASS_data\DEM\PERMANENT\edgeV.txt
   • カラーテーブルの設定：
     ラスター　→　カラー調整　→　カラーテーブル [r.colors]
     「terrain」のカラーテーブルを選択。

処理後の画像

3. 勾配の計算1：上記のedgeHとedgeVの両者から勾配を計算
   • ラスター　→　ラスターマップカリキュレータ [r.mapcalc]
   [入力項目]
   
   • 「出力」: dem_sobel
   • 「式」: int(sqrt(dem_edgeH@PERMANENT^2 + dem_edgeV@PERMANENT^2 ))
   • カラーテーブルの設定：
     ラスター　→　カラー調整　→　カラーテーブル [r.colors]
     「grey.eq」のカラーテーブルを選択。

処理後の画像

4. 勾配の計算2：GRASSの機能を利用
   • ラスター　→　地形解析 → 傾斜と方位 [r.slope.aspect]
   [入力項目]
   
   • 「必須」→ Name of input elevation raster map: dem@PERMANENT
   • 「出力」→ 出力する傾斜ラスターマップ名: slope
   • 「出力」→ 出力する傾斜方位ラスターマップ名: aspect

処理後の画像（勾配）, （方位）

5. 陰影起伏図
   • ラスター　→　地形解析 → Compute shaded relief [r.relief]
   [入力項目]
   
   • 「必須」→ 入力するラスターマップ名: dem@PERMANENT
   • 「必須」→ Name for output shaded relief map: shade45
   • 「Sun position」→ 地平線上の太陽高度の角度: 45
   • 「Sun position」→ 太陽の方位（東から反時計回り）: 45

処理後の画像
[太陽高度45°、東から45度の方向角）][（太陽高度45°、東から90度の方向角）]

流域地形解析

1. 流下方向の計算：最も急な下り勾配の方向を決定。1〜8の値が生成され、45を掛けると、東方向から反時計回りの方向角になる（0の場合は勾配なし）。
2. 累積ピクセル数の計算：当該ピクセルに流れ込むピクセル数を計算。
3. stream networkの抽出：累積ピクセル数が一定の閾値より大きいピクセルを取り出して、stream networkを生成
4. basinの抽出：各streamセグメントに対し、basinを抽出し、番号を付与する。
（上記の4つを同時に処理）
• ラスター　→　水文モデリング　→　分水界解析 [r.watershed]
[入力項目]

• 「Input」→ Name of input elevation raster map: dem@PERMANENT
• 「Input」→ Minimum size of exterior watershed basin: 100
• 「出力」→ Name for output accumulation raster map: accumulation
• 「出力」→ Name for output drainage direction raster map: drainage
• 「出力」→ Name for output basins raster map: basin
• 「出力」→ Name for output stream segments raster map: stream

処理後の画像
[流下方向] [累積ピクセル数] [stream network] [basin]

5. stream networkのラスター画像をベクターファイルに変換
   • ラスター　→　ラスターマップの作成 → Null値の管理 [r.null]
   [入力項目]
   
   • 「必須」→ Name of raster map for which to edit null values: stream@PERMANENT
   • 「変更」→ [複数] NULL値をセットするセルの一覧：0
   • ラスター　→　フィーチャーの変形 → 縮小 [r.thin]
   [入力項目]
   
   • 「必須」→ 入力するラスターマップ名：stream@PERMANENT
   • 「必須」→ 出力するラスターマップ名：stream_thin
6. ラスター　→　マップタイプの変換 → 「ラスター→ベクター」 [r.to.vect]
[入力項目]

• 「必須」→ 入力するラスターマップ名：stream_thin@PERMANENT
• 「必須」→ 出力するベクターマップ名：stream_vect
• 「必須」→ Output feature type: line

処理後の画像
[stream network（ベクター）]
[stream network（ベクター）をDEM画像に重ね合わせた結果（「マップレイヤー」のタグで表示させるデータを選択して重ね合わせられる）]

演習課題

指定日に各班の代表が発表します。それまでに各班で共同で作業をするように。

1. 各班で異なる標高データ（2次メッシュ）をダウンロードし、河川ネットワークを作成し、陰影起伏図の上に描画せよ。
2. 斜面勾配が40°以上となる地域を抽出して、陰影起伏図の上に描画せよ。
3. 降水量や河川の流量データ等、水文に関係するデータと、今回の学習で得られた地形解析結果を関連づけて解析したり、あるいはGRASSに実装されている他の水文関連の関数を使った解析したりし、その成果を報告せよ。

カラーテーブルを変更したDEM画像	水平方向に続くエッジの抽出結果
鉛直方向に続くエッジの抽出結果	水平・鉛直方向のエッジから計算した勾配
r.slope.aspectで求めた勾配	r.slope.aspectで求めた方位
陰影起伏図 （太陽高度45°、南東方向）	陰影起伏図 （太陽高度45°、南方向）
流下方向	累積ピクセル数
stream network	basin
stream network（ベクター）	stream network（ベクター）を DEM画像に重ね合わせた結果

[GRASSのインストール、標高データを用いた地滑り危険度マップの作成]
[植生指数 (NDVI) の計算、表示]
[標高データ (SRTM)の表示、植生指数 (NDVI) の3次元表示]
[反射率、輝度温度、標高データを用いた土地被覆分類]
[QGIS, Rを用いた公示地価データの空間統計分析：空間的自己回帰モデル]
[QGIS, Rを用いた公示地価データの空間統計分析：静的な時空間モデリング]

[講義「リモートセンシングと地理情報システム(GIS)」のトップページ]

[須崎純一のページへ]