平成２２年度 宇宙利用促進調整委託費～衛星利用の裾野拡大プログラム～ 委託事業

FOSS4G を活用した衛星データ利用のためのオープン・リソースの開発
GIS・衛星データ統合利用プラットフォームを利用した衛星データ利用のチュートリアルの作成

中解像度衛星画像を使用した画像解析チュートリアル目次
１．はじめに...............................................................................................................................................................1
２．解析の準備...........................................................................................................................................................2
２－１．マスクの作成.........................................................................................................................................3
２－２．画像幾何補正.......................................................................................................................................12
３．都内緑地の抽出................................................................................................................................................23
３－１．対象地域と使用データ概要...............................................................................................................24
３－２．各種植生指標の算出...........................................................................................................................26
３－２－１．NDVI........................................................................................................................................27
３－２－２．GVI...........................................................................................................................................33
３－３．領域の決定...........................................................................................................................................37
４．土地利用の評価................................................................................................................................................47
４－１．対象地域と使用データ概要...............................................................................................................49
４－２．教師付き分類による土地利用区分...................................................................................................50
４－２－１．トレーニングエリアの選定.................................................................................................51
４－２－２．最尤法による分類.................................................................................................................67
４－２－３．最大事後確率法による分類.................................................................................................75
４－３．教師なし分類による土地利用区分...................................................................................................84
４－３－１．クラスタリングの実行.........................................................................................................85
４－３－２．分類結果に対するクラスの定義.........................................................................................91
４－４．分類後の処理.......................................................................................................................................95
４－４－１．分類精度の評価.....................................................................................................................96
４－４－２．分類結果の編集...................................................................................................................100
５．災害被災地の抽出..........................................................................................................................................103
５－１．対象地域と使用データ概要.............................................................................................................104
５－２．２時期の差分画像による抽出.........................................................................................................105
５－３．抽出結果を表示.................................................................................................................................111
６．GPSを利用した現地検証...............................................................................................................................120
６－１．緑地抽出における現地検証.............................................................................................................121
６－２．土地被覆分類におけるトレーニングエリアの直接取得............................................................127
６－３．被災地抽出結果を用いた踏査地点の出力....................................................................................136
７．時系列データによる広域土地利用評価......................................................................................................140
７－１．対象地域と使用データの概要.........................................................................................................141
７－２．決定木による土地被覆分類.............................................................................................................142
８．参考文献...........................................................................................................................................................153

１．はじめに本チュートリアルは、「FOSS4G を活用した衛星データ利用のためのオープン・リソースの構築」
プロジェクトの 2010 年度事業の一部として作成されたもので、これから FOSS4G を利用して衛星デー
タを処理・解析を行おうとしている学生や技術者に対して、FOSS4G における具体的な解析作業の実際を説明することを目的としています。読者は本チュートリアルを元に実際の衛星データ解析を行うことで FOSS4G における衛星データ解析の基礎を習得し、実際の業務や高度な研究への適用を行うことができるようにステップアップすることが望まれます。
本チュートリアルでは、以下のような読者を対象として想定しています。

2009 年度に作成された「中解像度衛星画像の変換および表示チュートリアル」の内容を理解し、さらに実際の解析作業の習得を目指している学生や技術者
現在商用アプリケーションを利用しているが、近い将来にオープンソースアプリケーションの導入を検討している学生や技術者及び研究者
衛星解析の基礎的な内容の指導にオープンソースアプリケーションの利用を検討している教育者

また、本チュートリアルで説明する内容の概要は以下のとおりです。

幾何補正やマスクなど衛星解析の準備作業
バンド間演算を用いた解析および画像分類
現地調査結果などを考慮した解析結果の分析手法

これらの内容は衛星データ解析では最も一般的な知識であり、既に商用ソフトウェアなどで利用経験がある読者にとっては、それまで使用していたソフトウェアとの操作性の比較や実用の可能性を検討することに役立てることができるでしょう。
本チュートリアルは上記のような目的から、「中解像度衛星画像の変換および表示チュートリア
ル」と比較すると幾分高度な内容になっています。したがって読者は以下に挙げる知識や技術を習得していることが望まれます。

「中解像度衛星画像の変換および表示チュートリアル」の内容全般
初歩的な画像解析アルゴリズム

なお、本チュートリアルは２０１１年３月に作成されたものですが、衛星観測やソフトウェアなどの技術は日々進化を続けています。読者は常に最新の情報を入手することを心がけておくことが望まれます。

1２．解析の準備
解析を行う前には解析範囲を決定したり、ゲイン補正、幾何補正、大気補正、地形効果補正など
様々な補正を行う等の準備を行う必要があります。
解析範囲の決定は、GRASS ではリージョンを設定して範囲を限定することができるほか、マスクを設定して行政界などのポリゴンで範囲を限定することができます。
各種補正は解析の目的や手法によって行うべきものを選択する必要がありますが、プロダクトを入手した時点で既に実施されている補正もあり、また実施する解析手法によっては補正の必要が無い場
合もあります。必要以上の補正処理は極力避けるようにします。
本章では、GRASS におけるマスクの作成と幾何補正の手法について説明します。

2

２－１．マスクの作成
解析を行う範囲の設定は、リージョンを定義することで、ある程度限定することができます.しかし、市町村単位で解析を行うときなどは、市町村のポリゴン内だけを解析対象にする必要があります。このような複雑な設定はリージョンでは不可能です。このような場合はマスクを作成して範囲を限定します。本章では基盤地図情報（縮尺レベル 25,000）の行政界ポリゴンからマスクを作成する例を示します。
基盤地図情報の行政界ポリゴンは、あらかじめシェープファイルに変換されているもので、座標系
は UTM54N が定義されたデータを用います。基盤地図情報（縮尺レベル 25,000）の行政界データは県ごとに作成されているため、対象とする市区町村のポリゴンを選択してからマスクを作成します。対象とする市区町村を選択するにはその属性値を調べる必要があります。そこで、v.in.ogr コマンド注 2-1-1.
を用いて行政界シェープファイルを GRASS にインポートした後、目的の市町村の属性を調べます。

・属性テーブルから調べる

レイヤーツリーに行政界ポリゴンデータを追加し、「レイヤー設定」ボタンを押して表示されるポップアップメニューから「属性データを表示」を選択します。

「属性データを表示」を選択

図２－１－１．「属性データを表示」メニューを選択

基盤地図情報のシェープファイルは日本語の属性を含んでいるので、日本語の属性値を表示するようにあらかじめ設定しておく必要があります。日本語の属性表示の方法はＡ－１章を参照
してください。
選択すると図２－１－２のようなダイアログが表示されます。

注2-1-1.メニューの[File]-[Import vector data using OGR]3

図２－１－２．属性テーブル

図２－１－２では、"市町村コード"のフィールド名を"AdmCode"と変更しています。市町村を一意に識別するためには必ずこの市町村コードを利用するようにします。属性テーブルの"id"は図形に対して割り振られる値なので、飛び地や離島がある市町村では同一の市町村に２つ以上のid を持ちます。このテーブルを見ると、柏市の市町村コードは 12217 であることがわかります。

・Map Display から調べる

MapDisplay に市町村データを表示させた後、ツールボタンの「ラスター／ベクトルマップの検索」ボタンを押し、表示されるポップアップメニューから[ラスター／ベクトルマップの検索（表示モード）]を選択します。

4

「ラスター／ベクトルマップの検索」
ボタンを選択「ラスター／ベクトルマップの検索（表示モード）」を選択

図２－１－３．Map Display から属性を調べる

ツールを選択後、MapDisplay 上の要素をクリックすると、その図形に対する情報がレイヤー
マネージャの"コマンド出力"タブに出力されます。

5

図２－１－４．MapDisplay から選択して属性を表示した結果対象とする市町村コードを調べたら、次にベクタデータをラスタライズします。ラスタライズを行
うには、メニューの[ファイル]-[マップタイプの変換]-[ベクトル→ラスター]注2-1-2.を選択します。

[ファイル]-[マップタイプの変換]-[ベクトル→ラスター]
メニューを選択

図２－１－５．v.to.rast コマンドの選択

選択すると図２－１－６のようなダイアログが表示されます。まず、"必須"タブで入力ベクタレイ
ヤーと出力するラスタレイヤー名を指定します。
注2-1-2.GRASS コマンド v.to.rast と同値

6

入力・出力レイヤー名を指定

図２－１－６．v.to.rast 必須タブ

続いて"属性"タブを開き、"属性パラメータの列名"の欄に市町村コードを示す列名（今回の例ではAdmCode としている）を指定します。

ラスタライズ時の値に使用する属性を選択

図２－１－７．v.to.rast 属性タブ

7

さらに、"オプション"タブで"ラスター値のソース"に"attr"が指定されていることを確認します。

“ラスター値のソース”が”attr”となっていることを確認

図２－１－８．v.to.rast オプションタブ

以上の設定が完了したら「実行」ボタンを押します。出力ファイルが図２－１－９のように表示されていれば成功です。

図２－１－９．ラスタライズされた行政界ポリゴンデータ

8

最後にラスタライズされた行政界データから、目的の市町村のマスクを作成します。マスクの作成
は、レイヤーマネージャのメニューから[ラスター]-[マスク]注2-1-3.を選択します。

[ラスター]-[マスク]
メニューを選択

図２－１－１０．r.mask コマンドの選択

選択すると図２－１－１１のようなダイアログが表示されます。
"MASK を使用するためのラスターマップ"欄にはマスクの作成元となるラスタマップを選択します。
今回の例では先ほど作成した行政界ポリゴンをラスタライズしたデータを指定します。また、"MASKを使うカテゴリ値"に、本章の最初に調べた市町村コードを入力します。

注2-1-3.GRASS コマンド r.mask と同値

9MASK作成元のラスタレイヤー名を入力

カテゴリ値を指定

図２－１－１１．r.mask オプションタブ

指定が完了したら「実行」ボタンを押します。"コマンド出力"タブに図２－１－１２のように出力されれば成功です。

図２－１－１２．r.mask 実行結果作成されたデータの確認は、図２－１－９の状態からビューを再レンダリングします。図２－１－

10

１３のように指定した市町村だけが表示され、右下に「マスク」と表示されればマスクが有効になっています。

指定した市町村のみが表示される

赤字で「マスク」と表示される

図２－１－１３．マスク適用結果

マスクの実体は「MASK」という名称を持ったラスタレイヤーで、r.mask コマンドはこの「MASK」レイヤーを作成するコマンドであるといえます。したがって、マスクの解除は g.remove コ
マンドで MASK レイヤーを消去すれば可能です注 2-1-4.。また、マスクを一時的に解除して後でまた使用したい場合などは、g.rename コマンドで MASK レイヤーの名前を変更するとよいでしょう。
マスクを作成したら、後続の解析で出力ラスタレイヤーを作成する際に不要に大きなラスタレイ
ヤーが作成されないように、マスク範囲がおさまる程度の広さにリージョンを設定しておくとよいでしょう。

注2-1-4.本来マスクの解除は r.mask -r で行いますが、WinGRASS では動作しません。

11

２－２．画像幾何補正

画像幾何補正とは、衛星データがもつ位置的な誤差を修正することを指します。衛星データは衛星
軌道のわずかなズレやセンサの姿勢制御、センサそのものに内在する様々な要因や、さらには大気屈折率の変化や地形などによって位置的な誤差が必ず存在します。これらの誤差は、解析の目的や要求精度によっては許容範囲内である場合もあるので必ずしも必要な処理ではありませんが、必要と判断される場合はあらゆる解析を行う前に実施しておく必要があります。例えば、５章で紹介する２時期の画像を比較する場合などは、両者の位置を正確に合わせておく必要があります。
一般的に GIS における幾何補正の手順は、変換元と変換先で参照点を取得し、それらを元に変換式を決定して行います。とくに、変換先が実際の地物や観測データである場合は、参照点のことを
GCP（Ground Control Point）と呼びます。すなわち、幾何補正は GCP を真位置として補正を行う画像
上での GCP の位置が真位置に移動するように幾何変換を行います。
画像幾何補正を行うには、GRASS の幾何補正ツールを用いる方法と、QGIS の GeoReferencer プラグインを用いる方法があります。両者ともに活発にアップデートが行われているので、本チュートリア
ルの内容はすぐに古いものになってしまう恐れがありますが、基本的な手順は上記のとおりです。
既に GRASS にインポートされたデータを変換するには GRASS の Georectifier ツールを利用したほうがよいでしょう。データは GRASS にインポートされたものではない場合は QGIS の Georeferencer プラ
グインを用いるほうがよいでしょう。本チュートリアルでは GRASS の幾何補正ツールを用いる手法について説明します。

・幾何補正の準備

幾何補正を行うには、位置を合わせる対象とするデータをあらかじめインポートしておく必要があります。インポートしたデータはマップディスプレイに表示しておきます。以下の例では 25,000 分の
１基盤地図情報を用いて ALOS AVNIR-2 に対して幾何補正を行います。

図２－２－１．幾何補正前の状態

12

まず、幾何補正を行う一連の画像に対してグループを作成します注 2-2-1.。グループとしてまとめられている一連のレイヤーは同じ幾何補正が適用されるようになります。同一のシーンのレイヤーをグループとしてまとめておくことで、一括して処理することができます。
画像グループの作成は、レイヤーマネージャのメニューから[画像]-[画像とグループの作成]-[グループの作成/修正]注2-2-2.を選択します。

[画像]-[画像とグループの作成]-[グループの作成/修正]を選択

図２－２－２．i.group コマンドの選択

選択すると図２－２－３のようなダイアログが表示されます。"必須"タブでは作成するグループの名称を入力します。

注2-2-1.幾何補正ツールの起動時に作成することもできる。
注2-2-2.GRASS コマンド i.group と同値。

13

グループ名を入力

図２－２－３．i.group 必須タブ

続いてオプションタブを開きます。オプションタブではグループに含めるラスタレイヤーとサブグ
ループが必要な場合はサブグループ名を入力します。ここではサブグループは作成せずにレイヤーの登録だけを行います。

レイヤーを選択
続けて選択することで追加していく

図２－２－４．i.group オプションタブ

設定が完了したら「実行」ボタンを押します。

14

図２－２－５．i.group 実行結果

グループの作成が完了したら、幾何補正ツールを起動します。幾何補正ツールはレイヤーマネー
ジャのメニューから[ファイル]-[幾何補正処理]を選択します。

[ファイル]-[幾何補正処理]を選択

図２－２－６．幾何補正ツールの選択

起動すると、図２－２－７のようなダイアログが表示されます。まず、"幾何補正するマップタイプ"では"ラスター"にチェックを入れ、幾何補正を行うレイヤーがインポートされているロケーションとマップセットを選択します。設定が完了したら「次へ」を押します。

15

マップタイプを選択

ロケーション・マップセットを選択

図２－２－７．幾何補正のセットアップ（１）

続いて図２－２－８のようなダイアログが表示されます。"グループ選択"には先ほど作成したグ
ループを指定します。グループを作成していない場合は「グループを作成／編集」ボタンでグループを作成することができます。"マップ出力の拡張子"とは、変換後のレイヤー名の後に付加される文字列です。全ての設定が完了したら「次へ」ボタンを押します。

グループを選択

変換後のレイヤーサフィックスを指定

図２－２－８．幾何補正のセットアップ（２）最後に参照点の取得時に表示するレイヤーを選択します。任意のレイヤーを選択できますが、グループ内のレイヤーを選択するのが妥当です。選択が完了したら、「終了」ボタンを押します。

16

表示レイヤーを指定

図２－２－９．幾何補正のセットアップ（３）「終了」ボタンを押すと、図２－２－１０のように新しいマップディスプレイと参照点管理ダイア
ログが表示されます。マップディスプレイには先ほど選択したレイヤーが表示されます。

図２－２－１０．新マップディスプレイと参照点管理ダイアログ

新しいマップディスプレイのツールボタンは通常とやや異なっています。GCP の取得を行うときは「GCP のセット」ボタンを押してからマップ上をクリックします。

GCPのセットボタン

図２－２－１１．GCP のセットボタン

