ここでは、GRASSというGISソフトウェアを使った簡単な演習を説明しています。GRASSではリモートセンシング画像も表示、処理できます。このページは、GRASSのGUIでの操作を説明していきます。
[GRASSのインストール、標高データを用いた地滑り危険度マップの作成]
[植生指数 (NDVI) の計算、表示]
[標高データ (SRTM)の表示、植生指数 (NDVI) の3次元表示]
[反射率、輝度温度、標高データを用いた土地被覆分類]
[QGIS, Rを用いた公示地価データの空間統計分析：空間的自己回帰モデル]
[QGIS, Rを用いた公示地価データの空間統計分析：静的な時空間モデリング]
本講義用の衛星画像のファイル（p110r036_7t20011015_z53_nn30.tifなど）を、自分のGRASS用のディレクトリ（例えば、C:\GRASS_data\DEM\PERMANENTの中にコピーする。
※ETM+ (Enhanced Thematic Mapper Plus) 画像のダウンロード：本演習では既にダウンロード済みで上記URLにコピーしてあるが、今後は下記のWebサイトからダウンロードするといい。
QGISが既にインストールされているという前提で、話を進めていく。広範囲の衛星画像に対し、DEM画像で示された範囲に対応する衛星画像を切り抜いていく。
これでLandsat ETM+画像のバンド3, 4, 5, 61, 62のTIFF画像を切り抜いていきます。
|
|
|
| QGIS上で表示したDEM画像 | ポリゴン化して生成されたマスク画像 |
|
|
| クリッピング処理 | クリッピングされたLandsat ETM+画像 |
DEMの範囲に合わせて切り出したLandsat ETM+画像の、バンド3、4、5、6（61, 62の2種類）を、それぞれ「b3」「b4」「b5」「b61」「b62」という名前でインポートする。下記にコマンドを示すが、GUI上で処理する。ちなみに、ETM+画像のバンド3、4、5、6は、それぞれ緑、赤、近赤外、熱赤外の波長に対応している。
衛星画像に記録されているのはDN (Digital Number) と呼ばれる単なる数値である。DN自体には意味がなく、DNから変換される放射輝度、反射率は地表面の反射特性を表す物理量であり、地表面の種類や状態を推定するのに活用される。また、植生の密生度合いを反映する植生指数の一つであるnormalized difference vegetation index (NDVI) は放射輝度や反射率から計算するのが一般的である。また放射輝度あるいは反射率は、そのままでは大気上端での放射輝度や反射率を意味するので、大気補正を実施してからNDVIを計算することもある。ここではまずDNから放射輝度を求める手順を説明する。
/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
Landsat ETM+の場合、以下の式に従ってDNから（大気上端の）放射輝度へ変換する（一般的に下記の式のようにgainとoffsetを用いた線形式で記述される）。
上記の式をGRASS上で記述し、変換する。DN値が0の時は放射輝度（L3やL4, L5）を0に設定する条件処理を行う。
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」という関数を用いて処理する。
Normalized Difference Vegetation Index (NDVI) という植生の密生度合いを反映した指標を計算する。NDVIは-1から1の値を取り、植生が多いほど1に近い値を取る。
この処理後に、GRASSに既に用意された、NDVIに対応したカラーテーブルを適用する。
色付けの変更：
最終的に得られる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: 2hc2/λ5 (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モードで取得されたことが分かる。ファイルを開いて確認しよう。
色付けの変更：
| 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ではその前処理として、リサンプリングを行い、解像度を統一しておく必要がある。
まず、解像度の設定を確認する。
この結果、下記のように、DEMの解像度（11 m）が表示されるか確認する。
nsres: 11.0952222 ewres: 11.08798005 rows: 840
もし大きく異なる解像度が表示されれば
を使って、DEMやその派生画像を指定して、解像度を正しく設定する。
ここで、衛星画像の反射率データ（例えば、r3, r4, r5）をリサンプリングする。GRASSには複数のリサンプリング方法があるが、ここではbsplineを選んで処理してみる（あまり深い意味はない）。
この結果、r3に対してr3_resampleのようにリサンプリングした画像が得られたとして話を進める。
クラスタリングする前に、グループを作成する必要がある。下記ではサブグループも作り、slope, r3_resample, r4_resample, r5_resample, T61_resampleを指定する。
その後、下記のようにクラスタリングする前準備を行う。
この前準備の結果を用いて、分類する。
傾斜とr3, r4, r5（バンド3から5の反射率）、輝度温度を用いたクラスタリングの結果を下記に示す。画像が示すように、河川が分離できていない状態が陸域の市街地と分離できていない様子が分かる。
上記とは異なり、教師データを事前に与えて、それに従う分類手法、教師つき分類を実行する。その前に、今回もグループ、サブグループを作成する必要がある。教師データ、検証データ共に、画像上で特定の地域、Region of Interest (ROI) を指定する方法が簡便で効率的である。このROI選択時に反射率画像では識別しづらかったので、ここでは, b3_clip, b4_clip, b5_clipを指定する。例えば、グループ名、サブグループ名共に、「b345」とする。
その後、下記のようにROIを指定する画面を起ち上げる。
ここでは、5クラスを用意するとする。左下図のうち、赤線で囲んだアイコンをクリックすることで、5クラスの名称と色を指定できる。
左には、上側の画像で5クラスのROI（複数箇所）を選択した様子を示している。「Digitize new area」というアイコンをクリックし、範囲を選択し、右クリックでROIを決定できる。左図のうち、左側のグラフは「Histograms」を選択し、三角印のアイコンをクリックすると表示される。また下側の「Preview Display」では、「urban_results」等を選択することで、各クラスの分類結果のPreviewが表示される。
この前準備の結果を用いて、分類する。分類の良好さを画像で確認するために、「Name for output raster map holding reject threshold results:」にもファイル名を入力する。
分類を終えると分類結果画像とreject threshold results画像が表示されるが、分類結果はROIを指定した時の色合いと異なっている。そのため、下記のようにカラーテーブルを再度指定する。
上記のreject threshold results画像を使うことで、分類精度を大まかに把握することはできる。
下記の結果の内、例えば「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係数を計算するのが一般的である。今回、そこまで手が回らなかったが、意欲のある人はそのような定量的な検証結果を出してほしい。
指定日に各班の代表が発表します。それまでに各班で共同で作業をするように。
[GRASSのインストール、標高データを用いた地滑り危険度マップの作成]
[植生指数 (NDVI) の計算、表示]
[標高データ (SRTM)の表示、植生指数 (NDVI) の3次元表示]
[反射率、輝度温度、標高データを用いた土地被覆分類]
[QGIS, Rを用いた公示地価データの空間統計分析：空間的自己回帰モデル]
[QGIS, Rを用いた公示地価データの空間統計分析：静的な時空間モデリング]