とらりもん - 土地被覆分類図の作成 Diff
- Added parts are displayed like this.
- Deleted parts are displayed
like this.
2012/02/07 筑波大学 田中健太郎
!1. はじめに
[[前回|植生分類図]], 環境省自然環境保全基礎調査によって作成された植生分類図を扱った. 自然環境保全基礎調査では, 衛星データや空中写真, 現地調査の情報を元に作成されいる. 本章では, 実際にこうした土地被覆分類図を衛星データを用いて作成する方法を紹介する.
!2. 教師つき分類と教師なし分類
土地被覆分類方法の選択肢として, まず教師つき (supervised) と教師なし (unsupervised) の2つの方法がある.
教師つき分類とは, あらかじめ土地被覆分類項目がわかった訓練領域 (=教師データ) をヒントに, 未知の領域の分類を行う方法である.
そのため事前に訓練領域と呼ばれる地域を設定し, 訓練領域内の土地被覆項目 (ex)都市, 森林, 水田など) を把握し, デジタル化しておく必要がある. 一方で教師なし分類では, 画像値の類似性等を機械的に分析し, 分類を行う. そのため教師なし分類によって作成された分類項目に具体的な意味はなく, 地上検証などを用いて, 分類項目と実際の土地被覆項目 (ex)都市, 森林, 水田など)を対応させる必要がある.
つまり,
*教師付き分類は, ラベル付きデータ (事前に分類項目のわかったデータ) を用いて行う分類である.
*教師なし分類は, ラベル無しデータ (事前に分類項目のわからないデータ) を用いて行う分類である.
!3. 教師なし分類
本章では, まず教師なしの土地被覆分類について学ぶ. 教師なし分類の方法についても, ディシジョンツリー法, 最尤法, k-mean法, ISODATA法など, 様々な方法が存在する.
今回は, GRASSにおいて実装されている最尤法による分類を行う. 最尤法による分類では, (複数の) 衛星データのピクセル値 (画素値) が近いものを同じ項目として分類していく手順を,設定した項目数以下になるまで繰り返していく方法である.
用いる衛星データは, [[高分解能衛星Landsat: 植生指標NDVI, 温度画像]]で使用したデータと同じものとする.
GRASSを起動し,
$ grass
下記を選択する.
LOCATION: UTM54
MAPSET: PERMANENT
GRASSにおいて教師なし分類は, i.group, i.cluster, i.maxlikの3種類のコマンドで実装されている.
まずi.groupによって, 分類に使用する衛星データをグループ化する必要がある.
今回は, LandsatのBand1~4のデータを使用する.
for i in 1 2 3 4
do
r.null B$i setnull=0 #事前に0をnullにする.
done
g.region rast=B1
i.group group=landsat sub=main input=B1,B2,B3,B4
groupは, 衛星データをグループ化した際のグループ名である. subgroupを指定することによって, グループ内にも複数のグループを作成できる. 上記コマンド例では, グループ名 : Landsat, サブグループ名 : main とした.
Landsat ---- main ----- B1
|
-- B2
|
-- B3
|
-- B4
次にi.clusterとグループ化したデータを用いて, ピクセル値の類似性を分析し, 分類器 (分類のルール) を作成する.
i.cluster group=landsat subgroup=main sigfile=sig classes=10 report=class_report.txt
sigfileには, 分類器の情報が出力される.(追記:GRASS7ではsignaturefileとする.) classesは, 分類項目数であり, ここでは10とした. reportには, 分類 (クラスラリング ) の質を評価した情報が出力される.
最後にi.maxlikコマンドとsigfileを用いて, 対象地域全体の分類を行う.
i.maxlik group=landsat subgroup=main sigfile=sig class=LC_map_uns reject=LC_map_uns_re
classは, 分類結果のラスターデータ名である. rejectは, 分類の信頼水準を出力するラスターデータ名である.
※classと打ち込んでも動くが警告文が出てくるので、output=LC_map_unsと書いたほうが良い。(追記:2018/02/12 菊島)
さっそく分類結果を表示してみよう.
d.rast.leg LC_map_uns #d.rast.leg は, d.rastとd.legend を組みあわせた素晴らしいコマンドである. 消すときは, d.erase -f
{{attach_view(class1.png)}}
現時点で分類項目に意味はないが, 対象地域が10項目に分類されていることが確認できる.
続いて, 各分類項目と土地被覆項目とを対応させる. 方法は, 現地調査や空中写真判読など様々な方法があるが, 手軽に試すには衛星データのトゥルーカラー画像が良いだろう. まず対象地域が広域なので, 一部地域にズームする.
g.region n=4002967.5 s=3982476 w=412338 e=434511
d.rgb r=B3 g=B2 b=B1
d.mon x1
d.rast LC_map_uns
{{attach_view(class2.png)}}
{{attach_view(class3.png)}}
{{attach_view(class4.png)}}
2枚の画像を見比べてみると, 分類結果の項目番号が
1 水域
2 水域
3 植生
4 植生
5 植生
6 植生
7 都市
8 植生
9 都市
10 裸地
の様に対応していることが読みとれる. (ざっくりとした判読なので正しいとは限らないが, 話を進める.)
次に新しい項目番号を次の様に定義し, 分類結果の項目番号を対応させよう.
古い分類番号 新しい分類番号
1, 2 ---> 11 (水域)
3, 4, 5, 6, 8 ---> 12 (植生)
7, 9 ---> 13 (都市)
10 ---> 14 (裸地)
上記の対応を分類結果であるラスターデータに対して行うには, r.reclassを使う.
cat << EOF > reclass.txt
1 2 = 11 water
3 4 5 6 8 = 12 vegetation
7 9 = 13 urban area
10 = 14 bare
EOF
r.reclass input=LC_map_uns out=LC_map_uns_reclass rules=reclass.txt title=LC_map_uns_reclass
最後に作成した土地被覆分類図の色合いを設定し, 表示してみよう.
cat << EOF > color.txt
11 0:0:255
12 0:100:0
13 255:0:0
14 200:200:0
EOF
r.colors LC_map_uns_reclass color=rules < color.txt
g.region rast=LC_map_uns_reclass
d.rast.leg LC_map_uns_reclass
{{attach_view(class5.png)}}
以上がGRASSによる教師なし分類方法の解析例である.
!参考 : 分類結果の信頼水準
i.maxlikで作成した分類結果の信頼水準を表示してみよう.
d.rast.leg LC_map_uns_re
{{attach_view(class6.png)}}
黒い箇所ほど信頼水準が低い地点である.
r.report LC_map_uns_re un=h
+-----------------------------------------------------------------------------+
| RASTER MAP CATEGORY REPORT |
|LOCATION: UTM54 Tue Feb 7 22:47:32 2012|
|-----------------------------------------------------------------------------|
| north: 4103757.75 east: 546957.75 |
|REGION south: 3876042.75 west: 290714.25 |
| res: 28.5 res: 28.5 |
|-----------------------------------------------------------------------------|
|MASK:none |
|-----------------------------------------------------------------------------|
|MAP: Rejection Probability for LC_map_uns (LC_map_uns_re in PERMANENT) |
|-----------------------------------------------------------------------------|
| Category Information | |
| #|description | hectares|
|-----------------------------------------------------------------------------|
| 1|0.1% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 17,110|
| 2|0.5% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 18,002|
| 3|1% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 15,519|
| 4|2% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 25,472|
| 5|5% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 67,386|
| 6|10%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 112,386|
| 7|20%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 253,592|
| 8|30%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 289,088|
| 9|50%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 664,842|
|10|70%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 756,468|
|11|80%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 387,851|
|12|90%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 397,786|
|13|95%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 201,025|
|14|98%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 113,870|
|15|99%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 29,634|
|16|100% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 32,008|
| *|no data. . . . . . . . . . . . . . . . . . . . . . . . . . . . | 2,453,011|
|-----------------------------------------------------------------------------|
|TOTAL | 5,835,049|
+-----------------------------------------------------------------------------+
を実行すると, 値が10以上の地点が信頼水準70%以上であることがわかる.
そこで信頼水準70%以上の箇所だけを表示してみよう.
r.mapcalc "LC_map_uns_mask = if (LC_map_uns_re>= 10 , LC_map_uns, null())"
d.rast.leg LC_map_uns_mask
{{attach_view(class7.png)}}
表示してみるとだいぶ白い箇所 (信頼水準が70%以下の地点) が存在していることがわかる.
!1. はじめに
[[前回|植生分類図]], 環境省自然環境保全基礎調査によって作成された植生分類図を扱った. 自然環境保全基礎調査では, 衛星データや空中写真, 現地調査の情報を元に作成されいる. 本章では, 実際にこうした土地被覆分類図を衛星データを用いて作成する方法を紹介する.
!2. 教師つき分類と教師なし分類
土地被覆分類方法の選択肢として, まず教師つき (supervised) と教師なし (unsupervised) の2つの方法がある.
教師つき分類とは, あらかじめ土地被覆分類項目がわかった訓練領域 (=教師データ) をヒントに, 未知の領域の分類を行う方法である.
そのため事前に訓練領域と呼ばれる地域を設定し, 訓練領域内の土地被覆項目 (ex)都市, 森林, 水田など) を把握し, デジタル化しておく必要がある. 一方で教師なし分類では, 画像値の類似性等を機械的に分析し, 分類を行う. そのため教師なし分類によって作成された分類項目に具体的な意味はなく, 地上検証などを用いて, 分類項目と実際の土地被覆項目 (ex)都市, 森林, 水田など)を対応させる必要がある.
つまり,
*教師付き分類は, ラベル付きデータ (事前に分類項目のわかったデータ) を用いて行う分類である.
*教師なし分類は, ラベル無しデータ (事前に分類項目のわからないデータ) を用いて行う分類である.
!3. 教師なし分類
本章では, まず教師なしの土地被覆分類について学ぶ. 教師なし分類の方法についても, ディシジョンツリー法, 最尤法, k-mean法, ISODATA法など, 様々な方法が存在する.
今回は, GRASSにおいて実装されている最尤法による分類を行う. 最尤法による分類では, (複数の) 衛星データのピクセル値 (画素値) が近いものを同じ項目として分類していく手順を,設定した項目数以下になるまで繰り返していく方法である.
用いる衛星データは, [[高分解能衛星Landsat: 植生指標NDVI, 温度画像]]で使用したデータと同じものとする.
GRASSを起動し,
$ grass
下記を選択する.
LOCATION: UTM54
MAPSET: PERMANENT
GRASSにおいて教師なし分類は, i.group, i.cluster, i.maxlikの3種類のコマンドで実装されている.
まずi.groupによって, 分類に使用する衛星データをグループ化する必要がある.
今回は, LandsatのBand1~4のデータを使用する.
for i in 1 2 3 4
do
r.null B$i setnull=0 #事前に0をnullにする.
done
g.region rast=B1
i.group group=landsat sub=main input=B1,B2,B3,B4
groupは, 衛星データをグループ化した際のグループ名である. subgroupを指定することによって, グループ内にも複数のグループを作成できる. 上記コマンド例では, グループ名 : Landsat, サブグループ名 : main とした.
Landsat ---- main ----- B1
|
-- B2
|
-- B3
|
-- B4
次にi.clusterとグループ化したデータを用いて, ピクセル値の類似性を分析し, 分類器 (分類のルール) を作成する.
i.cluster group=landsat subgroup=main sigfile=sig classes=10 report=class_report.txt
sigfileには, 分類器の情報が出力される.(追記:GRASS7ではsignaturefileとする.) classesは, 分類項目数であり, ここでは10とした. reportには, 分類 (クラスラリング ) の質を評価した情報が出力される.
最後にi.maxlikコマンドとsigfileを用いて, 対象地域全体の分類を行う.
i.maxlik group=landsat subgroup=main sigfile=sig class=LC_map_uns reject=LC_map_uns_re
classは, 分類結果のラスターデータ名である. rejectは, 分類の信頼水準を出力するラスターデータ名である.
※classと打ち込んでも動くが警告文が出てくるので、output=LC_map_unsと書いたほうが良い。(追記:2018/02/12 菊島)
さっそく分類結果を表示してみよう.
d.rast.leg LC_map_uns #d.rast.leg は, d.rastとd.legend を組みあわせた素晴らしいコマンドである. 消すときは, d.erase -f
{{attach_view(class1.png)}}
現時点で分類項目に意味はないが, 対象地域が10項目に分類されていることが確認できる.
続いて, 各分類項目と土地被覆項目とを対応させる. 方法は, 現地調査や空中写真判読など様々な方法があるが, 手軽に試すには衛星データのトゥルーカラー画像が良いだろう. まず対象地域が広域なので, 一部地域にズームする.
g.region n=4002967.5 s=3982476 w=412338 e=434511
d.rgb r=B3 g=B2 b=B1
d.mon x1
d.rast LC_map_uns
{{attach_view(class2.png)}}
{{attach_view(class3.png)}}
{{attach_view(class4.png)}}
2枚の画像を見比べてみると, 分類結果の項目番号が
1 水域
2 水域
3 植生
4 植生
5 植生
6 植生
7 都市
8 植生
9 都市
10 裸地
の様に対応していることが読みとれる. (ざっくりとした判読なので正しいとは限らないが, 話を進める.)
次に新しい項目番号を次の様に定義し, 分類結果の項目番号を対応させよう.
古い分類番号 新しい分類番号
1, 2 ---> 11 (水域)
3, 4, 5, 6, 8 ---> 12 (植生)
7, 9 ---> 13 (都市)
10 ---> 14 (裸地)
上記の対応を分類結果であるラスターデータに対して行うには, r.reclassを使う.
cat << EOF > reclass.txt
1 2 = 11 water
3 4 5 6 8 = 12 vegetation
7 9 = 13 urban area
10 = 14 bare
EOF
r.reclass input=LC_map_uns out=LC_map_uns_reclass rules=reclass.txt title=LC_map_uns_reclass
最後に作成した土地被覆分類図の色合いを設定し, 表示してみよう.
cat << EOF > color.txt
11 0:0:255
12 0:100:0
13 255:0:0
14 200:200:0
EOF
r.colors LC_map_uns_reclass color=rules < color.txt
g.region rast=LC_map_uns_reclass
d.rast.leg LC_map_uns_reclass
{{attach_view(class5.png)}}
以上がGRASSによる教師なし分類方法の解析例である.
!参考 : 分類結果の信頼水準
i.maxlikで作成した分類結果の信頼水準を表示してみよう.
d.rast.leg LC_map_uns_re
{{attach_view(class6.png)}}
黒い箇所ほど信頼水準が低い地点である.
r.report LC_map_uns_re un=h
+-----------------------------------------------------------------------------+
| RASTER MAP CATEGORY REPORT |
|LOCATION: UTM54 Tue Feb 7 22:47:32 2012|
|-----------------------------------------------------------------------------|
| north: 4103757.75 east: 546957.75 |
|REGION south: 3876042.75 west: 290714.25 |
| res: 28.5 res: 28.5 |
|-----------------------------------------------------------------------------|
|MASK:none |
|-----------------------------------------------------------------------------|
|MAP: Rejection Probability for LC_map_uns (LC_map_uns_re in PERMANENT) |
|-----------------------------------------------------------------------------|
| Category Information | |
| #|description | hectares|
|-----------------------------------------------------------------------------|
| 1|0.1% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 17,110|
| 2|0.5% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 18,002|
| 3|1% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 15,519|
| 4|2% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 25,472|
| 5|5% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 67,386|
| 6|10%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 112,386|
| 7|20%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 253,592|
| 8|30%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 289,088|
| 9|50%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 664,842|
|10|70%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 756,468|
|11|80%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 387,851|
|12|90%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 397,786|
|13|95%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 201,025|
|14|98%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 113,870|
|15|99%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 29,634|
|16|100% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 32,008|
| *|no data. . . . . . . . . . . . . . . . . . . . . . . . . . . . | 2,453,011|
|-----------------------------------------------------------------------------|
|TOTAL | 5,835,049|
+-----------------------------------------------------------------------------+
を実行すると, 値が10以上の地点が信頼水準70%以上であることがわかる.
そこで信頼水準70%以上の箇所だけを表示してみよう.
r.mapcalc "LC_map_uns_mask = if (LC_map_uns_re>= 10 , LC_map_uns, null())"
d.rast.leg LC_map_uns_mask
{{attach_view(class7.png)}}
表示してみるとだいぶ白い箇所 (信頼水準が70%以下の地点) が存在していることがわかる.