17参照先のレイヤーが表示されているマップディスプレイはツールボタンからポインタを選択しておきます。
初期状態では、"ラスターの多項式変換方法の選択"では"第１選択"となっています。これは変換に用いる式の次数を表しています。この場合、下の基準点リストには未観測の基準点（各座標値が 0.0 となっている）が３点ある状態になっています。
必要な点数は変換に用いる計算式によって異なります。GRASS では計算式の次数として１次～３次式までを選択することができます。１次式はアフィン変換による線形変換で、以下の式によって表されます。

{

x '=a 1 xa 2 ya 3
y '=b 1 xb2 yb 3

・・・

式１

２次式は以下の式によって変換されます。

{

2

2

x '=a 1a2 xa 3 ya 4 x a 5 y a 6 xy
22
y '=b 1b2 xb3 yb 4 x b5 y b 6 xy

・・・

式２

2

3

３次式は以下の式によって変換されます。

{

2

2

3

2

x '=a 1a2 xa 3 ya 4 x a 5 xya 6 y a 7 x a8 x ya 9 xy a10 y
2
2
3
2
2
3
y ' =b 1b 2 xb3 yb 4 x b5 xyb6 y b7 x b 8 x yb 9 xy b 10 y

・・・

式３

各次数の計算式における係数は１次式は６、２次式は１２、３次式は２０です。参照点を取得する
ことでこれらの係数を求めるので、それぞれ６、１２、２０個の未知変量を求める必要があります。
参照点１組につき２つの方程式ができるので、最小限必要な点数はそれぞれ３点、６点、１０点となり、それ以上の点数を観測した場合は最小二乗法によって係数を求めます。参照点が間違っている可能性を考慮すると、参照点は必要数よりも若干多めに取得しておいたほうが安全です。参照点を追加するには管理ダイアログのツールボタンから「新規 GCP の追加」ボタンを押します。

18

「新規GCPの追加」ボタンを押す

図２－２－１２．GCP の追加ボタン

基準点を観測するには管理ダイアログから未観測の基準点の行を選択し、両方のマップディスプレ
イ上で、図２－２－１３のように参照元、参照先の同一地点をそれぞれクリックします。

同一と思われる箇所をクリック
図２－２－１３．参照点の観測

観測すると、図２－２－１４のように点が記録されます。

19

図２－２－１４．参照点が記録されたところ

参照点を追加するときは、リストの未観測の行を選択してから観測を開始します。そのまま観測す
ると、以前に観測した点を上書きしてしまいます。
同様の作業を必要な点数分追加していきます。参照点を追加するには、管理ダイアログのツールボ
タンから「新規 GCP の追加」ボタンを押します。
基準点数を必要数よりも多く観測し、最小二乗法によって係数を求めた場合、各参照点は変換後残差の分だけずれた位置に変換されるようになります。残差が大きい場合は参照点が間違っている可能性があるので、再度観測を行う必要があります。残差を計算するには、管理ダイアログのツールボタンから「RMS 再解析エラー」ボタンを押します。

「RMS再解析エラー」
ボタン
図２－２－１５．RMS 再解析エラーボタン

ボタンを押すと、図のように各点の右側に残差の計算結果が表示されます。"上位のエラー"とは順変換時の残差、"下位のエラー"とは逆変換時の残差を表しています。

20

図２－２－１６．残差計算結果

残差を確認し、変換を実行するには「幾何補正処理」注 2-2-3.ボタンを押します。

「幾何補正処理」
ボタンを押す

図２－２－１７．幾何補正の実行

幾何補正処理の実行はやや時間がかかります。処理が完了すると図２－２－１８のようなダイアロ
グが表示されます。

図２－２－１８．幾何補正完了通知ダイアログ

続いて図２－２－１９のような警告のダイアログが表示されます。これはベクタデータの属性テー
ブルに関するもので、ここでは関係ありません。

注2-2-3.GRASS コマンド i.rectify と同値

21

図２－２－１９．ベクタデータに関する警告

幾何補正結果を確認し、位置が合っていることを確認してみましょう。ズレが依然として大きい場
合は、参照点を追加したり、既に観測した点を削除するなどして編集する必要があります。

図２－２－２０．幾何補正後

22

３．都内緑地の抽出
都市部における環境評価や景観分析、防災計画の策定などの目的で、都市緑地の分布状況の調査が実施されることがよくありますが、緑地調査に航空写真や衛星画像を用いることによって、広域の緑地分布の概況を効率的に把握することができます。また、過去の航空写真や衛星画像と比較すること
で、緑地の変遷を調べることも可能です。このようなことから、緑地抽出は航空写真や衛星画像が用いられることが多い分野であるといえます。
都市部の緑地を抽出するには、画像分類による方法と、植生指数を用いて抽出する方法が考えられます。本章では、植生指数を用いて衛星画像から緑地を抽出する手法について説明します。

23

３－２－１．NDVINDVI（正規化植生指数）は、クロロフィルが赤色の波長帯（0.64～0.67）を強く吸収し、葉の細胞構造によって近赤外域の波長帯を強く反射することを利用した植生指数です。NDVI は以下の式で算出されます。

NDVI =NIR−REDNIRRED・・・・

式３－２－１－１

NDVI の値は-1～1 の値をとります。正規化することによって観測スペクトルの絶対強度の差異を無視することができるので、異なるシーンとの比較を行うことができるようになります。また、NDVI は無単位の比の値なので、地形効果や大気散乱効果などが各バンドで全て同じであれば、理論的にはそれらの値を除去することができます。
LANDSAT ETM+センサにおける NDVI は以下の式によって算出されます。

NDVI ETM ＋ =band4−band3
band4+band3

・・・・

式３－２－１－２

GRASS を用いて NDVI の算出を行うには Map Calculator注3-2-1-1.というツールを使用します。ここではあらかじめ下準備として、基盤地図情報（縮尺レベル 25,000）の行政界ポリゴンからマスクを作成し、その後で NDVI を算出する例を示します。なお、座標系は UTM54N に設定しています。
MapCalculator を起動するには、レイヤーマネージャのメニューから[ラスター]-[マップカリキュレー
タ]を選択します。

注3-2-1-1.Map Calculator は GRASS コマンドの r.mapcalc の GUI フロントエンドにあたる27[ラスター]-[マップカリキュレータ]メニューを選択

図３－２－１－１．マップカリキュレータの選択

マップカリキュレータを起動すると以下のようなダイアログが表示されます。

図３－２－１－２．マップカリキュレータ

マップカリキュレータでは簡単な四則演算から、論理演算、ビット演算、さらには関数演算や場合分けなど、非常に高度な演算を行うことができます。マップカリキュレータを用いて NDVI を算出す
るには以下のようにします。

① "Name of new map to create"の欄に出力ファイル名を入力。

28

② "Insert mapcalc function"のドロップダウンリストから「float」を選択。float()関数はカッコ内のデータもしくは計算結果を浮動小数に変換する関数です。本章で使用している ETM+画像のように 8bit の画像ではピクセルは符号なし整数として格納されていますが、整数値のまま計算を行うと演算結果も整数値になります。NDVI は上記の式のとおり除算が含まれており、-1～1 の値をとるので、分子と分母をそれぞれ浮動小数に変換しておく必要があります。

“Insert mapcalc function”から”float(x)”を選択
図３－２－１－３．関数の選択③ 画面下の入力欄に"float(x)"と表示されます。括弧内の"x"の文字を削除してカーソルをカッコ内に置き、"Insert existing map"のドロップダウンリストからバンド４のラスタデータを選択する（レイヤー名や＠以降のマップセット名に"+"等の文字があると演算子と重複してしまうので、名称を変更するか、マップセット名を削除しておきます。もしくは、ダブルクォーテーション
でレイヤー名を囲みます）。

29

①カッコ内にカーソルを置く②バンド４のレイヤーを選択
図３－２－１－４．レイヤーの選択④ ③で追加したレイヤー名の後に"-"演算子を追加し、さらにその後にバンド３のレイヤーを追加します。

図３－２－１－５．NDVI 式の分子部分を入力した状態⑤ 同様の要領で演算子とレイヤーを追加し、NDVI の式を完成させる。

30

図３－２－１－６．NDVI 式を完成させた状態⑥ 「Run」ボタンを押し、演算を実行する。

r.mapcalc はリージョンとマスクの設定を反映します。計算が完了したら図３－２－１－７のようなダイアログが表示されます。

図３－２－１－７．r.mapcalc 計算完了ダイアログ

作成された NDVI 画像を確認してみましょう。r.colors コマンドで出力ラスタレイヤーに"ndvi"カラーテーブルを割り当て、カラーバーとともに表示した結果が図３－２－１－８です。

31

図３－２－１－８．NDVI 計算結果

32

３－２－２．GVI

GVI は、タッセルドキャップ変換という重回帰分析で経験的に得られたモデル式によって算出される値のうち、グリーンネスに相当する値です。
タッセルドキャップ変換は LANSAT ETM+の場合はバンド１，２，３，４，５，７の値に対して決まった係数をかけて算出するもので、グリーンネスのほかにブライトネス、ウェットネス、ヘイズの値が算出されます。主成分分析に似ていますが、主成分分析における各成分の意味は明確ではないのに対して、タッセルドキャップ変換によって算出される値は明確に意味づけがされたものです。また、
主成分分析とは異なり、各成分は互いに直交するとは限りません。
GRASS には LANDSAT データに対するタッセルドキャップ変換を行うコマンド i.tasscap が用意されています。GRASS コマンド i.tasscap における、LANDSAT ETM+に対するタッセルドキャップ変換式は以下のとおりです。

Tass=C 1∗band1C 2∗band2C 3∗band3C 4∗band4C 5∗band5C 7∗band7

指数\ 係数C1

C2

C3Greenness-0.3344-0.3544-0.4556Brightness

0.35610.3972Wetness

0.2626Haze

0.0805C4C5C70.6966-0.0242-0.2630

0.39040.69660.2286

0.15960.2141

0.0926

0.0656-0.7629-0.5388-0.0498

0.1950-0.1327

0.5752-0.7775GRASS コマンド i.tasscap を実行するには、レイヤーマネージャのメニューから[画像]-[タッセルド
キャップ植生指数]を選択します。

33

[画像]-[タッセルドキャップ植生指数]を選択

図３－２－２－１．i.tasscap コマンドの選択

選択すると図３－２－２－２のようなダイアログが表示されます。"必須"タブには図のように、各バンドに対応したレイヤーを指定し、出力レイヤー名のプレフィックスを入力します。出力レイヤー
名は prefix.1、prefix.2...という連番が付与されます。

各バンドに対応した
レイヤーを指定

出力レイヤー名のプレフィックスを入力

図３－２－２－２．i.tasscap 必須タブ

34

"オプション"タブでは入力したデータに該当するセンサにチェックを入れます。ETM+センサの場合は"LANDSAT-7 の変換ルールを使用"にチェックを入れます。

“LANDSAT-7の変換ルールを使用”にチェックを入れる

図３－２－２－３．i.tasscap オプションタブ

以上の入力が完了したら「実行」ボタンを押します。"コマンド出力"タブに図のように変換結果が表示されれば変換成功です。

図３－２－２－４．i.tasscap 実行結果

35

変換結果を確認してみましょう。GVI レイヤーは".2"とサフィックスが付与されたレイヤーです。図３－２－２－５は GVI レイヤーに凡例をつけて表示したものです。

図３－２－２－５．GVI 計算結果

36

３－３．領域の決定

植生指数算出の目的は、これら指数と相関が高い植物の資源量や植物の活性を間接的に推定することと、土地被覆としての植生の領域を決定することが挙げられます。本章では、算出された植生指数
から植生域を決定する手法を説明します。
領域の決定の手順は、①植生指数画像を閾値により２値化、②（必要であれば）ノイズの除去、③ラスタ→ベクトル変換という流れで行います。
２値化とはグレースケール画像に対して、ある閾値以上を１、閾値未満を０に変換することによって画像を２種類の領域に分割することです。理想的には、図３－３－１のようなヒストグラムがあった場合、図の黄色の矢印のようなヒストグラムの谷の部分を閾値とするのがよいとされていますが、ヒストグラムによっては谷が明確に現れない場合もあります。

頻度

閾値に適した画素値

画素値

図３－３－１．ヒストグラムの例
ここでは、ヒストグラムを目視して閾値を決定注 3-3-1.したとして、その値による２値化から説明します。
ヒストグラムを表示するには、レイヤーマネージャのメニューから[画像]-[ヒストグラム]を選択す
るか、レイヤーツリーに対象とするラスタレイヤーが追加されている状態で「レイヤー設定」ボタン
を押し、ポップアップメニューから[ヒストグラム]を選択します。

注3-3-1.閾値の決定方法は p-タイル法、モード法、判別分析法などの手法があります。

37

[画像]-[ヒストグラム]を選択

図３－３－２．メニューから d.histogram コマンドの選択「レイヤー設定」ボタンから[ヒストグラム]を選択

図３－３－３．「レイヤー設定」ボタンから d.histgram コマンドの選択

レイヤーマネージャのメニュー選択すると図３－３－４のようなウィンドウが表示されます。ここ
で、ツールバーの一番左の「画像／ラスターファイルのヒストグラム作成」ボタンを押します注 3-3-2.。

注3-3-2.GRASS コマンド d.histogram と同値

38

「画像／ラスタファイルからヒストグラムを作成」
ボタンを押す図３－３－４．ヒストグラム表示ウィンドウ

ボタンを押すと図３－３－５のようなダイアログが表示されます。"必須"タブでは、ヒストグラムを作成するラスタレイヤーを指定します。

ヒストグラムを作成するラスタレイヤーを指定

図３－３－５．d.histogram 必須タブ

39

"オプション"タブではヒストグラム表示に関する様々な設定を行うことができます注 3-3-3.。
設定が完了したら、「OK」ボタンを押すと、ヒストグラムがウィンドウに表示されます。

図３－３－６．ヒストグラム表示

表示されたヒストグラムを参考に閾値を決定したら、２値化を行います。２値化はマップカリキュ
レータで行います（マップカリキュレータの起動については３－２－１章を参照）。マップカリキュ
レータの式を以下のように入力します。

if ( "NDVI レイヤー名" > 閾値, 1, 0 )上式は NDVI レイヤーのピクセル値が閾値以上ならば１、未満であれば０を出力するという意味になります。上式のように条件が真なら 1、偽なら 0 という場合は、後半の「1, 0」は省略することができます。

注3-3-3.本チュートリアル執筆時点ではオプションタブの設定内容は反映されません。

40

図３－３－７．マップカリキュレータ画面

実行結果の例を図３－３－８に示します。赤い部分が２値化により植生と判別された部分です。

図３－３－８．２値化実行結果

図３－３－８の結果である程度領域が決定していますが、この状態では孤立した小領域が多く見られる場合があります。このような小領域が多いことが好ましくない場合はノイズとして除去する処理を行う必要があります注 3-3-4.。小領域を除去するには、まず２値画像のラスタレイヤーに対して r.null コ
マンドで値が"0"のピクセルを"データ無し"に設定します注 3-3-5.。設定後、レイヤーマネージャのメ
注3-3-4.もっとも、本例で使用した ETM+のようにピクセルの地上サイズがある程度大きい場合は小領域の除去は必ずしも妥当な結果をもたらすとはいえません。ノイズ除去処理はピクセルの地上サイズや、ある程度大きな領域だけを抽出したい場合に行うべきです。
注3-3-5.マップカリキュレータの if 式で"データ無し"を出力することもできます。

41

ニューから[ラスター]-[カテゴリー値・ラベルの変更]-[サイズで再分類]を選択します注 3-3-6.。

[ラスター]-[カテゴリー値・ラベルの変更]-[サイズで再分類]を選択

図３－３－９．r.reclass.area コマンドの選択

選択すると図３－３－１０のようなダイアログが表示されます。"必須"タブには入力ラスタレイ
ヤーおよび出力ラスタレイヤーを指定します。

注3-3-6.GRASS コマンド r.reclass.area と同値

42

入力・出力ラスタレイヤーをそれぞれ指定

図３－３－１０．r.reclass.area 必須タブ

続いて"オプション"タブを開きます。今回の場合のように小領域を取り除きたい場合は、"大きい値のオプション>=エリアサイズ"の欄に除去対象となる領域の面積をヘクタール単位で入力します。本章で使用している ETM+画像の場合は地上解像度が 30m なので、1 ピクセルの面積は 0.09 ヘクタールです。逆に 1 ヘクタールはおよそ 3x3 ピクセルよりやや大きい程度です。

除去する対象面積をヘクタール単位で入力

図３－３－１１．r.reclass.area オプションタブ

処理の前後で画像を比較して見ましょう。

43

図３－３－１２．r.reclass.area 適用前後（左が適用前、右が適用後）最後に領域をポリゴン化してみましょう。ラスタデータをベクタデータに変換するには、レイヤー
マネージャのメニューから[ファイル]-[マップタイプの変換]-[ラスター→ベクトル]を選択します注 3-3-7.。

[ファイル]-[マップタイプの変換]-[ラスター→ベクトル]を選択

図３－３－１３．r.to.vect コマンドの選択

選択すると図３－３－１４のようなダイアログが表示されます。"必須"タブでは入力ラスタレイ
ヤーと出力ベクタレイヤー名を指定し、"地物フィーチャ"にはベクトル化する要素の種別を指定します。
この例の場合は"area"を指定します。

注3-3-7.GRASS コマンド r.to.vect と同値。

44

入力・出力レイヤー名を指定

出力ベクタの要素種別を指定

図３－３－１４．r.to.vect 必須タブ

オプションタブではベクトル化に対する様々なオプションがありますが、領域ポリゴンを作成する
場合に重要なオプションとしては"エリアフィーチャの平坦コーナー"があります。このオプションにチェックを入れない場合、作成されるポリゴンはピクセルの四角形の外周をトレースしたものになりますが、チェックを入れることによって階段状に並んだピクセルは斜めにトレースするようになります。

図３－３－１５． エリアフィーチャの平坦コーナーオプションの挙動（左がなし、右があり）

設定が完了したら「実行」ボタンを押します。

45

図３－３－１６．ラスタ→ベクタ変換実行結果抽出した結果が妥当なものであるかどうかは、植生図などとの比較でみることができます。抽出結果の評価については４－４－１章で詳しく説明します。

46

４．土地利用の評価
衛星画像データの特徴の一つに、広域にわたる画像を一度に取得できるという広域性が挙げられます。すなわち、衛星観測を用いることによって、ある時点での地表の現況を広域にわたって瞬時に取得することができるといえます。この性質は土地利用の調査を行う際に最も威力を発揮します。
土地利用の調査とは、現在ある地点の土地がどのように利用されているかを調べることで、都市計
画や国土計画を立案する際の基本的な情報を与えるものです。また、地球温暖化対策や生物多様性保全対策を行ううえでは、森林や農耕地の現況把握は特に重要になります。
このような調査においては、まず衛星画像を用いて土地被覆分類を行い、分類結果と現地調査やその他の資料とを統合的に扱うことで分析を進めていく手法が多く見られます。このように、近年では衛星画像を用いた土地被覆分類は基本的な調査手法として定着しつつあります。
土地被覆分類とは、地上の対象物の分光反射特性や空間的配置をもとに、何らかの分類基準にのっとって画像上のピクセルや領域を複数のクラスに分類していく作業のことです。土地被覆分類の手法
は、大きく分けて教師付き分類と教師なし分類に分けられます。また、ピクセルの値を基に分類を行う場合と、近傍画素との関連をテクスチャ特徴量として定義し、分類に利用する場合もあります。いずれにしても、まずオペレータが土地被覆をどのレベルまで分類するかを決定する必要があります。
分類クラスの決定は調査の目的によって当然異なり、また地域によって出現する可能性がほぼありえないクラス（例：熱帯地域における永久凍土など）もあります。参考として USGS の National Land
Cover Dataset 1992 （http://landcover.usgs.gov/natllandcover.php）において用いられた分類クラスを表４－１に示します。

表４－１．USGS 階層的分類クラス
大項目
水域
開発地

荒廃地

森林

細目
11 水域

12 永久凍土、雪氷

21 低密度居住地

22 高密度居住地

23 商業地域、工業地域、交通機関
31 露岩、砂地、地面

32 採石所、露天掘鉱山、砂利採取場

33 移行地域
41 落葉樹林

42 常緑樹林

43 混交林潅木地

51 潅木地樹林性農作地61 果樹園、ぶどう園、その他

野生／半野生草本

71 草地

農耕地

81 牧草地82 列状に植栽された作物

83 小麦

84 休閑地85 娯楽施設の草地
湿地91 樹林性湿地

92 抽水性草本湿地

実際には衛星データの空間分解能とスペクトル分解能の制約から、上記の表のように細かく土地被覆を分類することはできません。分類クラスの決定は、衛星データの特性や対象地域の特性などを考慮して適切に設定するべきです。

47

分類クラスが決定したら、実際に土地被覆分類を実行します。本章の以降の節では、具体的な分類の手順および分類精度の評価手法について説明します。

48

４－１．対象地域と使用データ概要

解析に使用するデータは Terra ASTER、3A01 プロダクトの画像データです。撮影時の諸元を表４－
１－１、撮影範囲を図４－１－１に示します。

表４－１－１．使用データ諸元
グラニュール ID

諸元

撮影日
AST3A10801060139550912150497 シーン中心

2008 年 1 月 6
北緯 34.691872 、東経 138.200608

図４－１－１．使用した ASTER シーンの撮影範囲対象とする地域は、静岡県牧之原市とし、分類に用いるチャンネルは、VNIR ３チャンネルと SWIR
６チャンネルを選択しました。

49

４－２．教師付き分類による土地利用区分

教師付き分類とは、オペレータが分類の基準となるデータを用意し、それに基づいて対象範囲のピ
クセルや領域を分類していく手法です。最尤法による分類が最も多く用いられていますが、多次元レ
ベルスライス法やデシジョンツリー法などのノンパラメトリック分類法や、広義にはニューラルネッ
トワークやサポートベクターマシンなどの学習システムによる分類処理も含まれます。本チュートリ
アルでは、最尤法および最大事後確率法による分類について説明します。
オペレータが用意する分類の基準となるデータのことをトレーニングデータと呼びます。衛星画像
分類においては、通常トレーニングデータは同一シーン内の領域を指定してそのピクセル値をトレー
ニングデータとして用いることが多いので、特にこの領域のことをトレーニングエリアと呼びます。
教師付き分類では、現地の状況に精通したオペレータがトレーニングエリアを選定することによって、分類精度が高く安定した分類結果を得ることができるという利点がありますが、逆にオペレータ
によって分類結果が異なってしまうということもいえます。すなわち、教師付き分類では良好なト
レーニングエリアを選定することができるかどうかが、分類の良否を左右するといえるでしょう。
以降の節では、実際に教師付き分類を行う手順を説明します。

50

４－２－１．トレーニングエリアの選定

教師付き分類ではまずトレーニングエリアを選定する作業からはじめます。GRASS ではトレーニン
グエリアを選定するコマンド i.class がありますが、このコマンドは Windows 版 GRASS では使用できません。したがって、代替手法として v.digit コマンドを用いてトレーニングエリアを選定する方法注 4-2-1-1.
と、QGIS でデジタイズ後 GRASS にインポートする方法があります。いずれの場合でも、最終的にはトレーニングエリアをラスタライズして使用します。ベクタレイヤーのラスタライズについては２－
１章を参照してください。

・v.digit コマンドによるトレーニングエリアの選定

v.digit コマンドは、レイヤーマネージャのメニューから[ベクトル]-[ベクトルマップの作成]-[デジタイズ]を選択します。

[ベクトル]-[ベクトルマップの作成]-[デジタイズ]
を選択

図４－２－１－１．v.digit コマンドの選択

選択すると図４－２－１－２のようなダイアログが表示されます。"必須"タブでは描画対象のベク
ターレイヤーを指定します。既存のベクターレイヤーを指定するか、新規作成するベクタレイヤー名を入力します。以下では新規作成する場合について説明します。

注4-2-1-1.描画をラスタデータとして保存する r.digit というコマンドもありますが、これも Windows 版 GRASS では使用できません。

51

描画するレイヤーを指定

図４－２－１－２．v.digit 必須タブ

続いてオプションタブでは、デジタイズ画面に背景として表示するためのコマンドラインを指定します。新規作成をする場合は、"存在しない場合は、新規にファイルを作成"にチェックを入れます。

新規作成の場合はチェック

背景表示コマンドを指定

図４－２－１－３．v.digit オプションタブ

背景表示コマンドがわからない場合はレイヤーマネージャのプロパティ画面を開いてコマンドライ
ンを調べてみましょう。

52

現在表示されているレイヤーのプロパティを選択

ステータスバーにコマンドラインが表示されている

図４－２－１－４．表示コマンドを調べる

設定が完了したら「実行」ボタンを押します。図４－２－１－５のように新規マップキャンバスとツールボックスが表示されます。

53

図４－２－１－５．v.digit 起動画面

起動したら、まず新規に作成されたベクタレイヤーの属性テーブルを編集します。v.digit ツール
ボックスの「Open settings」ボタンを押します。

「Open settings」
ボタンを押す

図４－２－１－６．「Open settings」ボタンの選択

ボタンを押すと各種設定を行うダイアログが表示されます。ここで"Table"タブを開くと図４－２－
１－７のような画面が表示されます。ここで、"Add new column"ボタンを押します。

54

「Add new column」
ボタンを押す

図４－２－１－７．Settings ダイアログの"Table"タブ

ボタンを押すと図４－２－１－８のように新しいフィールドが作成されるので、ここにフィールド名を入力します。ここではトレーニングエリアクラス番号を記録したいので、"Class"という名前の整数型のフィールドを作成します。フィールドの入力が完了したら「Create table」ボタンを押します。一
度「Create table」ボタンを押すと、以後はこの画面ではフィールドを編集することができなくなります。
再度フィールドの編集を行うには、v.db.addcol コマンド注 4-2-1-2.などで行う必要があります。

新しく作成されたフィールドに名称とフィールドの型を入力完了したら「Create table」ボタンを押す

図４－２－１－８．フィールドを追加したところ
注4-2-1-2.レイヤーマネージャのメニューから[データベース]-[データベース接続]-[列追加]を選択。

55

設定ダイアログでは、頂点スナップの設定や、背景表示コマンドの変更などを行うことができます。
設定が完了したら、右下の「OK」ボタンを押して設定を保存します。
トレーニングエリアのデジタイズは以下の要領で行います。まず、「Digitize new boundary」ボタンを押します。

「Digitize new boundary」
ボタンを押す

図４－２－１－９．「Digitize new boundary」ボタン

選択後、図４－２－１－１０のようにトレーニングエリアを囲みます。

右クリックで終了

左ボタンで頂点を追加

図４－２－１－１０．Boundary のデジタイズ

56

右クリックでデジタイズを終了すると、図４－２－１－１１のように属性の入力を要求するダイアログが表示されます。しかし、領域要素の境界のデジタイズを行う場合は、この時点では属性を入力する必要はありません。

図４－２－１－１１．属性入力ダイアログ

この時点では境界線が閉じていません。境界線を閉じて多角形にするにはツールボックスの「Movevertex」ボタンを押し、両方の端点をつなぐように一方の頂点を移動します。両端点が設定ダイアログのスナップ間隔よりも近づけば接続とみなされ、図形が閉じられます。

57

「Move vertex」ボタンを押す

左クリックで頂点を選択し、移動後右クリックで確定する
緑色で表示されれば、図形が閉じられている

図４－２－１－１２．Move vertex ツールで図形を閉じる

境界線が描画されましたが、まだこの時点では領域要素になっていません。トレーニングエリアは領域要素で記録する必要があります。領域要素を作成するには、作成された境界要素内に重心要素を追加する必要があります。
重心要素を追加するには、ツールボックスの「Digitize new centroid」ボタンを押し、境界要素内をクリックします。

58

「Digitize new centroid」ボタンを押す

左クリックで境界要素内の任意の地点をクリック

図４－２－１－１３．重心要素の追加

重心要素を追加すると、前出の図４－２－１－１１のダイアログが表示されます。領域要素の属性
は重心要素の属性が使用されるので、ここで分類クラス番号を入力します。

59

①クラスの番号を入力②「submit」ボタンを押す③記録が成功したことを通知する
メッセージが表示される

図４－２－１－１４．属性の入力以上でトレーニングエリアが１つ作成されました。同様の要領で各クラスのトレーニングエリアを作成していきます。トレーニングエリアの指定が完了したら、ベクタ→ラスタ変換を行います。これについては後述します。
デジタイズを終了するときは、ツールボックスの「Save and exit」ボタンを押します。

「Save and exit」
ボタンを押す

図４－２－１－１５．「Save and exit」ボタン

ボタンを押すと、元の v.digit ダイアログの"コマンド出力"タブにデジタイズによって作成された要素の情報が出力されます。領域要素が作成されていることを確認してください。

60

図４－２－１－１６．v.digit コマンド出力タブ

・QGIS によるトレーニングエリアの選定

ここでは QGIS でトレーニングエリアのシェープファイルを作成する方法を説明します。作成したシェープファイルは GRASS にインポートして利用します。なお、QGIS では GRASS のように３つのレイヤー（ファイル）をそれぞれ RGB に割り当ててカラー表示をする機能はないので、あらかじめGDAL ユーティリティコマンド等でカラー画像（またはカラー合成を指定する VRT ファイル）を作成しておいたほうがいいでしょう。
まず、トレーニングエリアのファイルを作成します。ツールボタンの「新規ベクタレイヤ」ボタンを押します。

「新規ベクタレイヤ」
ボタンを押す
図４－２－１－１７．「新規ベクタレイヤ」ボタン

ボタンを押すと図４－２－１－１８のようなダイアログが表示されます。ここでは新規作成されるシェープファイルの設定を行います。トレーニングエリアのシェープファイルを作成するには、まず"タイプ"で"ポリゴン"を指定します。また、「CRS の指定」ボタンを押して空間参照を指定します。さらに、"新しい属性"の欄に属性名（ここでは「class」と入力）を指定、属性のタイプを"整数値"に設定した後、「属性リストに追加する」ボタンを押します。

61

“タイプ”は”ポリゴン”を指定

属性の名称と型を指定後「属性リストに追加する」ボタンを押して属性を追加する空間参照を選択する

４－２－１－１８．新規シェープファイルの作成

設定が完了したら「OK」ボタンを押します。ファイルの保存ダイアログが表示されるので、ファイル名を入力して「保存」ボタンを押します。保存すると、新規作成されたシェープファイルがレイヤに追加されます。

62

図４－２－１－１９．新規レイヤーが追加されたところ

トレーニングエリアの作成を行うには、レイヤーが選択された状態で、ツールバーから「編集モード変更」ボタンを押します。

「編集モード変更」
ボタンを押す各種編集ツールが選択可能になる

図４－２－１－２０．「編集モード変更」ボタン

新規にトレーニングエリアを作成するには、「ポリゴンを入力」ボタンを押し、図４－２－１－２１のようにポリゴンを入力します。

63

「ポリゴンを入力」
ボタンを押す左ボタンで頂点を追加し、最後の頂点で右クリック

図４－２－１－２１．新規ポリゴンの入力

最後の頂点で右クリックすると、ポリゴンの入力が完了し、ポリゴンの属性を入力するダイアログが表示されます。クラス番号を入力後「OK」ボタンを押すと、ポリゴンが作成されます。

64

図４－２－１－２２．属性入力ダイアログ

図４－２－１－２３．ポリゴン入力が完了したところ

異常の要領でトレーニングエリアを作成していきます。ただし、この段階ではまだファイルに書き65

込まれていません。ファイルに書き込むには再度「編集モード変更」ボタンを押します。図４－２－
１－２４のような確認ダイアログが表示されるので、「保存」ボタンを押してファイルに書き込みます。

図４－２－１－２４．書き込み確認ダイアログ

QGIS でトレーニングエリアを作成する利点としては、参照する画像やベクタデータ、WMS 等を手軽に追加できることが挙げられます。同様の機能は GRASS にもありますが、v.digit コマンドで表示コ
マンドを複数指定する必要があり、やや面倒です。また、領域要素の作成も QGIS のほうが簡単に行うことができます。GRASS の v.digit コマンドはトポロジーを考慮したデジタイズを行うことができます
が、トレーニングエリアの選定ではトポロジーは特に必要としないため、トレーニングエリアの選定を行う場合は QGIS のほうが便利でしょう。

66

４－２－２．最尤法による分類

最尤法による分類は画像の教師付き分類では最もよく用いられる手法です。最尤法は、誤分類した
ときの損失の期待値を最小にするという、ベイズの決定則に基づく方法です。具体的には分類対象 x と各クラス c に対して尤度と呼ばれる値を両者の類似度として計算し、尤度が最大であるクラスに分類するというものです。
画像分類で用いられる最尤法においては、分類クラス c と分類対象 x とのマハラノビス距離dM(x,c)、分類クラス c の分散共分散行列を cV とすると、尤度 L(x,c)は以下の式で表されます。

L  x , c=

1
K /2

2 

{

}

1 2
exp − d M  x , c 
2
∣cV ∣
1/ 2

・・・

式４－２－２－１

ただし、K は特徴空間の次元数、|cV|は cV の行列式の値です。マハラノビス距離 dM(x,c)とは分類ク
ラス c の分散共分散行列で x と c とのユークリッド距離を正規化したもので、クラス c の平均ベクトル
を M=(M1, M2, ・・・MK)、x の特徴ベクトルを X=(X1, X2, ・・・XK)とすると、以下の式によって表されます。

d 2M  x , c = X −M cV −1  X −M t・・・

式４－２－２－２式４－２－２－２はユークリッド距離の式に c の分散共分散行列の逆行列が挟まっている形です。
このことから、c の分散共分散が大きい、すなわち広く分布しているクラスのほうが相対的にユーク
リッド距離よりも近い距離の値となります。
式４－２－２－１および式４－２－２－２からわかるように、尤度とは平均がクラス c の平均ベク
トル、標準偏差が cV の行列式の値の正規分布に当たります。一般には対数尤度関数の形に変換し、さらに変形した式を用いて計算量を減らす工夫がなされることが多いです。
以上のことから、最尤法による画像分類を行ううえで必要な情報は、各クラスの平均ベクトルと分散共分散行列です。トレーニングエリアの選定はこれらを決定する作業であるといえます。
GRASS における最尤法による画像分類は以下の手順で行います。準備として、４－２－１章で作成したトレーニングエリアのデータはクラスの値でラスタライズされているものとします。
まず、分類に使用する画像を画像グループを作成して登録します。使用する画像数は特徴ベクトルの次元数 K に該当します。グループの作成は２－２章で行った方法と同様ですが、後述するコマンド
ではサブグループを作成して登録する必要があります。サブグループの作成は図４－２－２－１のよ
うに i.group のオプションタブでサブグループ名を入力して作成します。本章の例では ASTER VNIR ３
バンドと SWIR ６バンドを含むグループを作成しています。

67

“画像サブグループ名”にサブグループの名前を入力

図４－２－２－１．i.group オプションタブ

次に、トレーニングエリアに含まれるピクセルの値から、各クラスの平均ベクトルおよび分散共分
散行列を計算します。レイヤーマネージャのメニューから[画像]-[画像分類]-[最尤法による分類法監視下の入力]を選択します注 4-2-2-1.。

[画像]-[画像分類]-[最尤法による分類法監視下の入力]を選択

図４－２－２－２．i.gensig コマンドの選択

選択すると図４－２－２－３のようなダイアログが表示されます。図のように、ラスタライズされたトレーニングエリア、画像グループ名、サブグループ名、出力ファイル名を指定します。

注4-2-2-1.GRASS コマンド i.gensig と同値。

68

ラスタライズされたトレーニングエリア
画像グループ名
サブグループ名
出力ファイル名

図４－２－２－３．i.gensig 必須タブ

設定が完了したら、「実行」ボタンを押します。正しく作成されたら図４－２－２－４のように情報が出力されます。

図４－２－２－４．i.gensig コマンド出力タブ作成されたファイルを見てみましょう。ファイルはマップセットフォルダ以下の group/"グループ名69

"/subgroup/"サブグループ名"/sig/"sig ファイル名"です（図４－２－２－５）。

図４－２－２－５．sig ファイル保存先

#
#
3032
40.7573 20.4142 49.3493 18.0983 12.2718 12.7622 10.3918 9.68602 10.8206
10.3145
7.14249 7.16651
25.2841 17.2993 183.777
8.86715 7.38421 33.4546 15.2421
4.01386 3.4664 7.71732 5.10756 3.68659
5.25272 4.45222 14.4496 7.26521 3.13954 5.7689
4.18488 3.56644 9.24185 5.40391 2.58071 3.33273 4.12092
3.74329 3.23305 8.07275 5.14238 2.30608 3.07378 2.52656 3.5675
4.23093 3.65337 5.75487 5.08626 2.69279 3.2856 2.74865 2.6058 4.49304#
1304
47.046 25.1679 70.1787 30.2753 18.3681 19.8106 16.6327 15.51 15.1526
6.05774
4.11736 5.31559-4.88543 -6.63939 51.5667
0.704898 0.5393 5.50411 2.84433
0.798092 0.860619 -0.216245 0.972258 2.09617
1.07626 0.976579 1.01315 1.53339 1.09127 2.74153
0.821212 0.678776 -0.185274 0.697523 0.562025 0.782248 1.79052
1.17529 1.1238 -1.01752 0.849517 0.725414 1.02069 0.517479 1.76813
0.356749 0.427152 -0.276713 0.535084 0.582309 0.734992 0.370759 0.367242
1.89611
図４－２－２－６．sig ファイルの内容（一部抜粋）

sig ファイルの中身は、クラスのラベル、クラスの平均ベクトルおよび分散共分散行列の要素が記録

70

されています。分散共分散行列は対称行列なので、下三角成分のみが保存されています。
sig ファイルの作成が完了したら、分類を実行します。リージョンやマスクを設定する場合はここで設定しておくとよいでしょう。これ以前にマスクを設定するとマスクの外にあるトレーニングエリアの情報を取得できなくなってしまいます。分類を行うにはレイヤーマネージャのメニューから[画像][画像分類]-[最尤法による分類(MLC)]を選択します注 4-2-2-2.。

[画像]-[画像分類]-[最尤法による分類（MLC）]
を選択

図４－２－２－７．i.maxlik コマンドの選択

選択すると図４－２－２－８のようなダイアログが表示されます。"必須"タブでは図のように入力、出力レイヤーを指定します。

注4-2-2-2.GRASS コマンド i.maxlik と同値。

71

画像グループ名
サブグループ名
i.gensigで作成されたsigファイル
出力ファイル名

図４－２－２－８．i.maxlik 必須タブ

"オプション"タブの"閾値結果の拒否を保留するラスターマップ名"を指定すると、分類の信頼度が低いピクセルがこのレイヤーに保存されます。このレイヤーの用途としては、分類結果を元に後続作業を行う際に、マスクとして利用するということが挙げられます。
最尤法による分類では、分類対象のピクセルは必ず提示されたクラスのいずれかであるという条件が前提であり、この前提に従っていない場合、例えば極端な例として画像内に森林があるにもかかわらず分類クラスとして森林を設定していないという場合は良好な分類を行う保障がありません注 4-2-2-3.。

注4-2-2-3.もっとも、そもそも最尤法とは「最もそれらしいクラスに分類する」手法であり、必ずしも最良の分類結果を与える手法ではありません。

72

信頼度の低いピクセルを保存したい場合は保存するラスタレイヤー名を指定

図４－２－２－９．i.maxlik オプションタブ

設定が完了したら「実行」ボタンを押します。完了すると図４－２－２－１０のように"コマンド出力"タブに実行結果が表示されます。このとき、各クラスのピクセル値が正規分布から大きく乖離している場合は警告が表示されることがあります。最尤法による分類では各クラスのピクセル値が正規分布に従っていることが前提となっており、トレーニングエリアによって取得されたピクセルの分布が正規分布に従っていない場合は分類精度は悪化します。最尤法ではそのほかにも様々な前提条件が満たされていることを仮定しています。

73

図４－２－２－１０．i.maxlik コマンド出力タブ

分類結果を確認してみましょう。図４－２－２－１１は牧之原市に対して図のようにクラスを定義して分類した結果です。わかりやすく表示するために r.category コマンドを用いてラベルを割り当てています。ラベルの割り当てなどの作業については４－３－２章で説明します。

図４－２－２－１１．分類結果の例

74

４－２－３．最大事後確率法による分類

４－２－２章で述べた最尤法による分類では、分類クラスの生起確率（事前確率）は全てのクラス
について一定であるという前提条件がありますが、最大事後確率法（Maximum a Posteriori : MAP）とは事前確率であるクラスの出現率を与えて計算を行い、分類クラスを決定する方法です。
最尤法によるクラスの推定は以下の式で表すことができます。

 ML  x =argmax f  x∣・・・

式４－２－３－１

式４－２－３－１の x は観測値でピクセル値に該当し、θ は母集団パラメータであり、分類するクラスに該当します。f(x|θ)はすなわち観測値 x がクラス θ から生じた確率（事後確率）です。これに対して最大事後確率推定は g(θ)を θ の事前分布として以下の式で表されます。

 MAP  x=argmax f  x∣ g ・・・

式４－２－３－２

θ の事前分布、すなわちクラスの出現率をどのようにして求めるかが問題になりますが、GRASS における最大事後確率推定では画像ピラミッドが用いられます。
画像ピラミッドとは入力画像に対して図４－２－３－１のように順次解像度を下げた一連の画像のことです。

8x8pix16x16pix

32x32pix

図４－２－３－１．画像ピラミッド

まず、全ての解像度において全ピクセルの対数尤度を計算し、最も解像度が低い画像のピクセルに対して最尤分類によるクラスを割り当てます。次に、最も解像度が低い画像のクラスと隣接ピクセルのクラスおよび１段階解像度が高い画像の各ピクセルの対数尤度とを元に１段階解像度が高い画像の75ピクセルに対してクラスを割り当てます。このとき、クラスの決定には EM アルゴリズム（期待値最大化法）が用いられます。解像度が低いピクセルのクラスを考慮することによって、そのピクセルに由来するピクセルは結果として同じクラスに分類される確率が高くなります。したがって、最終的には同一クラスのピクセルの塊ができやすくなります。このことは、「あるクラス C のピクセルに隣接するピクセルのクラスは C である確率が高い」という仮定を反映したものです。

図４－２－３－２．分類の概念図

各クラスの対数尤度の計算では、オペレータが定義した分類クラスの「サブクラス」が用いられます。ここでサブクラスとは、あるクラスは複数のサブクラスが混合したものであるとして、トレーニ
ングエリアのピクセル分布はそれらサブクラスの集合であると考えます。サブクラスの統計値は最尤法で用いられるものと同様に平均ベクトルと分散共分散行列ですが、トレーニングエリアからそれらを求める方法として混合ガウス分布が用いられます。混合ガウス分布は、ある分布は K 個のガウス分
布の重ね合わせであると考え、以下の式によって表されます。

Kp  x =∑  m N  x∣m ,  m ・・・

式４－２－３－３m=1

上式の π は混合係数、N は平均ベクトル μ、分散共分散行列 σ のガウス分布です。混合ガウス分布を用いた場合の各クラスの対数尤度はそのクラスのサブクラスに対する対数尤度と混合係数の積和になります。最尤法による推定では各クラスの分布が正規分布であることを仮定しているので、混合ガウ
ス分布を用いることによってより精度の高い尤度を計算することができます。
GRASS において最大事後確率法による分類を行うには、まず各クラスの混合ガウス分布を求めるコ
マンド i.gensigset コマンドを実行し、その結果を用いて i.smap コマンドを実行します。まず、４－２－
２章で説明した最尤法による分類と同様に画像グループとサブグループを作成した後、レイヤーマ
ネージャのメニューから[画像]-[画像分類]-[逐次最大事後確率推定法監視下の入力]を選択します。

76

図４－２－３－３．i.gensigset コマンドの選択

選択すると、図４－２－３－４のようなダイアログが表示されます。４－２－２章の i.gensig コマン
ドと同様に、ラスタライズされたトレーニングエリア、画像グループ名、サブグループ名、出力ファ
イル名を指定します。画像グループ名をドロップダウンリストから選択すると"@マップセット"の文字列が付加されますが、i.gensigset コマンドではこの文字列があるとコマンドが失敗しますので、削除しておきます。

ラスタライズされたトレーニングエリア
画像グループ名
サブグループ名
出力ファイル名

図４－２－３－４．i.gensigset 必須タブ

77

オプションタブの"クラスのサブシグネチャの最大値"とは１つのクラスを構成するサブクラスの最大数のことです。GRASS では最適なサブクラス分割を行うように内部で最適化を行っていますが、オ
ペレータがサブクラスの数を制限したい場合はここに値を指定します。デフォルトでは"10"になっています。ここに"1"を指定すると i.gensig コマンドと同じ内容のファイルが作成されます注 4-2-3-1.。

サブクラス数の上限値を指定

図４－２－３－５．i.gensigset オプションタブ

設定が完了したら「実行」ボタンを押します。図４－２－３－６のようにクラス及びサブクラスの作成状況が表示されます。

注4-2-3-1.ただし、ファイルのフォーマットは異なります。

78

図４－２－３－６．i.gensigset コマンド出力タブ作成されたファイルを見てみましょう。ファイルはマップセットフォルダ以下の group/"グループ名"/subgroup/"サブグループ名"/sigset/"sigset ファイル名"です（図４－２－３－７）。

図４－２－３－７．sigset ファイル保存先

79

title: Labels
nbands: 9class:
classnum: 1
classtitle:classtype: 1
subclass:
pi: 0.148039
means: 39.0693 18.9034 58.8659 17.5489 10.654 11.2758 9.39115 8.93056 9.20215
covar:
3.17259 1.45423 11.9045 3.43426 1.30889 1.37988 0.783878 1.61673 1.2808
1.45423 1.48531 6.87302 1.90279 0.821364 0.812477 0.414916 0.883327 0.74408
11.9045 6.87302 71.5036 16.3989 5.48347 6.31135 3.26297 6.90906 4.55985
3.43426 1.90279 16.3989 6.28974 2.18505 2.27059 1.54415 2.65508 2.06291
1.30889 0.821364 5.48347 2.18505 1.73343 1.0014 0.80918 0.919601 0.936699
1.37988 0.812477 6.31135 2.27059 1.0014 1.96658 0.606109 1.10539 0.919523
0.783878 0.414916 3.26297 1.54415 0.80918 0.606109 1.45525 0.760948 0.493273
1.61673 0.883327 6.90906 2.65508 0.919601 1.10539 0.760948 2.45328 1.15773
1.2808 0.74408 4.55985 2.06291 0.936699 0.919523 0.493273 1.15773 1.91969
endsubclass:subclass:
pi: 0.248113
means: 38.632 18.5128 37.6124 14.3255 11.2995 11.2624 9.26628 8.36933 9.77444
covar:
5.71718 2.94551 19.721 4.01799 1.36877 2.25589 1.72236 1.35246 1.0139
2.94551 2.40135 12.552 2.68758 0.89749 1.53809 1.15786 0.841427 0.658218
19.721 12.552 116.892 18.2564 5.348 11.2239 6.08775 2.84609 1.34304
4.01799 2.68758 18.2564 5.67791 1.80629 2.86688 2.00791 1.42506 1.10539
1.36877 0.89749 5.348 1.80629 1.75664 1.09729 0.89225 0.532995 0.432841
2.25589 1.53809 11.2239 2.86688 1.09729 3.12959 1.09257 0.962517 0.473657
1.72236 1.15786 6.08775 2.00791 0.89225 1.09257 2.03745 0.786286 0.635216
1.35246 0.841427 2.84609 1.42506 0.532995 0.962517 0.786286 1.86916 0.585013
1.0139 0.658218 1.34304 1.10539 0.432841 0.473657 0.635216 0.585013 1.83502
endsubclass:

図４－２－３－８．sigset ファイル内容（一部抜粋）

図４－２－３－８のように、sigset ファイルにはサブクラスの混合係数、平均ベクトルおよび分散共分散行列が保存されています。
sigset ファイルの作成が完了したら、分類を実行します。逐次最大事後確率による分類は比較的計算量を要する処理なので、できる限りリージョンを限定しておいたほうがよいでしょう。レイヤーマ
ネージャのメニューから[画像]-[画像分類]-[逐次最大事後確率推定法（SMAP）]を選択します。

80

[画像]-[画像分類]-[逐次最大事後確率推定法（SMAP）]
を選択

図４－２－３－９．i.smap コマンドの選択

選択すると図４－２－３－１０のようなダイアログが表示されます。"必須"タブでは図のように入力、出力レイヤーを指定します。

画像グループ名
サブグループ名
i.gensigsetで作成されたsigファイル
出力ファイル名

図４－２－３－１０．i.smap 必須タブ

オプションタブの"一度に処理するサブマトリックスのサイズ"とは、一度に処理を行う画像のサイズで、上記にあるピラミッドレイヤーの作成などの処理を行う単位に該当します。i.smap コマンドでは内部でスムージングなどの処理を行っていますが、それらのパラメータは処理単位ごとに異なる場合

81

があります。異なるクラスが多く混在する場合は小さい値にしておくとよいでしょう。デフォルトは"128"です。

処理単位の画像サイズを指定

図４－２－３－１１．i.smap オプションタブ

なお、i.smap コマンドでは NULL セルおよびマスクをサポートしていません。とくにマスクは解析
結果に影響を与える恐れがあるので、解除しておいたほうがよいでしょう。設定が完了したら「実行」ボタンを押します。図４－２－３－１２のように処理状況が表示されます。

82

図４－２－３－１２．i.smap コマンド出力タブ

最尤法による分類結果と比較してみましょう。通常の最尤法と比較して小領域が少ないことがわかります。

図４－２－３－１３．分類法の比較（左が最尤法、右が最大事後確率法）

解析結果をマスクで切り取るには、マスク設定後 r.mapcalc コマンドを通してそのままコピーするとよいでしょう注 4-2-3-2.。

注4-2-3-2.r.mapcalc "解析結果のマスクレイヤー = 解析結果レイヤー"とする。

83

４－３．教師なし分類による土地利用区分

教師なし分類とは、与えられたピクセルをクラスタリング等によって機械的に分類する手法です。
機械的に分類を行うため、分類によって作成された各クラスがどのような土地被覆なのかは自明ではなく、したがって分類結果に対してユーザが意味づけを行う必要があります。
クラスタリングの利点としては、分類においてユーザが入力するデータは対象となる画像だけなのでユーザが分類処理に介入する部分がなく、したがって客観的な分類を行うことができるということが挙げられますが、上記のとおり分類結果が必ずしも意味を持っているわけではないので、結果の解釈が難しくなる場合があります。
以降の節では、実際に教師なし分類を行う手順を説明します。

84

４－３－１．クラスタリングの実行

クラスタリングとは与えられたデータを比較的均質であると考えられるグループに分割していく処理のことです。クラスタリングアルゴリズムは、データ数＝クラス数の状態から開始して順次クラスを統合していき、最終的に求めるクラス数になるまで統合を繰り返していく階層的クラスタリングと、
クラスタリングを反復して最適なクラス分割を求める非階層的クラスタリングに大きく分けることができます。GRASS におけるクラスタリングアルゴリズムは後者に含まれる ISODATA 法と呼ばれるア
ルゴリズムです。
ISODATA 法は非階層的クラスタリングの k-means 法を拡張した手法です。k-means 法では図４－３－１－１の手順によってクラスタリングを実行します。

種子点
クラス２

境界線
サンプルデータ
クラス１①特徴空間内に種子点と呼ばれるクラスの重心を無作為にばら撒く②各サンプルを最も近い種子点のクラスに分類する

クラス２

クラス２

クラス１

クラス１

③分割したクラスの重心を計算しなおす④新しい重心を元に再度分類しなおす③～④を再分類が起きなくなるまで繰り返す

図４－３－１－１．k-means 法概要

85

ISODATA 法は k-means 法に以下の条件を加えて分類結果を安定させるように改善したものです。

1.

サンプル数が少ないクラスタはクラスを作成しない

2.

クラスタ間の距離が閾値未満の場合はクラスを統合する

3.

クラスタ内分散が閾値以上の場合はクラスを分割する

GRASS では、全てのピクセルについてクラスタリングを行うと計算量が膨大になってしまうので、まず画像内から適当にサンプルのピクセルを抽出してクラスタリングを行い、各クラスの平均ベクト
ル及び分散共分散行列を計算したのち、最尤法による分類を行います。
以下では GRASS で以上の処理を実行する手順を説明します。クラスタリングを行うには、４－２－２章で説明した最尤法による分類と同様に画像グループとサブグループを作成した後、レイヤーマ
ネージャのメニューから[画像]-[画像分類]-[非監視下分類用クラスタリング入力]注4-3-1-1.を選択します。

[画像]-[画像分類]-[非監視下分類用クラスタリング入力]を選択

図４－３－１－２．i.cluster コマンドの選択

選択すると図４－３－１－３のようなダイアログが表示されます。必須タブでは画像グループとサブグループ、及び出力ファイルと初期クラス数を指定します。

注4-3-1-1.GRASS コマンド i.cluster と同値。

86

画像グループ名
サブグループ名
出力ファイル名
初期クラス数

図４－３－１－３．i.cluster 必須タブ

"設定"タブには ISODATA 法における様々なパラメータを指定します。"サンプリング区間"とはクラスタリングに使用するピクセルのサンプル間隔です。何も指定しない場合は画像内からおおよそ
10,000 ピクセルを取得するように間隔が決定されます。"最大反復数"とは図４－３－１－１の③～④の繰り返しの最大数で、ここで指定した繰り返し数に達した場合は計算を打ち切ります。"収束率"とは図４－３－１－１の③～④の過程でクラスが変化しなかったピクセルの全体に対する割合を指定します。
デフォルトでは 98%ですが、これはクラスが変化したサンプルが 2%未満であれば収束したものとみなして計算を終了します。"クラスターセパレーション"は上記の ISODATA 法の条件２の閾値で、デフォルトでは 0.0、すなわち統合を行いません。"クラスのピクセル最小値"は ISODATA 法の条件１のサンプル数に該当する値です。
これらの値は相互に関連がありますので、設定には多少注意を要します。極端な例としては、サン
プリング間隔に広い値を指定し、クラスのピクセル最小値を高い値に設定すると全てのクラスで条件を満たさなくなる可能性があります。

87サンプリング間隔（10,10のように指定）
収束計算の反復数
収束率
統合するクラス間距離
クラスの最小サンプル数

図４－３－１－４．i.cluster 設定タブ

"オプション"タブの"初期シグネチャを含むファイル名"とは、図４－３－１－１の①の過程で使用される初期クラスの位置に当たるもので、i.cluster の出力結果もしくは４－２－２章で説明した i.gensig
コマンドの出力結果を指定することができます。"最終レポートを含む出力ファイル名"はクラスタリン
グの計算過程を記録したレポートファイルを指定します。レポートファイルは最終的な分類結果をある程度予測するのに必要な情報を提供するので、指定しておいたほうがよいでしょう。

初期値として使用するシグネチャファイル

レポートファイル

図４－３－１－５．i.cluster オプションタブ

88

全ての設定が完了したら「実行」ボタンを押します。"コマンド出力"タブに計算結果が表示されます。出力結果を見る上で重要なのは反復数で、反復数が図４－３－１－４で指定した最大反復数に達している場合は収束条件を満たしていない可能性があります。詳細な情報はレポートファイルを作成しておくと確認することができます。

図４－３－１－６．i.cluster コマンド出力タブ

出力ファイルは４－２－２章で説明した i.gensig コマンドの出力結果と同じファイル形式です。最終
的な分類を行うには i.maxlik コマンドで図４－３－１－３で指定した出力ファイルを指定して分類を実
行します。

89

i.clusterで作成されたsigファイルを指定する

図４－３－１－７．i.maxlik 必須タブ

図４－３－１－８．教師なし分類実行結果

90

４－３－２．分類結果に対するクラスの定義

教師なし分類では機械的に分類を行った後で分類結果の解釈を分析者が行い、各分類クラスに意味づけを行います。クラスに対する意味づけは特に画像に対して何らかの処理を行うわけではなく、具体的には前出の図４－２－２－１１のように凡例を表示するために利用します。
GRASS ではラスタデータのピクセル値をカテゴリとして定義することができ、各カテゴリにカテゴ
リラベルを定義することができます。カテゴリラベルを定義するには、レイヤーマネージャのメニューから[ラスター]-[レポートと統計]-[カテゴリー情報の管理]を選択します注 4-3-2-1.。

[ラスター]-[レポートと統計]-[カテゴリー情報の管理]を選択

図４－３－２－１．r.category コマンドの選択

選択すると図４－３－２－２のようなダイアログが表示されます。まず"必須"タブではカテゴリを編集する対象となるラスタレイヤーを指定します。

注4-3-2-1.GRASS コマンド r.cateogy と同値。

91

カテゴリを編集するラスタレイヤーを指定

図４－３－２－２．r.category 必須タブ

オプションタブは実際にカテゴリを定義していく設定を行います。カテゴリラベルを定義するにはルールファイルを作成して適用します注 4-3-2-2.。ルールファイルは柔軟な書式をサポートしており、様々な形でラベルを定義することができます。ここでは分類結果に対してカテゴリラベルを設定する
例を示します。

0 : nodata
1 : water
2 : forest(shadow)
3 : forest
4 : tea plantation
5 : developed
6 : field
図４－３－２－３．r.category ルールファイル

ルールファイルを作成したら、r.category のオプションタブの"カテゴリラベルルールを含むファイル"の欄に指定します。

注4-3-2-2.msys 等のコマンドラインが利用できる環境であれば、コンソールからの入力も可能です。

92ルールファイルを指定

図４－３－２－４．r.category オプションタブ

カテゴリラベルを設定したら、凡例として表示してみましょう。d.legend コマンドで凡例を表示してみましょう。図４－３－２－５の表示では、カテゴリ番号を非表示にしています。

93

図４－３－２－５．教師なし分類結果に対するカテゴリラベルの表示“カテゴリ番号の非表示”にチェックを入れる

図４－３－２－６．d.legend アドバンスタブ

94

４－４．分類後の処理

土地被覆の分類処理はあくまで計算によって得られた結果であり、必ずしも分類精度が保障されているものではありません。分類処理がどの程度正しく行われたを評価するには植生図や国土数値情報など既存の主題図と比較する方法もありますが、統計的な検定による評価を行うこともできます。分
類結果が思わしくない場合は、教師付き分類であればトレーニングエリアを再度取得し直すなどの作業を行う必要があります。また、教師なし分類では無意味にクラスが分けられていることもあるので、それらを統合して一つのクラスにまとめることも必要になる場合があります。
本章では分類精度を統計的に検定する手法およびカテゴリの編集方法について説明します。

95

４－４－１．分類精度の評価

教師付き分類ではトレーニングエリアを真値として利用することができます。すなわち、トレーニ
ングエリア内のピクセルが分類結果の画像内でそのとおりに分類されたかどうかを調べることで、分類精度を調べることができます。
このような分類制度の評価を行うコマンドとして r.kappa があります。r.kappa を実行するには、レイ
ヤーマネージャのメニューから[画像]-[レポートと統計]-[カッパー分析]を選択します。

[画像]-[レポートと統計]-[カッパー分析]を選択

図４－４－１－１．r.kappa コマンドの選択※
r.kappa コマンドは GRASS データベースへのパスにスペースが含まれていると動作しません。インストール直後のデフォルトの状態では、GRASS データベースは"C:\Users\ユーザ名\Documents\GIS
Database"となっているため、動作しないので注意してください。

選択すると図４－４－１－２のようなダイアログが表示されます。"分類結果を含むラスターマップ
名"には分類結果のレイヤー、"リファレンスクラスを含むラスターマップ名"にはトレーニングエリア
をラスタライズしたレイヤーをそれぞれ指定します。ただし、それぞれのレイヤーのカテゴリは一致している必要があります。したがって、分類結果のラスタレイヤーに対して４－３－２章で説明した
ラベル付けを行った場合は、トレーニングエリアのラスタレイヤーも同様にラベル付けをしておく必要があります。

96

分類結果レイヤーを指定
トレーニングエリアレイヤーを指定

図４－４－１－２．r.kappa 必須タブ

オプションタブではレポートをファイルに書き出すかどうかを指定することができます。ファイルに書き出す場合は"エラーマトリックスと kappa を含む出力ファイル名"に出力ファイル名を指定します。
また、"ワイドレポート"にはチェックを入れておきましょう注 4-4-1-1.。

“ワイドレポート”にチェック出力するレポートファイルを指定

図４－４－１－３．r.kappa オプションタブ

注4-4-1-1.このオプションはレポートの折り返し行数を多くするものです。指定しない場合、カテゴリが５以上あるとエラーマトリックスの途中で折り返されてしまいます。

97

設定が完了したら「実行」ボタンを押します。図４－４－１－４はレポートの例です。

ACCURACY ASSESSMENT
LOCATION: UTM54NFri Jun 25 01:33:05 2010
MASK: none
MAPS: MAP1 = Labels (training_area@makinohara in makinohara)
MAP2 = Maximum Likelihood Classification (mlc@makinohara in makinohara)
Error Matrix
Panel #1 of 1
cat#
M 1
A
2
P
3
2
4
5
6
Col Sum
Cats
1
2
3
4
5
6

1
2992
30
0
1
0
0
3023

MAP1
2
3
19
19
1264 10
3
1536
0
0
0
92
0
0
1286 1657% Commission
1.319261
3.067485
4.119850
1.233046
10.887720
2.826087

4
0
0
0
1602
0
0
1602

5
2
0
60
12
2881
13
2968

6
0
0
3
7
260
447
717

Row Sum
3032
1304
1602
1622
3233
460
11253

% Ommission Estimated Kappa
1.025471
0.981962
1.710731
0.965367
7.302354
0.951688
0.000000
0.985623
2.931267
0.85211937.656904
0.969816

エ
ラ
ー
マ
ト
リ
ッ
ク
ス

オミッション・コミッション
エラー

Kappa
0.940608

Kappa Variance0.000006

カッパ係数・分散

Obs Correct
10722

Total Obs
11253

正答率

% Observed Correct
95.281258

MAP1 Category Description
1: (no description)
2: (no description)
3: (no description)
4: (no description)
5: (no description)
6: (no description)
MAP2 Category Description
1: (no description)
2: (no description)
3: (no description)
4: (no description)
5: (no description)
6: (no description)

図４－４－１－４．r.kappa レポート例

上から順に説明していきます。エラーマトリックスはトレーニングエリアと分類結果の対応表で、
例えば縦方向に見ると、トレーニングエリアのクラス１では正しく分類されたピクセルが 2992、クラス２に誤分類されたピクセルが 30、クラス４に誤分類されたピクセルが 1 であることがわかります。
横方向に見ると、クラス１に分類されたピクセルのうち、本来はクラス２とクラス３であるピクセル
がそれぞれ１９ピクセル、クラス５であるピクセルが２ピクセルあったことがわかります。

98エラーマトリックスの結果を総合的に評価した値が次のオミッション・コミッションエラーの数値です。コミッションエラーは横方向のエラー率で、全ピクセル中別のクラス由来のピクセルの割合を

19192
×100=1.319261 % となります。オミッションエ
3032

示しています。クラス１の例では、

ラーは縦方向のエラー率で、トレーニングエリアのピクセルが他のクラスに分類されてしまった比率

301
×100=1.025471 % となります。分子を誤分類ではな3032を示しています。クラス１の例では

く正しく分類されたピクセル数で計算した値は、それぞれプロデューサ精度、ユーザ精度と呼ばれ、前者はトレーニングエリアの精度、後者は分類画像の信頼度をそれぞれ表しています。さらに、各カ
テゴリの"Estimated Kappa"の値は横方向に対して計算されるカッパ係数で、

K E=

クラス１の例では

Nxii − xi + x + i
Nxi +−x i + x + i

ただし、

N :総ピクセル数
x i + :クラスの横方向の総和
x + i :クラスの縦方向の総和
x ii :正解ピクセル数

2992∗11253−3023∗3032=0.981962 となります。
3032∗11253−3023∗3032

全体のカッパ係数は
M

K =

M

N ∑ x ii −∑  xi + x + i 
i =1

2

i=1

ただし M はクラス数

M

N −∑  x i + x + i 
i =1で計算されます。カッパ係数は無作為な割り当てをどの程度回避しているかを表す指標で、-1～1 の値をとり、値が高いほど精度がよいことを意味します。カッパ係数及びカッパ分散を用いると２つの分類結果に有意差があるかどうかを調べることができます。
本章の例はトレーニングエリアを真値として検定を行っていますが、真値として既存の主題図（植生図など）を利用することで同様の精度評価を行うことができます。ただし、その場合は真値として利用する主題図と分類項目を一致させておく必要があります。

99

４－４－２．分類結果の編集

教師なし分類では、意図しないクラスが作成されることがあります。たとえば図４－３－２－５の分類結果を見ると、森林部が日が当たっている部分と影の部分で別のクラスとして作成されています。
このような場合は２つのクラスを統合して１つにまとめる作業が必要になります。場合によっては、意図的に多めにクラス分割しておいて後で統合することで意図した分類結果を得るようにすることも考えられます。
クラスの統合を行うには、レイヤーマネージャのメニューから[ラスター]-[カテゴリー値・ラベルの変更]-[対話形式による再分類]を選択します注 4-4-2-1.。

[ラスター]-[カテゴリー値・ラベルの変更]-[対話形式による再分類]を選択

図４－４－２－１．r.reclass コマンドの選択

選択すると図４－４－２－２のようなダイアログが表示されます。分類結果を編集する対象のレイ
ヤーと、編集結果の出力レイヤーを指定後、"再分類用のルール入力"の欄に再分類ルールを入力していきます。

注4-4-2-1.GRASS コマンド r.reclass と同値。

100

入力・出力レイヤーを指定

クラス編集ルールを記述

図４－４－２－２．r.reclass コマンド画面

クラスの編集は非常に柔軟に行うことができます。ここでは主題図に対する編集ルールの例を示します。前出の図４－３－２－５の森林（影部）と森林のカテゴリを統合するには以下のようにルールを記述します。

図４－４－２－３．r.reclass ルール入力例

上記のコマンドはクラス２，３を２に統合し、それ以降をひとつずつずらしています。ここでルールを記述していないクラスはデータ無しとなってしまうので注意が必要です。
作成されたレイヤーは図４－４－２－４のようになります。

101

図４－４－２－４．クラス編集結果

102

５．災害被災地の抽出
衛星リモートセンシングの特徴として、同一地点を周期的に観測することができるという特徴があります。このことは地震、豪雨被害、山林火災などの大規模災害が発生したときにいち早く災害地域を特定することに適しており、被害の拡大を食い止めたり災害救助に役立てることができます。特に中解像度の衛星センサは観測範囲が広いことから、効率よく被災地を抽出することに適しています。
被災地抽出作業では、最終的には被災地と判断された部分を地図上にプロットすることが目的となりますが、災害救助や被害の拡大防止という最終目的を考慮すると、複雑な解析を行って正確に抽出するというよりはなるべく簡単な作業で迅速に被災地の概略地点を特定して報告することが重要です。
画像から被災地抽出の手法としては、実際は災害後の画像から目視による判断を行う場合が多いですが、上記のように衛星観測の周期性を利用して災害前後の画像を比較することで、迅速かつより精度が高い被災地抽出を行うことができる場合があります。
そこで本章では、災害前後２時期の ALOS AVNIR-2 画像から被災地の抽出を行う手法について説明
します。

103

５－１．対象地域と使用データ概要

本チュートリアルでは２００８年６月１４日に発生した岩手・宮城内陸地震を事例として取り上げます。岩手・宮城内陸地震は岩手県奥州市および宮城県栗原市において最大震度６強を観測し、両県を中心に被害が発生しました。この地震による被害の特徴としては、同程度の地震と比較すると建造物の被害は比較的少なく、山体崩壊やがけ崩れ、河道閉塞などの土砂災害が多く発生した事が挙げられます。
解析に使用した衛星データは震災前後の ALOS AVNIR-2、1B2-Georeference プロダクトの画像データです。撮影時の諸元を表５－１－１、また撮影範囲を図５－１－１に示します。

表５－１－１．使用データ諸元
シーン ID

ALAV2A038802790

ALAV2A129822820

諸元
撮影日

2006 年 10 月 17

シーン中心

北緯 39.022 、東経 140.777

ポインティング角

34.30

撮影日

2008 年 07 月 02

シーン中心

北緯 38.782 、東経 140.887

ポインティング角-2.50

図５－１－１．使用した AVNIR-2 シーンの撮影範囲

104５－２．２時期の差分画像による抽出

前記のとおり、解析対象とする地震では土砂災害が多く発生しています。山体崩壊やがけ崩れなどの土砂災害では、発生地点では植生が失われて地表がむき出しになっているものと考えられます。すなわち、２時期の観測データから土砂災害発生地点を特定するには、災害の前には植生があり、災害後に植生が失われた場所を抽出することで絞り込むことができます。そのようにして抽出した領域から、２時期の撮影時期の間に開発などによって植生が失われた箇所を除くことで、災害発生箇所を特定することができます。
植生が失われた箇所を抽出するには、２時期間で画像差分をとる方法と、両者の植生指数画像を作成してそれらの差分をとる方法が考えられます。いずれにしても、２時期の画像間で演算を行う場合
は両方の画像の位置合わせを事前に行っておく必要があります。

・単純差分画像による抽出

２時期の画像における差異を調べるに当たってまず考えられる手法は、単純に両者の差分を取ることです。当然両者の画像では撮影時期、特に季節が異なるため色調は異なりますが、単純な差分画像においては、地上の状況に変化がない箇所では色調の差だけが残った画像ということがいえるので、差分画像上において周辺と明らかに色調が異なる箇所ではやはり地上の状況に何らかの変化が起こっている箇所であるということがいえます。このことから、単純に差分をとる方法は手法の簡単さもあり、被災地抽出では有効な手段の一つです。
まずは、２時期のデータで差分を計算してみましょう。リージョンを設定後、マップカリキューレ
タでバンド２に対して「"期末レイヤー"-"期首レイヤー"」のように入力します。

図５－２－１．差分計算入力例

本章では、リージョンは重複部分が全て入る範囲に設定して計算を行っています。r.mapcalc ではこのような場合でも重複部分だけを計算対象にします。
計算結果を表示してみましょう。図５－２－２はバンド２に対する２時期の差分結果です。

105図５－２－２．バンド２に対する差分計算結果

この画像では総じて青色になっていますが、所々赤く表示されている部分があります。カラーバーに示されているように、赤い部分は高い値が記録されている箇所、すなわち期末の画像では高い画素
値で期首の画像では低い画素値の部分であることを示しています。土砂災害発生箇所は地表が露出するため高い反射を示すことが多いので、図で赤く表示されている箇所は災害発生箇所である可能性があります。
以下に、図５－２－２で赤く表示されている箇所と、期首、期末のトゥルーカラー画像を示します。
災害地域をある程度抽出することができますが、期末画像の雲なども同時に抽出してしまう欠点があ
ります。

106

表５－２－１．主な抽出箇所の画像比較
場所

差分画像

期首画像

荒砥沢ダム付近

駒の湯温泉付近

湯の倉温泉付近

107

期末画像

・NDVI 画像に対する差分画像による抽出

NDVI 画像の差分をとる方法は、単純な差分をとる方法よりも手間がかかりますが、-1～1 の間に正規化された値同士の差分なので差分画像は 2～-2 の値をとり、結果を定量的に取り扱いやすくなります。
また、雲の影響をある程度軽減することができる上、結果の画像は植生の喪失という意味づけが明確にできるので、解析の根拠も明白になります。
以下の例は、期首、期末それぞれで NDVI 画像を作成し、それらの間で差分をとった結果です。

図５－２－３．NDVI 差分画像

これではわかりにくいので、NDVI が 0.5 以上減少した箇所だけを抽出してみましょう。GRASS コ
マンド r.null を起動し、setnull オプション（コマンドダイアログでは"オプション"タブの"NULL 値をセットするセルの一覧"）のところに「-0.5-2」と入力して実行します。

108

“NULL値をセットするセルの一覧”の欄に”-0.5-2”と入力

図５－２－４．r.null オプションタブ入力例

さらに、カラールールで全ての色を単一の色になるように以下のように設定します。

「0%」「100%」と入力し、両方を赤色にする

109図５－２－５．カラールールの設定例

上記の手順で被災地を抽出し、トゥルーカラー画像と重ね合わせた結果が図５－２－６です。中央
やや北の部分にある雲を抽出していないことがわかります。

図５－２－６．栗駒山周辺と荒砥沢ダムの被災地抽出結果

110５－３．抽出結果を表示

被災地抽出作業の最終目的は、地形図等の地図上に抽出箇所を図示して被害状況を整理することです。GRASS 上では前節の図５－２－６のように表示することで被災箇所をわかりやすく表示することができますが、その他の GIS ソフトなどでも利用可能にするために適切な形でデータをエクスポート
しておいたほうがよいでしょう。また、そうすることによって現地調査結果との統合も容易になります。
本節ではラスタデータで出力後 QGIS 上で基盤地図情報 WMS にオーバーレイする方法と、ベクタデータで出力後 Google Earth にオーバーレイする方法を説明します。

・QGIS 上で基盤地図情報 WMS とともに表示

前節の図５－２－６は表示上の設定で２値画像のように表示しているだけで、実際のデータは浮動小数の値が記録されています。このままエクスポートしても QGIS で表示することはできますが、２値化してピクセルもバイト型にしてエクスポートすることで、データ量を圧縮できる上、取り扱いが容易になります。
リージョンを設定後、r.mapcalc で２値画像を作成します。２値化の要領は３－３章と同様ですが、
既に NULL 値が設定されているので、r.mapcalc で以下のように入力することで２値化することができます。

if( not( isnull( "NDVI 差分レイヤー" ) ) )上式は NULL 値に設定されているピクセルは０、それ以外は１を出力します。この時点で再度
NULL 値を設定することもできますが、GRASS ではエクスポート時に NULL 値を出力することができません注 5-3-1.。したがって、必ずしもここで NULL 値を出力する必要はありません。

注5-3-1.r.out.gdal コマンドには"nodata"オプションがありますが、本チュートリアル執筆時点では出力ファイルにnodata の値を設定することができません。

111

図５－３－１．２値化結果作成された２値画像を GeoTIFF 形式でエクスポートします。r.out.gdal コマンドの"必須"タブでエク
スポートするレイヤーと出力ファイル名を指定し、"オプション"タブの"ファイルフォーマット"を"GTiff"、"ファイルタイプ"に"Byte"をそれぞれ指定します。

フォーマットを”GTiff”

ファイルタイプを”Byte”

５－３－２．r.out.gdal オプションタブ作成されたファイルを QGIS 上で表示してみましょう。QGIS に出力ファイルを読み込んだ状態で、112ラスタレイヤプロパティダイアログの"透過性"画面を開き、"データ無しとする値"に"0"を入力します。

② “データ無しとする値”に”0”を指定① “透過性”タブを開く

図５－３－３．QGIS ラスタレイヤプロパティダイアログの透過性画面

出力ファイルを基盤地図情報にオーバーレイして表示してみましょう。出力ファイルは GRASS の
ロケーションと同じ空間参照が定義されているので、基盤地図情報をオーバーレイする際には空間参照に注意してください。

113

図５－３－４．基盤地図情報 WMS と被災地抽出結果のオーバーレイ

・Google Earth 上に表示

Google Earth 上に表示するには、上記のラスタデータをグラウンドオーバーレイとして表示する方
法と、ラスタ→ベクトル変換後に KML 形式に変換して表示する方法があります。本節では後者の手法について説明します。
ラスタ→ベクトル変換は３－３章で説明したとおりですが、ここでもやはり差分レイヤーを２値化する必要があります。また、ラスタ→ベクトル変換を行うためには２値化の際に抽出部分以外を NULL
値に設定する必要があります。そこで、r.mapcalc には以下のような式を入力します。

if( "NDVI 差分レイヤー" < -0.5, 1, null() )上式で作成したラスタレイヤーは抽出箇所が１、それ以外が NULL になっています。この状態で３－３章で説明したノイズ除去を行うこともできますが、ベクトル化した後にノイズ除去を行うこともできます。ベクトルデータのノイズ除去を行うには、レイヤーマネージャのメニューから[ベクトル][ベクトルマップの作成]-[円滑化または単純化]注5-3-2.を選択します。

注5-3-2.GRASS コマンド v.generalize と同値

114

[ベクトル]-[ベクトルマップの作成]-[円滑化または単純化]
を選択

図５－３－５．v.generalize コマンドの選択

選択すると図５－３－６のようなダイアログが表示されます。小領域を削除する場合は、"一般化アルゴリズム"で"remove_small"を選択し、"最大許容差値"に削除し対象領域の面積を指定します注 5-3-3.。

入力・出力レイヤーを指定

アルゴリズムを”remove_small”に設定

最大許容差値に削除する面積を指定

図５－３－６．v.generalize 必須タブ

注5-3-3.v.generalize コマンドは非常に多くのアルゴリズムを用いて図形の単純化、平滑化を行うコマンドですが、使用するアルゴリズムによって要求パラメータが異なります。そのため、日本語訳に違和感があるかもしれません。なお、"最大許容差値"は threshold パラメータと同値です。

115

図５－３－７．v.generalize 実行結果
（赤色の部分はもとの差分ラスタレイヤー）

小領域が多いと描画に時間がかかる恐れがあるので、できる限り削除しておいたほうがいいでしょう。
ベクトル化が完了したら、エクスポートを実行します。レイヤーマネージャのメニューから[ファイル]-[ベクトルマップのエクスポート]-[OGR を使用した複数フォーマットのエクスポート]を選択します注 5-3-4.。

注5-3-4.GRASS コマンド v.out.ogr と同値。

116

[ファイル]-[ベクトルマップのエクスポート]-[OGRを使用した複数フォーマットのエクスポート]を選択

図５－３－８．v.out.ogr コマンドの選択

選択すると図５－３－９のようなダイアログが表示されます。"必須"タブではエクスポートするベ
クタレイヤーと出力ファイル名を指定します。出力ファイル名の指定では「参照」ボタンがありません。従ってフルパスを手入力する必要があります。

入力レイヤー、出力ファイルを指定

図５－３－９．v.out.ogr 必須タブ

続いて入力タブを開きます。ここでは出力する要素形式を指定します。今回の例では抽出した領域のポリゴンだけをエクスポートしたいので、"area"だけにチェックを入れます。

117

“area”だけにチェック

図５－３－１０．v.out.ogr 入力タブ

続いて"作成"タブを開きます。一番上の"OGR レイヤー名"の欄を指定した場合、"必須"タブで入力したファイル名はフォルダ名となり、その中にこの欄で指定したファイルが作成されます。空欄に下場合は"必須"タブで指定したファイル名がそのまま使用されます。
OGR フォーマットは OGR がサポートしているファイルフォーマットを選択することができます。
この中には KML も含まれていますが、この時点で KML を作成しても座標値が UTM の値になっているので無効な KML ファイルとなってしまいます。ポリゴン形式を保存することができるフォーマットであれば何でも構いませんが、ここではとりあえず"ESRI_Shapefile"を選択しておきます。

出力フォーマットを選択

118

図５－３－１１．v.out.ogr 作成タブ

以上の設定が完了したら「実行」ボタンを押します。シェープファイルが作成されたら、作成されたシェープファイルの座標値を経緯度に変換しつつ KML ファイルを作成します。GDAL ユーティリ
ティコマンド（FWTools に同梱）の ogr2ogr を使用して、コマンドプロンプトから以下のように入力します。

C:\Data> ogr2ogr -f "KML" -t_srs EPSG:4326 "出力ファイル名" "入力ファイル名" ⏎入力ファイルと出力ファイルの指定順序に注意してください。作成された KML ファイルは Google
Earth がインストールされているマシンで閲覧することができます。

119

６．GPS を利用した現地検証
衛星画像解析による分析結果は主観的、客観的に関わらずあくまで推定結果であり、必ずしも現地の状況を確実に表している保障はありません。一方、現地調査は調査対象の現況を調べる最も確実な手法ですが、広域を一度に現地調査を行うことは現実的ではありません。そこで、両者を併用することで広域の現況をより正確に把握することができます。
一言に現地調査といっても調査項目、手法は様々です。植生調査であれば現況写真を撮るだけの場
合もあれば、詳細な毎木調査を行う場合もあるでしょう。そこで、これらの調査と平行して GPS による位置情報を取得することで、それらの現地調査と衛星解析結果を統合して分析を行うことができる
ようになります。また、処理を行う前にグラウンドトゥルースとして現地調査結果を入力することによって解析精度を向上させることができます。とくに幾何補正や教師付き分類では大きな威力を発揮します。
GPS 測量には様々な測位手法があり、最も精度の高いスタティック測量や RTK ではミリ単位の精度で測位を行うことができます。もちろん精度が高いことに越したことはありませんが、観測機材が高価であったり、観測そのものや後処理に手間と時間を要したりします。一方、最も精度が低い単独測
位では位置精度は数メートル程度で、条件が悪い場合は数十メートルほど誤差が出る場合があります。
しかし、測位作業及び後処理の手間はほとんどないため、他の現地調査と並行して測位を行ってもそれほど作業の負担になることはありません。また、最近ではハンディ型 GPS が多く普及してきており、安価で入手できるようになってきています。本チュートリアルで利用している中解像度衛星画像の地上解像度がおおよそ 10～30 メートル程度であることを考慮すると、単独測位による GPS でも充分であることも多いでしょう。
そこで、本章ではこれまでの章で述べた例と現地調査を組み合わせるための概念的な手法について説明します。本章で述べる手法はあくまで一例であり、実際の調査業務においては適切ではない可能性もあります。

120

６－１．緑地抽出における現地検証

３－２章で説明した植生指数はあくまで指数であり、植物の現存量や活性度に対する絶対的な指標を与えるものではありません。算出された植生指数の値と現存量や活性度などの実際の値との関連付けを行うには、現地調査による実測値が欠かせません。また、解析の結果、植生域であると判断された箇所には本当に植生があるのかを確認するにはやはり現地調査が必要になるでしょう。
そこでこのような場合は対象地域から無作為に調査地点を選定して現地調査を行う必要があります。
空間的な無作為抽出を行う手法として、GRASS を用いてランダムに点をばら撒く方法があります。
ランダムに点をばら撒くには、レイヤーマネージャのメニューから[ベクトル]-[ポイントの作成]-[ランダムなポイントを作成]を選択します注 6-1-1.。

[ベクトル]-[ポイントの作成]-[ランダムなポイントを作成]
を選択

図６－１－１．v.random コマンドの選択

選択すると図６－１－２のようなダイアログが表示されます。作成するベクトルレイヤー名と発生させるランダム点数を入力します。v.random コマンドはリージョン内に対してランダム点を発生させるので、あらかじめ調査対象範囲にリージョンを設定しておきます。

注6-1-1.GRASS コマンド v.random と同値。

121

作成するレイヤー名
作成するポイント数

図６－１－２．v.random 必須タブ

図６－１－３．v.random 実行結果

122

v.random は図６－１－３のようにリージョン全体に点を作成してしまいます。対象範囲外の点を削除するには、レイヤーマネージャのメニューから[ベクトル]-[別のベクトルマップと共に検索]を選択します注 6-1-2.。

[ベクトル]-[別のベクトルマップと共に検索]を選択

図６－１－４．v.select コマンドの選択

選択すると図６－１－５のようなダイアログが表示されます。"必須"タブでは、図のように取り出す対象のレイヤー、検索条件とするレイヤー、および出力レイヤー名をそれぞれ指定します。上記の例ではベクトルマップ（Ａ）にはポイントデータ、ベクトルマップ（Ｂ）には行政界ポリゴンを指定します。

注6-1-2.GRASS コマンド v.select と同値。

123

抽出対象のレイヤー名
検索条件とするレイヤー名
出力レイヤー名

図６－１－５．v.select 必須タブ

"選択"タブでは両方のベクトルマップにおける対象要素種別とサブレイヤーの番号をそれぞれ指定
します。v.select コマンドには要素を属性で選択するオプションはないので、上記の例の場合はベクト
ルマップ（Ｂ）で指定したレイヤーは、行政界レイヤーのポリゴンデータから対象の市町村ポリゴン
だけを v.extract コマンド等であらかじめ取り出しておいたものを指定する必要があります。

図６－１－６．v.select 選択タブ124

さらに"オプション"タブでは実行する演算を選択することができますが、本チュートリアル執筆時
点では"overlap"だけが選択可能です。"overlap"は「ベクトルマップ（Ｂ）に重複する要素をベクトル
マップ（Ａ）から抽出」する演算です。

図６－１－７．v.select オプションタブ

設定が完了したら「実行」ボタンを押します。図６－１－８のように行政界内の点だけが抽出され
ます。

125

図６－１－８．v.select 実行結果

ランダム点を発生させ、行政界で抽出したら、さらに私有地内や水面に落ちている点を除去して調査対象点を選定します。植生域のみを現地調査対象とする場合はさらに対象点を絞り込む必要があります。最終的な対象点はかなり絞られることを考慮して、ランダム点はやや多めに設定しておいたほうがよいでしょう。
植生の領域が正しく抽出できたかどうかを調べる場合は、現地調査地点を多めにして植生／非植生を調べ、"正解率"を算出することができます。また、胸高断面積合計や立木密度など定量的な数値と
NDVI との関連を調べたい場合は、回帰直線や相関係数を求めることによって関連の度合いを調べることができるでしょう。

126

６－２．土地被覆分類におけるトレーニングエリアの直接取得

現地調査によるトレーニングエリアの取得は最も確実にトレーニングエリアを取得する方法です。
現地でトレーニングエリアを取得するには、GPS 移動計測で領域を囲むように踏査してもよいですが、より簡便に行うには点測定を行ってバッファを発生させ、トレーニングエリアにする方法があります。
まず、現地調査では GPS による点測定を行い、測定点に属性として分類項目（番号）と、その測定点の周囲何メートルまでが同じ分類項目かを記録しておきます。
例として、現地調査結果を図６－２－１のように取りまとめたとします。

XY
242510.77247406.63
245968.84
247185.1
243660.32
239313.43
241000.96
247782.21
248040.75250543.5
242987.74
243437.23
241827.18
243322.08
226845225270.21
244604.43
230153.11230809.36
230624.64
229199.21
226555.75
218302.12215829.58
219252.47

3854088.54
3849297.37
3849569.73
3847746.58
3846135.84
3846307.58
3840439.01
3850749.15
3847035.18
3852358.52
3855832.36
3857024.64
3858415.22
3859618.96
3851833.28
3853436.62
3834787.17
3837632.23
3838373.54
3837940.91
3859971.88
3869831.64
3843803.31
3848257.66
3849915.13

class distance
5
200
5
200
2
100
3
100
2
80
1
200
2
30
3
150
4
300
6
150
2
250
6
100
5
250
1
80
5
150
4
50
5
150
4
100
4
50
6
20
1
250
1
200
3
150
3
150
3
70

図６－２－１．GPS 観測取りまとめ例

上記のテキストファイルを GRASS コマンド v.in.ascii でインポートした後、レイヤーマネージャのメニューから[ベクトル]-[バッファリング]を選択します注 6-2-1.。

注6-2-1.GRASS コマンド v.buffer と同値。

127

[ベクトル]-[バッファリング]を選択

図６－２－２．v.buffer コマンドの選択

選択すると図６－２－３のようなダイアログが表示されます。必須タブにはバッファ発生元のレイ
ヤーと出力レイヤー名をそれぞれ指定します。

バッファ発生元のレイヤー
出力レイヤー名

図６－２－３．v.buffer 必須タブ

続いて"距離"タブを開き、"バッファー距離に使用する列名"欄に現地調査で記録した距離（この場
合は"distance"）を指定します。

128

バッファ距離を定義したフィールドを指定

図６－２－４．v.buffer 距離タブ以上の設定が完了したら「実行」ボタンを押します。実行結果を図６－２－５に示します。各点の周りにそれぞれ異なる距離のバッファが作成されていることがわかります。

図６－２－５．v.buffer 実行結果

129

ただし、v.buffer コマンドはバッファ発生元から属性を引き継ぎません。そこで、バッファ発生元から属性テーブルを明示的にコピーする必要があります。ただし、単純にコピーしただけではテーブルの各行との対応付けがうまくいきません。そこで、以下の手順でバッファの各要素に属性をコピーします。
まず、バッファ発生元の要素種別をポイントから重心要素に変換する必要があります。要素種別の変換を行うにはレイヤーマネージャのメニューから[ベクトル]-[ベクトルマップの作成]-[オブジェクト
タイプの変換]を選択します注 6-2-2.。

[ベクトル]-[ベクトルマップの作成]-[オブジェクトタイプの変換]を選択

図６－２－６．v.type コマンドの選択

選択すると図６－２－７のようなダイアログが表示されます。"必須"タブには入力レイヤーと出力
レイヤー名をそれぞれ指定します。

注6-2-2.GRASS コマンド v.type と同値

130

入力レイヤー、出力レイヤー名を指定

図６－２－７．v.type 必須タブ

続いてオプションタブの"Conversion"のドロップダウンリストから"point to centroid"を選択します。

“point to centroid”を選択

図６－２－８．v.typ オプションタブ

設定が完了したら「実行」ボタンを押します。作成される重心要素のレイヤーには属性情報が持ち越されています。次に、レイヤーマネージャのメニューから[ベクトル]-[ベクトルマップの重ね合わせ]-[ベクトルマップのパッチ]を選択します注 6-2-3.。

注6-2-3.GRASS コマンド v.patch と同値。

131

[ベクトル]-[ベクトルマップの重ね合わせ]-[ベクトルマップのパッチ]を選択

図６－２－９．v.patch コマンドの選択

v.patch コマンドは複数のベクタレイヤーを統合するコマンドです。バッファポリゴンと先ほど選択した重心要素を統合することで、バッファポリゴンが属性を持つようになります。選択すると図６－
２－１０のようなダイアログが表示されます。図のように、"入力するベクトルマップ名"の欄には統合する２つのレイヤー（この場合はバッファレイヤーと重心レイヤー）を指定し、出力レイヤー名を入力します。

合成する複数のレイヤーを指定
出力レイヤー名を指定

図６－２－１０．v.patch 必須タブ

132

実行後出力されるベクタレイヤーの各バッファポリゴンには重心要素の要素番号（cat の値）が付与されますが、重心要素の属性値は持ち越すことができません。最後に、既存の属性テーブルとバッ
ファポリゴンと重心要素の合成レイヤーとを接続します。レイヤーマネージャのメニューから[データベース]-[データベース接続]-[ベクトルマップーデータベース接続設定]を選択します注 6-2-4.。

[データベース]-[データベース接続]-[ベクトルマップ－データベース接続設定]を選択

図６－２－１１．v.db.connect コマンドの選択

選択すると図６－２－１２のようなダイアログが表示されます。"必須"タブには属性テーブルを接続させたいベクタレイヤーを指定します。

注6-2-4.GRASS コマンド v.db.connect と同値。

133

接続対象のレイヤーを指定

図６－２－１２．v.db.connect 必須タブ

続いて"オプション"タブを開きます。ここの"テーブル名"に既存のテーブル（ここではバッファ発
生元の観測点レイヤーのテーブル）を指定し、"キー列名"には"cat"を指定します。すなわち、"cat"フィールド（要素番号に対応）をキーとしてテーブルと要素を接続します。

既存のテーブルを選択

接続時のキーを指定

図６－２－１３．v.db.connect オプションタブ

134

設定が完了したら「実行」ボタンを押します。これでバッファポリゴンと元の観測点の属性値が接続されました。あとはラスタライズ時に"class"の値を属性値として使用すればトレーニングエリアのラ
スタデータが作成されます（ラスタライズの詳細は２－１章を参照してください）。

図６－２－１４．現地調査によるトレーニングエリア作成例

現地調査で得られたトレーニングエリアデータは最も信頼性が高いデータであるといえるので、４－４－１章で述べた精度評価に現地調査結果を利用すれば、現地検証結果との比較を行ったことと同義であるといえます。ただし、衛星データ撮影日と現地調査実施日が大きく離れている場合はその限りではないことに注意する必要があります。

135

６－３．被災地抽出結果を用いた踏査地点の出力

被災地抽出を行った結果を災害救助や復旧作業に用いることは被災地抽出作業の主目的であるといえます。抽出地点を地形図などの図面にオーバーレイすることで被災現場の位置を特定できますが、被災地が山奥であったり、または大規模な災害などで道路や建物などが損壊している状況では既存の地図では位置の特定が難しい場合も考えられます。そのような場合に GPS 機器にあらかじめ被災現場のおおよその位置を登録しておくことで、被災現場へのナビゲーションとして利用することができます。
ここでは QGIS を用いて被災地の代表点をデジタイズして取得し、GPX データに変換する簡単な方法を説明します。
まず、図６－３－１のように QGIS 上に AVNIR-2 の RGB 合成画像（VRT）と基盤地図情報の等高線、および被災地抽出結果が表示されているとします。この状態から被災地のおおよその地点を取得していきます。

図６－３－１．被災地抽出結果を表示した状態

まず、シェープファイルを新規作成します（QGIS 上でのシェープファイルの新規作成については４－２－１章参照）。新規作成ダイアログの"タイプ"の項目は"点"を選択し、「CRS の指定」ボタンでラスタレイヤーと同じ空間参照を選択します。さらに属性を追加しますが、ここでは仮に地名と災害種
別を属性として持たせることにします。

136

図６－３－２．新規ベクタレイヤー設定例

設定が完了したら「OK」ボタンを押してシェープファイルを作成しますが、日本語の属性値を入力し、なおかつこの後シェープファイルを GPX ファイルに変換する場合は、ファイルを保存するときのダイアログで"エンコード"を UTF-8 に設定する必要があります注 6-3-1.。

“エンコード”をUTF-8に設定する

図６－３－３．シェープファイル保存時のエンコード指定
注6-3-1.なお、フィールド名に日本語を使用することも可能ですが、GPX 変換時にエラーが出てしまいます。

137

作成されたら、キャンバス上から点を取得していきます。

「編集モード変更」
ボタンを押す「点を入力」
ボタンを押す

図６－３－４．ポイントシェープの編集開始

点の取得が完了したら、シェープファイルから GPX ファイルに変換します。FWTools のコマンドから、以下のように入力します。

C:\Data> ogr2ogr -f "GPX" -t_srs EPSG:4326 "出力 GPX ファイル名" "シェープファイル
名" ⏎上記のコマンドは GPX ファイルに変換を行いつつ、座標系を WGS84 経緯度に変換しています。入力するシェープファイルの空間参照が未定義の場合は、入力側の空間参照を指定する"-s_srs"オプションを指定するとよいでしょう。
作成された GPX ファイルを確認してみましょう。以下のような内容が出力されていれば成功です。

<?xml version="1.0"?><gpx version="1.1" creator="GDAL 1.7.0b2, FWTools 2.4.7"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xmlns="http://www.topografix.com/GPX/1/1"
xsi:schemaLocation="http://www.topografix.com/GPX/1/1
http://www.topografix.com/GPX/1/1/gpx.xsd"><metadata><bounds minlat="38.893122086476978"
minlon="140.766253072837150" maxlat="38.958999785579728"
maxlon="140.848041711914080"/></metadata><wpt lat="38.897930085189870" lon="140.848041711914080"><name>荒砥沢ダム</name><type>土砂崩れ</type></wpt><wpt lat="38.958999785579728" lon="140.812480251497250"><name>駒の湯温泉</name><type>土砂崩れ</type></wpt><wpt lat="38.893122086476978" lon="140.766253072837150"><name>湯の倉温泉</name><type>河道閉塞</type></wpt></gpx>

138

上記の例では各点をウェイポイントとして作成していますが、現場までのアクセスルートを地図上でデジタイズして GPX のルートレイヤーに変換することも可能です。

139

７．時系列データによる広域土地被覆評価
第４章や第５章の始めで述べたように、衛星画像データの特徴として観測範囲の広域性とデータ取得の周期性があります。なかでも NASA の地球観測衛星 TERRA と AQUA に搭載されている MODIS（中解像度画像放射計、MODerate resolution Imaging Spectroradiometer）は、観測幅が 2330km、地表面の全域の画像を１～２日毎という高頻度で取得しています。このように高頻度で観測される衛星画像
データは、周期的に変化する地表面の状況を検出するのに適しています。例えば、第３章で緑地の抽出に用いた NDVI（正規化植生指数）は、植生や土地利用の種類によって、特徴的な時系列変化を示します。
そこで本章では、この NDVI の時系列変化に基づき、広域の土地被覆を分類する方法を解説します。

140

７－１．対象地域と使用データの概要
解析対象地域は、中央アジアのカザフスタン、ウズベキスタン、トルクメニスタンを含む北緯 40 度から 47 度、東経 58 度から 65 度の範囲としました。この地域には縮小が問題になっているアラル海や、
それにそそぐアムダリヤ川、シムダリヤ川の河川沿いの灌漑地域、さらにカラクム砂漠キジルクム砂漠などの多様な土地被覆が含まれています。
解析に使用するデータは MODIS、MOD09Q1 プロダクトの画像データです。このデータは空間解像度 250m、8 日間の合成画像です。バンドは赤（620 - 670 nm）と近赤外（841 - 876nm）の 2 バンドからなります。このデータは NASA の LP DAAC（Land Processes Distributed Active Archive Center）からダウンロードが可能です。場所は h22v04 に含まれます。

図７－１－１

141

分析対象範囲７－２．決定木による土地被覆分類
決定木というのは、分類手法の一つで、画像だけでなく様々な場合に用いられます。画像分類の場
合もっとも単純なものは、第３章で示したように NDVI の値に閾値を設定し、ある値以上であれば緑地、それ以下であれば緑地ではないといった、二分法になります。
ここで、対象地域における代表的な土地被覆の NDVI の値を図７－２－１に示します。

図７－２－１

対象地域における典型的土地被覆の NDVI 値の変化

上図のように、水面の場合は通年通して 0 以下、砂漠などでは低い値を取りますが、農作物がある
場合は、特定の時期にピークを示したり、低い値を示します。例えば、小麦は春先に高い値を示しま
すが、綿花や水田は夏に高い値を示します。また水田は、６月初め頃に田面に水を張るため低い値を示します。このような特徴的な値に基づいて、分類木を作成し、各土地被覆を決定することとしました。

142

７－２－１．MODIS データの変換と GRASS へのインポート
MODIS のデータの変換は、昨年度のチュートリアルの４－４で示している方法で可能です。本分析では一年間のデータを参考にするため、47 枚と多くの画像を使用します。この変換を手動で行うと、単純作業であり、大変手間がかかります。
そこで、ここではスクリプトと呼ばれる簡単なプログラム言語を使って、変換とインポートを効率的に行う手順について説明します。使用するスクリプトは、以下のアドレスからダウンロードして下さい。
http://cse.niaes.affrc.go.jp/niwasaki/tutorial/script.zip
このファイルを解凍すると、中にはgrass_import.sh
hdf2img.shmodis_hdf2erdas_ll_wgs84.sh
ndvi_calc.sh
ndvi_class.sh
ndvi_series.shというファイルが入っています。これらのファイルの中身がスクリプトと呼ばれるファイルになります。
さて、まずは「GRASS 6.4.0 with MSYS」というアイコンをダブルクリックして GRASS を立ち上げます。

図 7-2-1-1 GRASS6.4.0 with MSYS のアイコン

すると、次の図の様な画面が立ち上がります。

143

図 7-2-1-2 GRASS6.4.0 with MSYS の起動画面

通常の起動画面に比べて、左上にクリーム色の背景の画面が出てきます。これが、MSYS の画面になります。ここから先は通常の GRASS での作業と同じように、ロケーションとマップセットを作成して下さい。ここでは、「C:\GRASS」というフォルダに「WGS84_LatLon」というロケーションと「NDVI」というマップセットを作成しました。以上の設定で GRASS を起動すると次のような画面になります。

図 7-2-1-2 GRASS6.4.0 with MSYS の起動画面

ここで、先ほどの GRASS フォルダの中に HDF というフォルダを作ってその中にダウンロードした
hdf ファイルをコピーします。また、同じく hdf フォルダの中にスクリプトをコピーして下さい。

144

その上で，MSYS の画面上で以下の様に入力して下さい。
cd /c/GRASS/hdf

ここで、CD とは、現在作業しているディレクトリの位置を変えるためのコマンドです。入力すると、先ほど作成した hdf フォルダに移動したことになります。ここで、
sh hdf2img.shと、入力して下さい。sh とはスクリプトを実行するために命令で、hdf2img.sh が実行するスクリプ
トになります。
このスクリプト実行されると自動的に hdf ファイルが変換され、MOD09Q1.A2010001.h22v04.b01.img
といった形式のファイルのが大量に作成されます。このうち最初の「MOD09Q1」は使用している
MODIS プロダクトの名前、「A2010001」はデータの日付、「h22v04」はタイルの名前、「b01」はレイヤー名になります。
この hdf2img.sh というスクリプトの中身は以下の通りです。
#!/bin/bash
LIST=(`ls *.hdf | tr '\n' ' '| sed 's/,$//g'`)for hdf_list in ${LIST[@]};doecho `sh ./modis_hdf2erdas_ll_wgs84.sh ${hdf_list}`
done

これは、フォルダの中にある全ての hdf ファイルについて、別のスクリプト
modis_hdf2erdas_ll_wgs84.sh を実行するというものです。この、modis_hdf2erdas_ll_wgs84.sh は、GDAL
コマンドを使って、hdf ファイルを緯度経度形式の img ファイルに変換するもので、GRASS の主開発者の一人である Markus Neteler 氏によって作成されたものです。
次に、
sh grass_import.shと入力して下さい。多少時間がかかりますが、これでインポートが完了します。
このスクリプトの中身は、以下の通りです。
#!/bin/bash
LIST=(`ls *.img | tr '\n' ' '`)for file_name in ${LIST[@]}; do
part=(`echo "$file_name" | tr '.' ' '`)r.in.gdal --overwrite input="$file_name" output="${part[1]}"."${part[3]}"
done

145

GRASS のコマンドである r.in.gdal コマンドを、全ての img ファイルについて行うというものになります。上記の変換作業とインポート作業は、それなりに時間がかかるかもしれませんが、スクリプトを使うことにより、作業の手間を減らすことが出来ます。
インポートされた衛星画像のカラーテーブルを Grey に設定し、赤バンドと緑バンドにバンド１、青にバンド２を設定して表示した例を図 7-2-1-3 に示します。

図 7-2-1-3 変換した画像の表示例（DOY225）

146

７－２－２．NDVI の計算
NDVI の計算方法については、本書の第３章で細かく書かれていますので、そちらをご参考下さい。
ただし、第３章であつかった LANDSAT では、バンド３(赤)とバンド４（近赤外）を使いましたが、MODIS ではバンド１が赤、バンド２が近赤外になるので、その点に注意して下さい。
まず、計算の範囲を設定します。メニューから「Config→Region→Set Region」を選択して下さい。

図 7-2-2-1 リージョンセッティングのメニュー

次に、境界タブをクリックして、リージョンの値を入力して下さい。ここでは、北端が 47、南端が40、東端が 65、西端が 58 になります。

図 7-2-2-2 範囲設定。

さらに、「解像度」タブをクリックし、グリッド解像度 2D に、「00:00:05.5」と入力して下さい。

147

図 7-2-2-3 解像度設定。

そして RUN をクリックすると、計算領域が設定されます。
次に、以下の様に入力して下さい。
sh ndvi_calc.sh

スクリプトの中身は、以下の通りです。
#!/bin/bash
DOY=2010001
while [ $DOY -le 2010365 ] # set datedor.mapcalc NDVI_A"$DOY" = "int(\
if(float(A"$DOY".b02 - A"$DOY".b01) / (A"$DOY".b02 + A"$DOY".b01) > 1, null(),\if(float(A"$DOY".b02 - A"$DOY".b01) / (A"$DOY".b02 + A"$DOY".b01) < -1,
null(),\
float(A"$DOY".b02 - A"$DOY".b01) / (A"$DOY".b02 + A"$DOY".b01))) * 10000)"
eval DOY=`expr "$DOY" + 8` # calc ODY
done

148

このスクリプトは、８日ごとの MODIS のバンド１のデータとバンド２のデータから NDVI を計算しています。ただし、NDVI は理論上-1 から 1 の間に収まるはずですが、時として主とした衛星画像のエラーがあるために、そこからはみ出す値を示す場合があります。そのため、ここでは、範囲外の値については除外しています。さらに、NDVI の計算は浮動小数点で行っていますが、データを保存する際にはサイズが大きくなるので、10000 倍して整数値として保存しています。
このプログラムを実施すると、合計 46 枚の NDVI 画像が生成されます。生成した NDVI の画像例は以下の通りです。アムダリヤ川とシルダリヤ川沿いに NDVI の高い地域が見られますが、これは灌漑農業地域です。また、全体として植被が少なく低い値を示し、そして縮小したアラル海を確認することが出来ます。

図 7-2-2-4 NDVI の分析例

149

７－２－３．決定木による画像分類の実施
本解析では、全期間での NDVI の全期間を通しての平均値や季節毎の最大値、最小値に基づき、決定木を用いた分析を行うこととしました。そのため、分析方法は以下のような二段階になります。
まず、NDVI の平均値、季節毎の最大値、最小値を算出します。次のようにコマンドを入力して下さい。
sh ndvi_series.sh

スクリプトの中身は以下の通りです。
#!/bin/bash

LIST=`g.mlist pattern=NDVI_A* | tr '\n' ','| sed 's+,$++g'`
r.series input=$LIST out=NDVI_MAX,NDVI_AVE method=maximum,average,minimumr.series input=NDVI_A2010105, NDVI_A2010113, \NDVI_A2010121, NDVI_A2010129, \NDVI_A2010137 method=maximum \out=NDVI_spring_max

r.series input=NDVI_A2010201,NDVI_A2010209, \NDVI_A2010217,NDVI_A2010225, \NDVI_A2010233,NDVI_A2010241\

method=maximum out=NDVI_summer_max

r.series input=NDVI_A2010129,NDVI_A2010137, \NDVI_A2010145,NDVI_A2010153, \NDVI_A2010161,NDVI_A2010169, \NDVI_A2010177,NDVI_A2010185\

method=minimum out=NDVI_june_minr.mapcalc "NDVI_diff_SprSum = NDVI_summer_max - NDVI_spring_max"
r.mapcalc "NDVI_ratio_SprSum = int((NDVI_diff_SprSum / NDVI_spring_max) * 100)"r.mapcalc "NDVI_ratio_JunSum = int(((NDVI_summer_max - NDVI_june_min) / \NDVI_june_min) * 100)"

こでは、r.series というコマンドを使って、全期間や季節毎の平均値、最大値、最小値を算出するとともに、さらに季節毎の値の比を計算しています。計算結果のうち、全期間の平均値および最大値を150

以下に示します。

図 7-2-3-1

全期間を通しての NDVI の平均値

図 7-2-3-2

全期間を通しての NDVI の最大値

実際の計算の実施に当たっては、以下の様に入力して下さい。
sh ndvi_class.sh

スクリプトの中身は、以下の通りです。
#!/bin/bash

r.mapcalc LANDCOVER = "(if(NDVI_AVE < 1, 0, \if(NDVI_MAX < 2000, 1, \151

if(NDVI_AVE < 2000, 2, \if(NDVI_MAX < 4000, 3, \if(NDVI_diff_SprSum < 0, if(NDVI_spring_max > 4000, 4, 5),\if(NDVI_ratio_SprSum >= 50, \if(NDVI_june_min < 1000, 6, if(NDVI_ratio_JunSum > 500, 7,8)), \
if(NDVI_june_min < 1000, 9, if(NDVI_ratio_JunSum > 400,
10, 11)))))))))"

