とらりもんHOME  Index  Search  Changes  Login

高分解能衛星Landsat: 植生指標NDVI, 温度画像

筑波大学農林工学系 奈佐原顕郎

はじめに

 Landsatは、1970年代から継続的に、代替りしながら運用されている、高分解能光学パッシブ民生衛星である。センサーは、初期にはMSS, 中期にはTM, そして最近はETM+というのを搭載している。時間分解能は数十日と低いが、空間分解能は30m程度と高く、しかも可視3バンド、近赤外1バンドなどを持っているので、トゥルーカラー画像やフォールスカラー画像を構成できる。植生指数も計算できる。初期のMSSセンサー以外は、中間赤外・熱赤外バンドも持っているので、熱画像解析(都市ヒートアイランドや蒸発散研究)などにも使われる。

Landsatのページ(本家)

The Worldwide Reference System (WRS)とは? ... Landsatデータ用の、世界各地の位置コード。→ WRSと地図の対応

Landsatデータの取得

 Landsatデータの多くは, 下記のサイトから無料!!!で貰える。

メリーランド大学Global Land Cover Facility

FTPサイト: ftp://ftp.glcf.umd.edu/glcf/Landsat/

このサービスを利用して、Landsatのデータを解析してみよう。

関東地方周辺の、2001年9月23日のLandsat-ETM+画像をダウンロードする(バンド1,2,3,4のみ):

$ wget -c ftp://ftp.glcf.umd.edu/glcf/Landsat/WRS2/p107/r035/p107r035_7x20010924.ETM-EarthSat-Orthorectified/*nn[1-4]0.tif.gz

もし転送速度が遅かったら、ここからダウンロードしてもよい。

 ダウンロードが終わったら、.gzファイルを解凍すること:

$ for i in *tif.gz; do gunzip $i; done

GRASSの起動と、UTMでのLOCATION設定

この画像は、UTM投影法のzone54で表現されているので、それに対応したLocationでGRASSを起動しよう:

$ grass -text
(ここでProjection valuesを選ぶ)
LOCATION:   UTM54
MAPSET:     PERMANENT
Please specify the coordinate system for location      c 
Do you wish to specify a geodetic datum for this location?(y/n) [y]   y
Please specify datum name          wgs84
Now select Datum Transformation Parameters        1
Enter Zone:   54
Is this South Hemisphere (y/n) [n]  n

DEFINE THE DEFAULT REGION
NORTH EDGE: 10000
SOUTH EDGE: 0
WEST EDGE:  0
EAST EDGE:  10000
GRID RESOLUTION: 15
(このあたりは適当でよい。どうせregionは再設定するので)

GRASSが起動し、LOCATIONの設定が終わったら、データを読み込もう。このLandsatのデータは、GeoTIFFというファイル形式で記録されている。GeoTIFFは、もともと写真などの画像に使われていたTIFF形式を、地理情報データにも拡張し、位置などのヘッダー情報を自分自身に含んでいる。このような自己記述型のファイル形式は、ヘッダー情報などを手動で読み込んで解析しようとすると大変だが、それを処理してくれるソフトが整っていれば、逆に非常に簡便に、正確な読み込みやデータ管理ができる。

 GRASSは、GDALという、さまざまな衛星データのファイル形式(もちろんGeoTIFFも含む)に対応するライブラリとリンクしているので、GeoTIFF画像は、以下のコマンドでダイレクトに読み込むことができる。

 r.in.gdal input=p107r035_7t20010924_z54_nn10.tif output=B1
 r.in.gdal input=p107r035_7t20010924_z54_nn20.tif output=B2
 r.in.gdal input=p107r035_7t20010924_z54_nn30.tif output=B3
 r.in.gdal input=p107r035_7t20010924_z54_nn40.tif output=B4

(このoutput=以下のデータ名は、GRASSで扱うときの名前なので、各自で任意につけてよい。)

これらの4つのラスターデータを、それぞれ表示してみよ。

 d.mon x0
 g.region rast=B1
 d.rast B1

カラー合成

 衛星センサーは人間の目と同様に、様々な波長の光を見分けているので、それらを適当に光の3原色(青、緑、赤)にわりあててディスプレーに表示すれば、目で見たかんじの画像になる。Landsatはちょうど、バンド1が青、バンド2が緑、バンド3が赤に対応するので、そのように色を割り当てて合成表示してみよう。

 d.rgb blue=B1 green=B2 red=B3
2001_0924_ETM_123.jpg
 d.zoom
2001_0924_ETM_123_Tsukuba.jpg

このようなカラー合成(人の肉眼と同様の見えかた)を、トゥルーカラー画像と呼ぶ。

 一方、植生は近赤外光を強く反射するので、植生を特に強調して表示したいときは、バンド2を青、バンド3を緑、バンド4(近赤外)を赤に割り当てて表示することが多い:

 d.rgb blue=B2 green=B3 red=B4
