土地の傾斜区分毎の面積算出方法~QGISとGRASS GISを使う~ |
記事中で示した「傾斜区分ごとの面積」の算出の方法について、備忘録を兼ねてやり方をアップしました。
使用したソフトウェアはフリー・オープンのGISソフトウェア(Foss4G)です。
地図の表示とベクタデータでの作業=QGISver1.8.0
ラスタデータでの作業=GRASS GIS ver6.4.3
*いずれもUbuntu版
1.分析方法の概要
手順としては、
(1)傾斜DEMの作成
(2)傾斜区分ごとに二値画像化
(3)二値画像化したラスタのベクタ化(ポリゴン化)
(4)ベクタデータのエクスポート(shape file化)
(5)流域ベクタの作成
(6)流域ベクタで傾斜区分ベクタをクリップ
(7)流域ごとの傾斜区分ベクタの面積を算出
となります。
(1)~(4)をGRASS GIS(ver6.4.3)で、(5)~(7)はQGIS(ver1.8.0)を使用しました。
2.使用データ
分析に用いたデータは、国土地理院発行の基盤地図情報数値標高モデルを利用しました。
国土地理院の「基盤地図情報ダウンロードサービス」から必要なデータをダウンロードできます。
3.傾斜区分図の作成
まず、標高の代わりに傾斜角度の入ったDEMデータを作成します。
標高データから傾斜を計算させるわけです。
道南エリアのDEMデータ(GeoTiff化した数値標高モデル)はGRASSにインポート済みとします。
「Raster」→「地形解析」→「傾斜と方位」
「必須」タブでインポート済みのDEMデータを指定します。
ここでは「merge3_utm54」というデータを指定しています。
「出力」タブの「出力する傾斜ラスターマップ名」に出力する新規データ名を指定します。
「keisha」という安直なネーミングとしました。
余談ですが、私はファイル名に無理して英語を使うぐらいなら「日本語ベタ打ちの方がまし」と開き直っています。
「実行」をクリックすると「keisha」という名前の傾斜ラスタが作成されます。
4.傾斜区分ラスタの二値画像化
ここでやりたいことは、傾斜角度が1/100(0.5729°)未満の部分とそれ以外の部分にデータを区分することです。
テキストファイルに以下のようなコマンドを書きます。
0 thru 0.5729 =1
0.5729 thru 100 =NULL
0から0.5729の値をもつ部分に1を代入し、0.5729から100 の値を持つ部分にNULL値を代入する、という意味です。
上記のテキストファイルに適当な名前をつけて保存します。
ここでは「reclass_u0.5729」としました。
「Raster」→「カテゴリー値・ラベルの変更」→「Reclassify」
「必須」タブの「再分類されるラスターマップ」に先ほど作成した「keisha」を指定します。
「出力するラスターマップ名」に「reclass_u0.5729」としました。
「オプション」タブの「再分類ルールを含むファイル」で、先ほど保存したテキストファイル「reclass_u0.5729」を指定します。
「実行」をクリックすると、0.5729未満の部分は1、それ以外のところは値なしとなったラスタデータ「reclass_u0.5729」が作成されます。
5.二値画像化した傾斜ラスタのベクタ化(ポリゴン化)
二値画像化した傾斜ラスタをポリゴンデータに変換します。
「ファイル」→「マップタイプの変換」→「ラスター→ベクトル」
「必須」タブの「入力するラスターマップ名」に二値画像化したラスタデータ「reclass_u0.5729」を指定します。
「出力するベクトルマップ名」は「u05729」とします。
「地物フィーチャー」は「area」を選択します。
「実行」をクリックすると傾斜角度1/100以下の範囲がポリゴンベクタデータとして出力されます。
6.ファイルのエクスポート
GRASSからデータをエクスポートします。
「ファイル」→「ベクトルマップのエクスポート」→「Common export formats」
「必須」タブの「入力するベクトルマップ名」に先ほどベクタ化したファイルである「u05729」を指定します。
「OGR出力のデータソース名」も同じ名称の「u05729」としました。
「入力」タブの「地物の型」を「area」に指定します。
「実行」をクリックすると「u05729」という名称のshape fileが作成されます。
7.流域のポリゴンデータの作成
ここからは、QGISを使用して作業を進めます。
まず、流域ごとの範囲を示すポリゴンデータを作成します。
今回はおおむね市町村界を通るようなラインを目測で引いて作成しました。
「assabu_area」という名称のshape fileで保存します。
8.流域ベクタで傾斜ベクタをクリップ
流域ごとに傾斜ベクタを切り分けます。
GRASSからエクポートしたベクタデータ「u05729」をQGISで表示します。
紫色の部分が傾斜1/100未満の土地です。
海域も含まれていますので、当然このままでは使えません。
(国土地理院発行基盤地図情報数値標高モデル10mメッシュ6239~6241・6339~6341を改変して使用)
「ベクタ」→「空間演算ツール」→「クリップ」
「ベクタレイヤーの入力」にクリップされる傾斜ベクタ「u5729」を選択。
「レイヤをクリップする」にクリップする流域ベクタ「assabu_area」を選択。
「出力Shapefile」を「u05729_assabu」として、新たに生成するファイル名を指定します。
「OK」をクリックすると厚沢部川流域でクリップされた新たな傾斜ベクタ(紫色の部分)が作成されます。
(国土地理院発行基盤地図情報数値標高モデル10mメッシュ6239~6241・6339~6341を改変して使用)
9.面積フィールドの作成
クリップされた傾斜ベクタ「u05729_assabu」に面積のフィールドを追加します。
作成したベクタデータ「u05729_assabu」を表示します。
「レイヤ」→「属性テーブルのオプション」
下に並んでいるアイコンから、鉛筆風のアイコンをクリックして属性テーブルを編集可能にします。
アイコン列右端の計算機風アイコンをクリックして「フィールド計算機」を開きます。
「新しいフィールドを作る」にチェックします。
「出力フィールド名」を「area2」とします。
「関数リスト」→「ジオメトリ」→「$area」を選択します。
「式」フィールドには単に「$area」と表示されます。
「OK」をクリックすると新たなフィールド「area2」が追加され、行(地物)ごとの面積が表示されます。
10.流域ごとの面積の算出
「ベクタ」→「解析ツール」→「基本統計」
「ベクタレイヤの入力」に「u05729_assabu」を指定します。
「ターゲットフィールド」に先ほど新規作成した「area2」を指定します。
「OK」をクリックすると「統計値出力」に各種統計量が表示されます。
「統計値出力」のどこでも良いのでクリックして、コピー(Ctrl+C)します。
表計算ソフトなどにペーストします。
全ての統計量がペーストされるので、目的とする統計量以外は削除します。
ここでは、「合計」の「27940935.0」だけを残します。
以上の操作を他の流域でも行い、流域ごとの、傾斜区分ごとの面積一覧表を作成します。
追記
今回紹介した面積算出の手法は、2013年7月6日に行われたFOSS4GHokkaido2013ハンズオンデイ「QGIS中・上級編 ラスタ解析」(講師 北海道大学農学研究院の三島哲雄さん)で習得したラスタ解析の手法を応用しました。
三島さんとFOSS4G2013Hokkaidoスタッフのみなさまに感謝です。