このスクリプトを実行すると LANDCOVER という名前のファイルが生成されます。この画像はクラス
0 から 11 までに分類されています。

図 7-2-3-3

分類木を用いた時系列画像分類結果

ここでは、まず平均値が 0 以下の場合はクラス 0 として、分類しています。この場合、水面に該当する部分になります。次に、最大値が 2000 以下、平均値が 2000 以下、そして最大値が 4000 以下の場
合をそれぞれクラス 1，2，3 としました。最大値が 2000 以下の場合は砂漠など、平均値が 2000 以下や最大値が 4000 以下の地域は植皮の少ない灌木地帯、市街地などが該当します。
クラス 4 と 5 は春期の NDVI から夏期の NDVI を引いたときに、負の値を取るクラスになります。そのうち、春期の最大値が 4000 より大きいものはクラス４で、これは小麦の一毛作などが考えられます。
4000 以下のものはクラス 5 になり、湿地や森林などの自然植生が考えられます。
クラス 6、7，8 は春期と夏期の比が 50 以上のもの、つまり夏期に一つのピークを持つ場合で、クラス 6 と 7 が 7 月の初旬に低い値を持つ、すなわち水田等を示し、クラス 8 は綿花などを含むと考えられます。
クラス 9，10，11 は，春期と夏期の最大値の比が 50 以下のもの、すなわち春と夏に二つのピークを152

