反射率、輝度温度、標高データを用いた土地被覆分類

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

    演習内容:
  1. GRASSのインストール、標高データを用いた地滑り危険度マップの作成
  2. 反射率、輝度温度、標高データを用いた土地被覆分類
  3. QGIS, Rを用いた公示地価データの空間統計分析
    本演習内容:
  1. 使用する衛星画像
  2. QGISを使った衛星画像の切り出し
  3. 切り出した衛星画像のインポート
  4. 放射輝度の計算
  5. 反射率の計算
  6. NDVIの計算
  7. 輝度温度の計算
  8. 画像解像度の統一:リサンプリング
  9. 教師なし土地被覆分類(クラスタリング)
  10. 教師つき土地被覆分類(最尤法)
  11. 演習課題

使用する衛星画像

本講義用の衛星画像のファイル(p110r036_7t20011015_z53_nn30.tifなど)を、自分のGRASS用のディレクトリ(例えば、C:\GRASS_data\DEM\PERMANENTの中にコピーする。

※ETM+ (Enhanced Thematic Mapper Plus) 画像のダウンロード:本演習では既にダウンロード済みで上記URLにコピーしてあるが、今後は下記のWebサイトからダウンロードするといい。

QGISを使った衛星画像の切り出し

QGISが既にインストールされているという前提で、話を進めていく。広範囲の衛星画像に対し、DEM画像で示された範囲に対応する衛星画像を切り抜いていく。

  1. QGIS上で、DEM画像(TIFFフォーマット)の表示:QGISの左上の窓内でDEM画像ファイルをダブルクリック
  2. マスクの作成
    「ラスタ」→「変換」→「ポリゴン化」
    ※注意:「フィールド名:DN」にチェックを入れる。
  3. 衛星画像のマスク領域の抽出
    「ラスタ」→「抽出」→「クリッパー」

    これでLandsat ETM+画像のバンド3, 4, 5, 61, 62のTIFF画像を切り抜いていきます。

QGIS上で表示したDEM画像 ポリゴン化してマスクとして生成されたDEM画像
QGIS上で表示したDEM画像 ポリゴン化して生成されたマスク画像
クリッピング処理 クリッピングされたLandsat ETM+画像
クリッピング処理 クリッピングされたLandsat ETM+画像

切り出した衛星画像のインポート

DEMの範囲に合わせて切り出したLandsat ETM+画像の、バンド3、4、5、6(61, 62の2種類)を、それぞれ「b3」「b4」「b5」「b61」「b62」という名前でインポートする。下記にコマンドを示すが、GUI上で処理する。ちなみに、ETM+画像のバンド3、4、5、6は、それぞれ緑、赤、近赤外、熱赤外の波長に対応している。

「ファイル」→「ラスターデータのインポート」→「r.in.gdal」
r.in.gdal -o input=p110r036_7t20011015_z53_nn30_clip.tif output=b3
r.in.gdal -o input=p110r036_7t20011015_z53_nn40_clip.tif output=b4
r.in.gdal -o input=p110r036_7t20011015_z53_nn50_clip.tif output=b5
r.in.gdal -o input=p110r036_7t20011015_z53_nn61_clip.tif output=b61
r.in.gdal -o input=p110r036_7t20011015_z53_nn62_clip.tif output=b62

放射輝度の計算

衛星画像に記録されているのはDN (Digital Number) と呼ばれる単なる数値である。DN自体には意味がなく、DNから変換される放射輝度、反射率は地表面の反射特性を表す物理量であり、地表面の種類や状態を推定するのに活用される。また、植生の密生度合いを反映する植生指数の一つであるnormalized difference vegetation index (NDVI) は放射輝度や反射率から計算するのが一般的である。また放射輝度あるいは反射率は、そのままでは大気上端での放射輝度や反射率を意味するので、大気補正を実施してからNDVIを計算することもある。ここではまずDNから放射輝度を求める手順を説明する。

  1. メタファイル(ヘッダファイル)から変換係数の抽出:

    /home/student/utm/PERMANENTの下にある「p110r036_7x20011015.met.txt」というファイルをテキストエディタで(Wordpadでも何でも)開いてみる。すると、下記の記述を見つけるだろう。

    GROUP = MIN_MAX_RADIANCE 
    (中略)
    	LMAX_BAND3 = 152.900
    	LMIN_BAND3 = -5.000
    	LMAX_BAND4 = 157.400
    	LMIN_BAND4 = -5.100
    	LMAX_BAND5 = 31.060
    	LMIN_BAND5 = -1.000
    	LMAX_BAND61 = 17.040
    	LMIN_BAND61 = 0.000
    	LMAX_BAND62 = 12.650
    	LMIN_BAND62 = 3.200
    (中略)
    END_GROUP = MIN_MAX_RADIANCE 
    GROUP = MIN_MAX_PIXEL_VALUE
    (中略)
    	QCALMAX_BAND3 = 255.0
    	QCALMIN_BAND3 = 1.0 
    	QCALMAX_BAND4 = 255.0
    	QCALMIN_BAND4 = 1.0
    	QCALMAX_BAND5 = 255.0
    	QCALMIN_BAND5 = 1.0
    	QCALMAX_BAND61 = 255.0
    	QCALMIN_BAND61 = 1.0
    	QCALMAX_BAND62 = 255.0
    	QCALMIN_BAND62 = 1.0
    (中略)
    END_GROUP = MIN_MAX_PIXEL_VALUE 
    
  2. DNから放射輝度への変換:

    Landsat ETM+の場合、以下の式に従ってDNから(大気上端の)放射輝度へ変換する(一般的に下記の式のようにgainとoffsetを用いた線形式で記述される)。

    衛星画像のDNから放射輝度への変換式

    上記の式をGRASS上で記述し、変換する。DN値が0の時は放射輝度(L3やL4, L5)を0に設定する条件処理を行う。

    「ラスター」→「r.mapcalc」
    r.mapcalc "L3=if(b3==0, 0, (b3-1.0)/(255.0-1.0)*(152.9+5.0)-5.0)"
    r.mapcalc "L4=if(b4==0, 0, (b4-1.0)/(255.0-1.0)*(157.4+5.1)-5.1)"
    r.mapcalc "L5=if(b5==0, 0, (b5-1.0)/(255.0-1.0)*(31.06+1.0)-1.0)"
    r.mapcalc "L61=if(b61==0, 0, (b61-1.0)/(255.0-1.0)*(17.04-0.0)+0.0)"
    r.mapcalc "L62=if(b62==0, 0, (b62-1.0)/(255.0-1.0)*(12.650-3.2)+3.2)"

反射率の計算

Landsat ETM+に限らず、全ての衛星画像において、以下の式に従って大気上端の放射輝度から大気上端の反射率へ変換できる。大気上端 (top of atmosphere) ということで、TOA radiance, TOA reflectanceと表現されることもある。

ここでは、バンド3, 4, 5の反射率r3, r4, r5を求めてみる。

放射輝度から反射率への変換式

バンド Eλ (W/m2/μm)
1 1970
2 1842
3 1547
4 1044
5 225.7
7 82.06

上記のパラメータを利用して、放射輝度と同様に「r.mapcalc」という関数を用いて処理する。

「ラスター」→「r.mapcalc」
r.mapcalc "r3=3.14159*L3*0.997197^2/1547/cos(46.68*3.14159/180) "
r.mapcalc "r4=3.14159*L4*0.997197^2/1044/cos(46.68*3.14159/180) "
r.mapcalc "r5=3.14159*L5*0.997197^2/225.7/ cos(46.68*3.14159/180) "

NDVIの計算

Normalized Difference Vegetation Index (NDVI) という植生の密生度合いを反映した指標を計算する。NDVIは-1から1の値を取り、植生が多いほど1に近い値を取る。

「ラスター」→「r.mapcalc」
r.mapcalc "ndvi=if(r3==0 && r4==0, -1, (r4-r3)/(r4+r3))"

この処理後に、GRASSに既に用意された、NDVIに対応したカラーテーブルを適用する。

色付けの変更:

「ラスター」→「カラー調整」→「カラーテーブル r.colors」
r.colors map=ndvi rules=ndvi

最終的に得られるNDVI画像を下記に示す。

Landsat ETM+から計算されたNDVI画像

輝度温度の計算

衛星画像のDN値から輝度温度を計算するには、まず放射輝度へ変換する必要がある。ここでは、既にb61のDNデータを放射輝度L61に変換してあるという前提で話を進める。

輝度温度 (Tb):観測物体と等しい放射エネルギーを放射する黒体の温度

黒体 (black body):観測物体と等しい放射エネルギーを放射する黒体の温度

プランクの放射法則:放射輝度と輝度温度の関係
Lλ 黒体の分光放射輝度 (W/m2/str/μm)
Tb 輝度温度(黒体の絶対温度)(K)
λ 波長 (μm)
c 光速度: 2.988 × 108(m/s)
h プランク定数: 6.626 × 10-34(J s)
k ボルツマン定数: 1.380 × 10-23(J/K)

式(1)がプランクの放射法則と呼ばれ、それを変形した式(2)を用いることで、黒体の絶対温度つまり輝度温度が得られる。


Landsat ETM+の熱バンドの波長帯

下表のK1, K2の値は以下のWebページを参照した。

Landsat 7 Science Data Users Handbook

K1: 2hc25 (W/m2/str/μm) K2: hc/k λ (K)
Landsat 7 666.09 1282.71

注意: バンドごとにHigh gainモードで観測されたか、Normal gainモードかは異なる。p110r036_7x20011015.met.txtを確認すると、128・129行目に、BAND61はバンド6をLow gainモード、BAND62はHigh gainモードで取得されたことが分かる。ファイルを開いて確認しよう。

「ラスター」→「r.mapcalc」
r.mapcalc "T61=if(L61==0, 0, 1282.71/log(666.09 /L61 + 1))"

色付けの変更:

「ラスター」→「カラー調整」→「カラーテーブル r.colors」
r.colors map=T61 rules=gray

Landsat ETM+ b6から計算された輝度温度画像

放射率(emissivity)を考慮する場合

放射率を考慮した物体の温度と輝度温度の関係
Lλ ある物体から観測される分光放射輝度 (W/m2/str/μm)
Lλ,p ある物体が持つ分光放射輝度 (W/m2/str/μm)
ε 放射率 (無次元:0から1の値を取る)

式(3)のように、ある物体から観測される分光放射輝度はその物体が持つ分光放射輝度に放射率を掛けた値に等しいと仮定される。式(1)を式(3)に代入すると式(4)と書き換えられる。これをTbについて解けば、式(5)が導ける(Tについて解くことも可能)。

波長が長いとレイリー・ジーンズ近似 (Rayleigh-Jeans approximation) によって式(6)のように分光放射輝度を近似できる。式(6)を式(3)に代入すれば式(7)が導ける。

(レイリー・ジーンズ近似が成立するのは、hc/λkT << 1が成り立つ時である。例えばT=280(K)の時、λ>>50 μ mであり、マイクロ波以上の波長の領域である)


画像解像度の統一:リサンプリング

今回は10 m解像度のDEMやその派生画像と、28.5 mのLandsat ETM+画像を組み合わせて土地被覆分類を行う。GRASSではその前処理として、リサンプリングを行い、解像度を統一しておく必要がある。

まず、解像度の設定を確認する。

「設定」→「領域」→「領域の表示 g.region -p」

この結果、下記のように、DEMの解像度(11 m)が表示されるか確認する。

nsres:      11.0952222
ewres:      11.08798005
rows:       840

もし大きく異なる解像度が表示されれば

「設定」→「領域」→「領域設定 g.region」

を使って、DEMやその派生画像を指定して、解像度を正しく設定する。

ここで、衛星画像の反射率データ(例えば、r3, r4, r5)をリサンプリングする。GRASSには複数のリサンプリング方法があるが、ここではbsplineを選んで処理してみる(あまり深い意味はない)。

「ラスター」→「ラスターマップの作成」→「Resample using bspline r.resample.bspline」

この結果、r3に対してr3_resampleのようにリサンプリングした画像が得られたとして話を進める。

教師なし土地被覆分類(クラスタリング)

クラスタリングする前に、グループを作成する必要がある。下記ではサブグループも作り、slope, r3_resample, r4_resample, r5_resample, T61_resampleを指定する。

「画像」→「画像とグループの作成」→「グループの作成/修正 i.group」
注意:ウィンドウが表れたら、グループ名(例えば、「r345_T61_slope」)を入力した後、サブグループ名を入力する前に「追加」ボタンを押して、該当する画像を追加する。その後、「Edit/create subgroup」にチェックを入れて、サブグループ名(例えば、「r345_T61_slope」)を入力する。(先にサブグループ名を入力すると、画像を選択できなくなる)

その後、下記のようにクラスタリングする前準備を行う。

「画像」→「画像分類」→「教師無し分類用クラスタリング入力 i.cluster」

この前準備の結果を用いて、分類する。

「画像」→「画像分類」→「最尤法による分類(MLC) i.maxlik」

傾斜とr3, r4, r5(バンド3から5の反射率)、輝度温度を用いたクラスタリングの結果を下記に示す。画像が示すように、河川が分離できていない状態が陸域の市街地と分離できていない様子が分かる。

グループの作成 傾斜とr3, r4, r5(バンド3から5の反射率)、輝度温度を用いたクラスタリングの結果

教師つき土地被覆分類(最尤法)

上記とは異なり、教師データを事前に与えて、それに従う分類手法、教師つき分類を実行する。その前に、今回もグループ、サブグループを作成する必要がある。教師データ、検証データ共に、画像上で特定の地域、Region of Interest (ROI) を指定する方法が簡便で効率的である。このROI選択時に反射率画像では識別しづらかったので、ここでは, b3_clip, b4_clip, b5_clipを指定する。例えば、グループ名、サブグループ名共に、「b345」とする。

「画像」→「画像とグループの作成」→「グループの作成/修正 i.group」

その後、下記のようにROIを指定する画面を起ち上げる。

「画像」→「画像分類」→「Interactive input for supervised classification g.gui.iclass」
クラスの名称と色の指定

ここでは、5クラスを用意するとする。左下図のうち、赤線で囲んだアイコンをクリックすることで、5クラスの名称と色を指定できる。

  1. urban(人工構造物等)
  2. agriculture(農地)
  3. forest(森林)
  4. bareland(裸地:学校の校庭等)
  5. waterbody(水域)

ROIの指定

左には、上側の画像で5クラスのROI(複数箇所)を選択した様子を示している。「Digitize new area」というアイコンをクリックし、範囲を選択し、右クリックでROIを決定できる。左図のうち、左側のグラフは「Histograms」を選択し、三角印のアイコンをクリックすると表示される。また下側の「Preview Display」では、「urban_results」等を選択することで、各クラスの分類結果のPreviewが表示される。


この前準備の結果を用いて、分類する。分類の良好さを画像で確認するために、「Name for output raster map holding reject threshold results:」にもファイル名を入力する。

「画像」→「画像分類」→「最尤法による分類(MLC) i.maxlik」

最尤法による分類用のウィンドウ1 最尤法による分類用のウィンドウ2

分類を終えると分類結果画像とreject threshold results画像が表示されるが、分類結果はROIを指定した時の色合いと異なっている。そのため、下記のようにカラーテーブルを再度指定する。

カラーテーブルの変更

最尤法による分類結果 最尤法による分類結果のreject threshold results画像

上記のreject threshold results画像を使うことで、分類精度を大まかに把握することはできる。

「ラスター」→「レポートと統計」→「ラスターマップをカテゴリ別に集計 r.report」
注意:r.reportのウィンドウ中の「Statistics」タブの中の、「area in square meters」「percente cover」にチェックを入れる。

下記の結果の内、例えば「90%」は90%の信頼区間を表す。

r.report map=mlm_b345_rst_reject@PERMANENT units=me,p                           
+------------------------+
|                         RASTER MAP CATEGORY REPORT |
|LOCATION: data             Mon May 16 16:48:29 2016|
|------------------------|
|        north: 3873299.99812    east: 568508.733604 |
|REGION  south: 3863980.01147    west: 557032.674252 |
|        res:      11.0952222    res:    11.08798005 |
|----------------------|
|MASK: none            |
|----------------------|
|MAP: Rejection Probability for mlc_b345_rst 
 (mlc_b345_rst_reject@PERMANENT in|
|----------------------|
| Category Information |     square|   %  |
| #|description        |     meters| cover|
|----------------------|
| 1|0.1% . . . . . . . |  7,695,249|  7.19|
| 2|0.5% . . . . . . . |  5,218,907|  4.88|
| 3|1% . . . . . . . . |  3,263,201|  3.05|
| 4|2% . . . . . . . . |  4,374,473|  4.09|
| 5|5% . . . . . . . . |  7,558,939|  7.07|
| 6|10%. . . . . . . . |  8,558,752|  8.00|
| 7|20%. . . . . . . . | 11,876,206| 11.10|
| 8|30%. . . . . . . . |  9,474,909|  8.86|
| 9|50%. . . . . . . . | 15,200,796| 14.21|
|10|70%. . . . . . . . | 13,275,969| 12.41|
|11|80%. . . . . . . . |  6,288,967|  5.88|
|12|90%. . . . . . . . |  6,416,173|  6.00|
|13|95%. . . . . . . . |  3,164,290|  2.96|
|14|98%. . . . . . . . |  1,987,815|  1.86|
|15|99%. . . . . . . . |    667,280|  0.62|
|16|100% . . . . . . . |    734,205|  0.69|
| *|no data. . . . . . |  1,200,587|  1.12|
|----------------------|
|TOTAL                 |106,956,720|100.00|
+----------------------+

本来は、検証用のROIを用意して、得られた分類結果の精度を評価しなくてはならない。講義で学習したように、User's accuracy, Producer's accuracy, Overall accuracyに加えて、Kappa係数を計算するのが一般的である。今回、そこまで手が回らなかったが、意欲のある人はそのような定量的な検証結果を出してほしい。

演習課題

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

注意:出力結果は、エクスポート機能を使って、画像として保存し、PPTファイルに貼り付けるように。
  1. 教師つき土地被覆分類(最尤法)では、画像輝度値(Digital Number: DN)をそのまま使用した。DNではなく、バンド3〜5の反射率を用いて同様に最尤法で分類せよ。DNを用いた結果と比較した違いを述べよ(上記のように定量的に示せなくても、違いが顕著な地域を明示する等、何らかの方法で示すように)。
  2. 上記のデータに加えて、温度だけ、温度と傾斜を加えた場合の分類結果との違いを確認せよ。
  3. 下記のように、同一範囲(同一パス、ロウ)において複数時期のLandsat TM, ETM+画像が利用できる。複数時期の画像を使ってできる土地被覆分類を考えて、その結果を発表せよ。

[GRASSのインストール、標高データを用いた地滑り危険度マップの作成]
[反射率、輝度温度、標高データを用いた土地被覆分類]
[QGIS, Rを用いた公示地価データの空間統計分析]
[須崎純一のページへ]
須崎純一 京都大学大学院 工学研究科社会基盤工学専攻 空間情報学講座