2001_0924_ETM_234.jpg

このようなカラー合成(人の肉眼と同様の見えかた)を、フォールスカラー画像と呼ぶ。

放射輝度への変換

 以上のような処理は、単に「絵を眺める」ためならOKだが、定量的な解析をするには、実はまだ足りないことがある。そのひとつは、放射輝度への変換である。放射輝度とは、衛星センサーの各バンドが観測対象から受け取った光の強さを、実際の物理的な単位(例えばW/m^2/sr/micro_mなど)で表現するものと考えれば良い。それに対して、これまでの解析で扱ったのは、各バンドの光の強さを、255階調の整数値(8ビット非負整数)として相対的に表現したものであった(それをデジタルカウント値とかデジタルナンバーと呼ぶ)。

 デジタルカウント値を放射輝度に変換するには、通常、各衛星データの仕様を調べねばならない。その仕様がどこにどのように記載されているかは、衛星センサーごとに違うし、同じ衛星センサーでもデータの配布元や配布形式によって違ったりするので、ケースバイケースで判断する必要がある(そのためには衛星リモートセンシングや情報工学の基礎的な知識が必要である)。

参考になる: DNから放射輝度(Radiance)を求める方法 Landsat-5 TM ... 川村さんのページ

 さて、我々のLandsatデータについては、このような仕様は、FTPサイトに独立ファイルとして存在する。それをダウンロードして見てみよう:

wget -c ftp://ftp.glcf.umd.edu/glcf/Landsat/WRS2/p107/r035/p107r035_7x20010924.ETM-EarthSat-Orthorectified/p107r035_7x20010924.met
cat p107r035_7x20010924.met

すると、以下のような記述が見つかるだろう:

   ...
   LMAX_BAND1 = 191.600 
   LMIN_BAND1 = -6.200 
   LMAX_BAND2 = 196.500
   LMIN_BAND2 = -6.400 
   LMAX_BAND3 = 152.900 
   LMIN_BAND3 = -5.000 
   LMAX_BAND4 = 241.100
   LMIN_BAND4 = -5.100
   ...
   QCALMAX_BAND1 = 255.0            
   QCALMIN_BAND1 = 1.0             
   QCALMAX_BAND2 = 255.0            
   QCALMIN_BAND2 = 1.0               
   QCALMAX_BAND3 = 255.0              
   QCALMIN_BAND3 = 1.0            
   QCALMAX_BAND4 = 255.0              
   QCALMIN_BAND4 = 1.0           
   ...

この、LMAXというのが、デジタルカウント値の最大値が対応する放射輝度、LMINというのが、デジタルカウント値の最小値が対応する放射輝度である。で、QCALMAXというのがデジタルカウント値の最大値、QCALMINというのがデジタルカウント値の最小値である。ここでは、バンド1からバンド4まで、全てのバンドについて、デジタルカウント値の最大は255, 最小は1であることがわかる。(8ビット非負整数は、0から255までなので、実は0もデジタルカウント値として存在するが、それはnull値に使われている)

 たとえばバンド1の場合、デジタルカウント値の255は191.6W/m^2/sr/micro_mという放射輝度に対応し、デジタルカウント値の1は-6.2W/m^2/sr/micro_mという放射輝度に対応する。

 1と255の間のデジタルカウント値は、それぞれ線形に放射輝度に対応する。従って、デジタルカウント値(DN)を放射輝度(L)に変換する式は下記のようになる:

         (DN - QCALMIN) * (LMAX - LMIN)
   L = --------------------------------- + LMIN
   	 QCALMAX - QCALMIN

(c.f. Neteler and Mitasova, "OpenSource GIS: A GRASS GIS Approach 2nd edition" pp. 239-240.)

では、この式に基づいて、実際にデジタルカウント値(DN)を放射輝度(L)に変換してみよう:

r.mapcalc "B1L=(B1-1)*(191.6+6.2)/(255.0-1.0)-6.2"
r.mapcalc "B2L=(B2-1)*(196.5+6.4)/(255.0-1.0)-6.4"
r.mapcalc "B3L=(B3-1)*(152.9+5.0)/(255.0-1.0)-5.0"
r.mapcalc "B4L=(B4-1)*(241.1+5.1)/(255.0-1.0)-5.1"

この放射輝度に白黒のカラーテーブルを適用し、上と同様に、トルーカラー画像やフォルスカラー画像を表示してみよ。DNで表示したときとは、微妙に色合いが違って見えるだろう(もちろん、放射輝度で表示するほうが正しい):

r.colors map=B1L color=grey
r.colors map=B2L color=grey
r.colors map=B3L color=grey
r.colors map=B4L color=grey
d.rgb blue=B1L green=B2L red=B3L

