[English page]

GRASSによるリモートセンシング画像処理:演習2

演習内容:

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

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

Landsat ETM+とSRTM画像から作成された鳥瞰図
前回まではUTM座標系のETM+画像を使用してきたが、ここから先は等経緯度座標系のSRTMという標高データを用いた話を展開する。

等経緯度座標系用のデータフォルダの作成

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

  1. LOCATION, MAPSET, DATABASEの指定:
    LOCATION:  latlon
    MAPSET:    PERMANENT
    DATABASE:  /home/student
    
    と入力したら、Escを押しながらEnter
    
  2. 座標系の指定:Latitude-Longitudeを指定
    Please specify the coordinate system for location <utm>
    
    A   x,y
    B   Latitude-Longitude
    C   UTM
    D   Other Projection
    RETURN to cancel
    
    >B
    
  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
    
    (種々の測地系が表示される)
    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
    
    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
    
    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
    
    
                      DEFINE THE DEFAULT REGION
    
    
                 ====== DEFAULT REGION ======
                 | NORTH EDGE: 36          |
                 |                          |
      WEST EDGE  |                          |EAST EDGE
      134    |                          |136
                 | SOUTH EDGE: 34          |
                 ============================
    
      PROJECTION: 3 (Latitude-Longitude)      ZONE: 0
    
                        GRID RESOLUTION
                            East-West:     0.0008
                          North-South:     0.0008
    
    
    と入力した後で、Escを押しながらEnter
    
  4. → GRASSのWindowが起動する

標高データの移動:SRTM (Shuttle Radar Topography Mission) 画像

/home/student/data/srtmの中にある全てのファイルを/home/student/latlon/PERMANENTの中に移動

注意:圧縮ファイルのままで、解凍する必要はない

※Windowsのエクスプローラーを使っても、あるいはCygwinのWindow内でコマンドで指定しても可能

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

CGIAR Consortium for Spatial Information, 90m DEM

標高データのインポートと表示

等経緯度座標系の標高データのSRTM画像を使用して演習を行う。

  1. 作業ディレクトリの変更:画像ファイルのあるディレクトリへ移動(GUIを使用する場合には不要)
    > cd /home/student/latlon/PERMANENT
    
  2. 画像ファイルのインポート:N35E135.hgt.zipファイルが必要(GUIの方が簡単に指定可能)
    > r.in.srtm input=N35E135  output=N35E135
    
  3. 画像ファイルの表示(GUIの方が簡単に指定可能)
    > d.mon x0
    > d.rast N35E135
    
    > r.colors map=N35E135 rules=srtm(SRTM用のカラーテーブルの読込)
    > d.rast N35E135
    
    ※N35E135以外にも、N35E134・N34E134・N34E135に対しても同様の処理を行ってみよう。

異なる範囲の複数の画像の合成(モザイク合成)

  1. 画像ファイルのインポート(「標高データの表示」で既にインポートしてあればこの処理は不要)
    > r.in.srtm input=N34E134  output=N34E134
    > r.in.srtm input=N34E135  output=N34E135
    > r.in.srtm input=N35E134  output=N35E134
    > r.in.srtm input=N35E135  output=N35E135
    
  2. 画像範囲の指定(「等経緯度座標系用のデータフォルダの作成」で北緯34度〜36度、東経134〜136度を指定していれば、この処理は不要)
    > g.region w=134 e=136 n=36 s=34 res=0.0008
    
  3. モザイク合成
    > r.mapcalc "area1=if(isnull(N34E134), -1, N34E134)"
    > r.mapcalc "area2=if(isnull(N34E135), -1, N34E135)"
    > r.mapcalc "area3=if(isnull(N35E134), -1, N35E134)"
    > r.mapcalc "area4=if(isnull(N35E135), -1, N35E135)"
    
    > r.mapcalc "areas=max(area1, area2, area3, area4)"
    
    > r.null map=areas setnull=-1
    
    > d.rast areas
    
    (SRTM用のカラーテーブルの読込)
    > r.colors map=areas rules=srtm
    > d.rast areas
    

投影法の変換:ETM+画像をUTM座標系から等経緯度座標系へ

  1. GRASSを立ち上げてUTM座標系を選択してから
    一度でもETM+画像のインポート(画像ファイル(Landsat ETM+)のインポートと表示で処理が終わっていれば、この処理は不要) )
    > r.in.gdal input=p110r036_7t20011015_z53_nn30.tif output=b3
    > r.in.gdal input=p110r036_7t20011015_z53_nn40.tif output=b4
    > r.in.gdal input=p110r036_7t20011015_z53_nn50.tif output=b5
    
  2. ※一度インポートしたら、/home/student/utm/PERMANENT/cell の中に記録されている。

  3. 一旦GRASSを閉じて、再度立ち上げて等経緯度座標系を選択してから
    画像範囲の指定(「等経緯度座標系用のデータフォルダの作成」で北緯34度〜36度、東経134〜136度を指定していれば、この処理は不要)
    > g.region w=134 e=136 n=36 s=34 res=0.0008
    
  4. 等経緯度座標系を選択したままで
    投影法の変換
    > r.proj input=ndvi location=utm mapset=PERMANENT output=ndvi_lat
    
    > d.mon x0
    > d.rast ndvi_lat
    
    ※ndvi以外にも、b3・b4に対しても同様の処理を行ってみよう。

標高とNDVI画像の重ね合わせ

  1. 等経緯度座標系を選択した状態で
    > d.mon x0
    > d.rast areas
    > d.rast -o ndvi_lat
    

3次元表示:鳥瞰図作成

  1. NVIZの起動(GUIの方が簡単に指定可能)
    > nviz elevation=areas color=ndvi_lat
    

    視線方向などの諸設定はGUI上で行う。

  2. 最終的に得られる鳥瞰図の一例を下記に示す。高さはSRTMの標高データを表し、表面の色はNDVIの値を表している。NDVI画像の範囲が限られているため、NDVI画像が無い地域は表面が灰色で示されている。

    Landsat ETM+とSRTM画像から作成された鳥瞰図の一例

演習:余力のある人へ

今回は、画像範囲を北緯34度〜36度、東経134〜136度に指定した(「等経緯度座標系用のデータフォルダの作成」)。そのため、Landsat ETM+画像の範囲はそれよりも広いにもかかわらず、一部が切り取られていることが上記の鳥瞰図でわかる。

そこで、対象範囲を北緯34度〜36度、東経134〜137度に変更して鳥瞰図を作成してみよう(SRTM画像ファイルは既にコピーしてある)。

参考にしたWebサイト


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