[English page]

GRASSによるリモートセンシング画像処理:演習1(GUIでの操作)

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

演習内容:

  1. UTM座標系の衛星画像(Landsat ETM+)の表示、植生指数(NDVI)画像の作成GUIを用いた操作の説明
  2. 等経緯度座標系の標高データ(SRTM)の表示、モザイク合成、3次元表示
  3. ASTER画像からのNDVI、輝度温度の計算、標高データとの3次元表示
  4. 植生・水指数の時系列変化を用いたトンレサップ湖周辺およびメコン川氾濫域の抽出(牧雅康先生)

この演習1で最終的に得られる画像

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

GRASSのインストール

  1. Cygwinのsetup.exeファイルのダウンロード:http://cygwin.com/より
  2. setup.exeの起動:インストールに必要なファイルのダウンロード
    setup.exeをダブルクリック → 「Download Without Installing」 → 保存先の指定 → ミラーサイトの指定
  3. (ダウンロード終了後)setup.exeの起動:インストール
    setup.exeをダブルクリック → 「Install from Local Directory」 → インストールディレクトリの指定 → ダウンロードしたファイルのディレクトリの指定 → インストールするパッケージの指定(分からなければ全てを選択)

GRASSの起動

  1. Cygwinの起動:
    スタート → プログラム → Cygwin
  2. X-windowの起動:
    $ startx
  3. GRASSの起動:
    $ grass62

UTM座標系用のデータフォルダの作成

GRASS起動後にWindow右下の「Projection values」をクリック

GRASS起動後

  1. LOCATION, MAPSET, DATABASEの指定:
    LOCATION:  utm
    MAPSET:    PERMANENT
    DATABASE:  /home/student
    
    と入力したら、Escを押しながらEnter
    
  2. 座標系の指定:UTMを指定
    Please specify the coordinate system for location <utm>
    
    A   x,y
    B   Latitude-Longitude
    C   UTM
    D   Other Projection
    RETURN to cancel
    
    >C(と入力してEnter)
    
  3. 測地系の指定:
    Do you wish to specify a geodetic datum for this location? (y/n) [y] y
    
    Please specify datum name
    Enter 'list' for the list of available datums
    or 'custom' if you wish to enter custom parameters
    Hit RETURN to cancel request
    
    >list(と入力してEnter)
    
    (種々の測地系が表示される)
    Short Name  Long Name / Description
    
    (中略)
    tokyo   Tokyo       (bessl ellipsoid)
    wgs72   WGS_1972    (wgs72 ellipsoid)
    wgs84   WGS_1984    (wgs84 ellipsoid)
    
    Please specify datum name
    Enter 'list' for the list of available datums
    or 'custom' if you wish to enter custom parameters
    Hit RETURN to cancel request
    
    >wgs84(と入力してEnter)
    
    Now select Datum Transformation Parameters
    Please think carefully about the area covered by your data
    and the accuracy you require before making your selection.
    
    Enter 'list' to see the list of available Parameter sets
    Enter the corresponding number, or <RETURN> to cancel request
    
    >list(と入力してEnter)
    
    Number Details
    ---
    1    Used in whole wgs84 region
         (PROJ.4 Params towgs84=0.000,0.000,0.000)
         Default 3-Parameter Transformation (May not be optimum for older datums;
     use this only if no more appropriate options are available.)
    ---
    
    (中略)
    
    
    >1(と入力してEnter)
    
    Enter Zone: 53(と入力してEnter)
    
    Is this South Hemisphere (y/n) [n] n(と入力してEnter)
    
    
                      DEFINE THE DEFAULT REGION
    
    
                 ====== DEFAULT REGION ======
                 | NORTH EDGE: 4000000     |
                 |                          |
      WEST EDGE  |                          |EAST EDGE
      380000    |                          |630000
                 | SOUTH EDGE: 3700000     |
                 ============================
    
      PROJECTION: 1 (UTM)                   ZONE: 53
    
                        GRID RESOLUTION
                            East-West:     28.5
                          North-South:     28.5
    
    
    と入力した後で、Escを押しながらEnter
    
    
  4. → GRASSのWindowが起動する

画像ファイル(Landsat ETM+)の移動

/home/student/data/etmの下にある全てのファイル(p110r036_7t20011015_z53_nn30.tifなど)を、Windowsのエクスプローラーを使って/home/student/utm/PERMANENTの中に移動

※ETM+ (Enhanced Thematic Mapper Plus) 画像のダウンロード:環境リモートセンシングの演習では既にダウンロード済みでローカルディスクにコピーしてあるが、今後は下記のWebサイトからダウンロードするといい。