植生指標

 衛星画像のいくつかのバンドを組み合わせると、植生の量や活性を反映するような指標を計算することができる。それを植生指標と呼ぶ。植生指標には様々なものがあるが、ここでは、一般に最もよく用いられる、正規化差分植生指標(Normalized Difference Vegetation Index; NDVI)を計算してみよう。NDVIは、次式で定義される:

            NIR - Red
   NDVI = -------------
            NIR + Red

ここでNIRは近赤外光の反射率または放射輝度(Landsatならバンド4)、Redは赤の反射率または放射輝度(Landsatならバンド3)である。

 それでは、前節で放射輝度に変換したLandsatデータを用いて、GRASS上でNDVIを計算してみよう:

r.mapcalc "NDVI = (B4L-B3L)/(B4L+B3L)"

NDVIは、原理的には-1から+1の間の値をとるが、マイナスになるのは水や氷など、植生がほとんど全く無いようなときだけである。それを勘案して、適当なカラーテーブルを用意しよう:

vi color_NDVI.txt
-1000 black
-0.1 black
0 red 
0.2 yellow 
0.5 green
0.7 blue
1.2 black
1000 black 
nv black
r.colors map=NDVI color=rules < color_NDVI.txt
d.rast NDVI
Landsat_NDVI.png
図: 植生指標(NDVI)の分布図。赤が低く、黄、緑、青の順に高くなる。9月下旬なので、川沿いの水田は収穫が終わっているのでNDVIは低い。筑波山などの森林のNDVIはまだ高い。

地表面温度

 以上は、可視・近赤外光に関する衛星データの処理であったが、Landsatは熱赤外線に関するデータも取得している。熱赤外線を使うと、地表面の温度(300K程度の温度)の分布を知ることができる。なぜならば、衛星で観測される熱赤外線のほとんどは、地表から熱放射されるものであり、その強度は地表面温度を強く反映するからである(プランクの法則)。

熱赤外線に相当するバンドは、Landsatではバンド6である。このデータをダウンロードして、上と同様に放射輝度に変換しよう:

wget -c ftp://ftp.glcf.umd.edu/glcf/Landsat/WRS2/p107/r035/p107r035_7x20010924.ETM-EarthSat-Orthorectified/*nn61.tif.gz
gunzip p107r035_7k20010924_z54_nn61.tif.gz
r.in.gdal input=p107r035_7k20010924_z54_nn61.tif output=B6
cat  p107r035_7x20010924.met

   ... LMAX_BAND61 = 17.040 LMIN_BAND61 = 0.000 ...


r.mapcalc "B6L=(B6-1)*(17.04-0.0)/(255.0-1.0)+0.0"

 さて、バンド6の放射輝度を地表面温度(摂氏)に変換するには、以下の式を使う:

r.mapcalc "TB6=1282.71/log(666.09/B6L + 1)-273.15"

(c.f. Neteler and Mitasova, "OpenSource GIS: A GRASS GIS Approach 2nd edition" pp. 239-240.)

これに、適当なカラーテーブルを当てて、表示してみよう:

vi color_temperature.txt
-10000 black
0 black
15 blue
20 green
25 yellow
30 red
40 black
10000 black
nv black 
r.colors map=TB6 color=rules < color_temperature.txt
d.rast TB6
Landsat_TB.png
図: 地表面温度の分布図。青が低く、緑、黄、赤の順に高くなる。水面(湖面や海面、川など)や森林で温度が低く、都市(土浦やつくば)や収穫後の水田などで温度が高い。

注意: 厳密には、Landsatなどの衛星センサーで得られた地表面温度を、そのまま信用することはできない。なぜなら、地表面温度と放射輝度の関係には、射出率というものが介在するからである。水やよく繁った植生は射出率はほとんど1.0であり、衛星センサーで得られた地表面温度は精度が高いが、裸地面や金属面などの射出率は低いことが多く、その場合、衛星センサーで得られた地表面温度は実際よりも低くなる。また、大気によるノイズなども、衛星センサーで得られた地表面温度に誤差を生じさせる。

UTM座標系

 UTM(Universal Transverse Mercator)座標系は、地球の北極から南極にむけて、まるでみかんの皮をむくように、細長い短冊状に地球表面をはぎとって平面に置いたような図法である。経度について6度づつ分割し、60領域(zone)で地表面を構成する。zoneは、日付変更線のすぐ東から始まって、アメリカ・ヨーロッパ・アフリカ・アジア・極東アジアの順に番号がついている。つくば周辺はzone54である。

 UTMは、Landsatクラスの分解能(数10m程度)か、それよりこまかい衛星画像について、非常によく用いられる図法である。対して、Landsatなどの中分解能か、それより荒い画像は、緯度経度直交座標系(Geographic)が用いられることが多い。

 これら2つの座標系の、それぞれの性質や規則を、よく覚えておこう。