持つ二毛作で、9 と 10 は小麦の水田の二期作、11 は小麦とその他の作物の二期作が考えられます。
ただしこの分類が妥当なものであるかどうかは、現地踏査や、既存の土地利用図などにより確認する必要があります。

153８．参考文献●

参考書籍

粟屋善雄，田中邦宏，2003．ランドサット TM データを用いた森林変化モニタリング指標の検討：
スギ林の伐採と成長のモニタリング，写真測量とリモートセンシング，Vol．42，no．5，pp．60-70．
C. Bouman and M. Shapiro 1994, "A Multiscale Random Field Model for Bayesian Image Segmentation",
IEEE Trans. on Image Processing., 3(2), 162-177"
長谷川昌弘，今村遼平，吉川眞，熊谷樹一郎，2002．ジオインフォマティックス入門，理工図書，
東京，256pp．
飯倉善和，横山隆三，1998．ランドサット TM 画像の正射投影とその評価，写真測量とリモートセンシング，Vol．37，no．4，pp．12-22．
伊理正夫，腰塚武志，1986．計算機科学と地理情報処理，共立出版，東京，294pp．
菅野正人，2007．ALOS 衛星による下川町民有林の風倒被害解析と被害把握への活用，日本森林学会大会学術講演集（118）．
加藤正人，2007．改訂森林リモートセンシング，日本林業調査会，東京，358pp．
小泉圭吾，谷本親伯，朴 春澤，2003．LANDSAT ５号 TM データを用いた竹林の抽出手法に関する研究，写真測量とリモートセンシング，Vol．42，no．6，pp．42-51．
小浪尊宏，小野博之，大谷 徹，藤良太郎，2009．航空レーザ計測データを活用した災害対応の有効性と課題の検討，写真測量とリモートセンシング，Vol．48，no．6，pp．374-380．
Markus Neteler，Helena Mitasova 著：植村哲士
第３版 日本語版，開発社，東京，467pp．