Global Land Cover Facility (GLCF) Earth Science Data Interface, University of Maryland

画像ファイル(Landsat ETM+)のインポートと表示

UTM座標系のLandsat ETM+画像を使用して演習を行う。

  1. 画像ファイル(TIFFフォーマット)のインポート:

    画像ファイルのインポート用の画面の起動

    次の3バンドのファイルに対し、下の画面の処理を繰り返す。


    画像ファイルのインポート(この例ではバンド3を処理している)

  2. インポートした画像の表示

    まずは、RGBの3層のレイヤーを指定し、redにb3、greenにb4、blueにb5を割り当てる。

    画像のRGB合成の指定

    次に、ディスプレイのウィンドウの上側にあるメニューバーの左端かその隣のアイコンをクリックして、画像を表示させる。

    ディスプレイのメニューバー:左端かその隣のアイコンをクリック

  3. 画像の拡大・縮小表示

    ディスプレイのウィンドウの上側にあるメニューバーの虫眼鏡のアイコンをクリックして、画像を拡大・縮小させる。

    画像の拡大・縮小表示

注意:上図のように画像が表示されないときは、「測地系の指定」での範囲の指定が間違っている可能性がある。その範囲を変更する手順を画像で示す手順を試すように。

画像の範囲の指定:手順1

画像の範囲の指定:手順2

画像の範囲の指定:手順3

上記の指定が終わったら、下図の左端のアイコンをクリックして新規にディスプレイを立ち上げて画像を表示させる。

新規にディスプレイの起動

画像の演算(NDVIの計算)

植生の粗密と相関のあるNDVI (Normalized Difference Vegetation Index) と呼ばれる植生指数を表示させる。Landsat ETM+画像の場合、バンド3と4から計算できる。

ただし衛星画像に記録されているのはDN (Digital Number) と呼ばれる単なる数値であり、DNから放射輝度、反射率を経由してNDVIを計算するのが一般的である。また放射輝度あるいは反射率は、そのままでは大気上端での放射輝度や反射率を意味するので、大気補正を実施してからNDVIを計算することもある。ここでは大気補正をせず、DNから放射輝度を求め、NDVIを計算する手順を説明する。

  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
    (中略)
    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
    (中略)
    END_GROUP = MIN_MAX_PIXEL_VALUE 
    
  2. DNから放射輝度への変換(GUIの方が簡単に指定可能):

    DN値が0の時は放射輝度(L3やL4, L5)を0に設定する条件処理を行う

    > r.mapcalc "L3=if(b3==0, 0, (b3-1)*(152.9+5.0)/(255.0-1.0)-5.0)"
    > r.mapcalc "L4=if(b4==0, 0, (b4-1)*(157.4+5.1)/(255.0-1.0)-5.1)"
    > r.mapcalc "L5=if(b5==0, 0, (b5-1)*(31.06+1.0)/(255.0-1.0)-1.0)"
    

    上の式をGUIで処理させるには、下記の手順に従って行う。

    手順1:「Map calculator」を選択

    演算処理1

    手順2:入力パラメータに例えば「b3」を、式の欄に「if(b3==0, 0, (b3-1)*(152.9+5.0)/(255.0-1.0)-5.0)」を入力する。

    演算処理2

    画像の端の領域が放射輝度が0になっているかを確認するには、画像を表示させた上で、 Icon of queryをクリックしてから画像上でクリックすれば、下図のように座標とRGBの値が表示される。

    画素値の表示

  3. 放射輝度からNDVIの計算:

    NDVIの計算

    一面緑色の意味のないような画像が表示されるだろう。しかし、画像上でクリックすれば-1から1までの値が入っていることが分かる。この範囲に合わせたカラーテーブルを割り当てよう。

  4. 色付けの変更:カラーテーブルを指定するウィンドウの起動

    NDVIのカラーテーブルの割り当て1

    NDVIのカラーテーブルの割り当て2

    上図の「Name of input raster map」で指定した「ndvi」は生成した画像を、「Name of predefined rules file」の「ndvi」は予め用意されたカラーテーブルを指定することを意味している。

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

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

参考にしたWebサイト


[演習1:UTM座標系の衛星画像(Landsat ETM+)の表示、NDVI画像の作成]
[演習2:等経緯度座標系の標高データ(SRTM)の表示、モザイク合成、3次元表示]
[演習3:ASTER画像からのNDVI、輝度温度の計算、標高データとの3次元表示]
[演習4:植生・水指数の時系列変化を用いたトンレサップ湖周辺およびメコン川氾濫域の抽出(牧雅康先生)]
[須崎純一のページへ]
須崎純一 京都大学大学院 工学研究科社会基盤工学専攻 空間情報学講座