訳，2009．オープンソース GIS グラス

アプローチ

日本地形学連合，1996．地形学から工学への提言，古今書院，東京，162pp．
日本リモートセンシング研究会，2001．改訂版
東京，325pp．

図解リモートセンシング，社団法人日本測量協会，

小倉久直，1978．物理・工学のための確率過程論，コロナ社，東京，266pp．
大津展之，栗田多喜夫，関田
P. G. Hoel 著：浅井巌，1996．パターン認識－理論と応用－，朝倉書店，東京，214pp．

晃，村上正康共訳，1963．初等統計学，培風館，東京，325pp．

島田静雄，2000．CAD・CG のための基礎数学，共立出版，東京，218pp．
財団法人画像情報教育振興協会，2006．ディジタル画像処理，財団法人画像情報教育振興協会，東京，382pp．
財団法人画像情報教育振興協会，1997．画像処理標準テキストブック，財団法人画像情報教育振興協会，東京，360pp．

●

ウェブページ153

Frank Warmerdam，2009．FWTools: Open Source GIS Binary Kit for Windows andLinux．http://fwtools.maptools.org/
GRASS Development Team，1999．GRASS GIS - The World Leading Free Software GIS．http://grass.itc.it/
Quantum GIS project，2002．Welcome to the Quantum GIS Project．http://qgis.org/
The Open Source Geospatial Foundation．GDAL - Geospatial Data Abstraction
Library．http://www.remotesensing.org/gdal/index.html
The Open Source Geospatial Foundation．PROJ.4 - Trac．http://trac.osgeo.org/proj/

154