衛星・GIS データ処理の基礎

衛星・GIS データ処理の基礎演習

海洋研究開発機構
小林秀樹
２０１9 年 4 月 21 日版

11

この文書は、クリエイティブ・コモンズライセンス CC BY NC SA に従う範囲で自由にお使いいただけます。

- 1 -衛星・GIS データ処理の基礎

目次
１章

はじめに

２章

QGIS によるデータの読み込みと表示
2.1 衛星データの仕様の確認
2.2 QGIS の起動と準備
2.3 衛星データの読み込みと表示
2.4 トゥルーカラー画像の作成
2.5 基盤地図情報の取得（WMS）

３章

反射率と植生指数の計算
3.1 反射率の算出
3.2 植生指数の計算

４章

地図投影変換と標高データ（DEM）の利用
4.1 ASTER GDEM データの読み込み
4.2 地図投影変換
4.3 地形解析

５章

ベクタデータを用いた解析
5.1 新規ベクタレイヤの作成
5.2 ポリゴン内の平均 NDVI と平均標高の計算
5.3 河川周辺の植生解析（バッファー解析）- 2 -衛星・GIS データ処理の基礎

１．

はじめに

本実習では GIS ソフトを利用して衛星データ処理を行います。利用するソフト QGIS はオープンソースのソフトです。QGIS は主要なオペレーティング・システムに対応しているため、ほとんどのコンピュータ環境で利用することができます。２０１9 年４月現在、最新の安定バージョンは 3.46 です。本実習では
QGIS3.4 の利用を前提として話を進めますがほとんどの機能は最新バージョンでも機能するはずです。
QGIS は GIS 処理のためのオープンソースの地理情報システムです。ソフトウ
エアは以下のサイトからダウンロードできます。Windows, Mac, Linux など複数の環境に対応しています。
http://www.qgis.org/ja/site/
本実習では、Landsat-8（ランドサット８号）の地表面反射率画像を利用して以下の処理を行います。

- 3 -衛星・GIS データ処理の基礎

２．

QGIS によるデータの読み込みと表示

2.1 衛星データの仕様の確認

データフォルダの中の、拡張子が xml ファイルを検索して、衛星データのメタ情報を取得する。

・ 衛星センサの名称：
・ 観測日：
・ 観測時間（日本標準時）：
・ 観測時の太陽天頂角：・ 観測時の太陽方位角：・ 測地系：
・ 投影法：
・ 撮影範囲Ø 中心緯度経度：
Ø 画像四隅の緯度経度：
Ø 面積・東西南北の広がり：
・ 空間解像度：
・ 観測波長（色）：・ 衛星データのファイルフォーマット：
・ 画像サイズ（縦横ピクセル数）：（注意）：GMT=グリニッジ標準時刻- 4 -衛星・GIS データ処理の基礎

2.２ QGIS の起動と事前準備

QGIS を起動後、今回の衛星データ解析用に環境に設定する。
①

プロセッシングツールの表示
ツールバーから(英語) [Processing]-[Toolbox]* ウインドウの右に Processing Toolbox が表示される。Toolbox はのちほど地域統計計算を行うときに利用する。

②

マックユーザーで Python がうまく作動しない場合の対応[Setting]-[Option]-[System]

Environment Enable “Use Custom Variable”＋マークをクリックして Prepend を選択、さらに Variable に PATH と記載して、以下の path を追加する。
"/Library/Frameworks/GDAL.framework/Programs:/Library/Frameworks/Python.frame
work/Versions/3.6/bin:"- 5 -衛星・GIS データ処理の基礎

2.３ 衛星データの読み込みと表示① 画像の読み込み

ツールバーから(英語)[Layer]-[Add layer]-[Add Raster layer…]

で読み込むファイルを選択する。ここでは、Landsat 画像の青、緑、赤、近赤外の波長のデータを選択する。各バンドの Landsat 画像のファイル名は、以下のとおりである。
青バンド（band 2）：LC81070352014215LGN00̲sr̲band2.tif
緑バンド (band 3)：LC81070352014215LGN00̲sr̲band3.tif
赤バンド (band 4)：LC81070352014215LGN00̲sr̲band4.tif
近赤外バンド(band 5)：LC81070352014215LGN00̲sr̲band5.tif
画像がウインドウ左下のレイヤパネルに追加される。もし、レイヤパネルが表示されていない場合には、(英語) [View]-[panel]-[Layer panel]を選択しレイヤパネルを追加する。
画像ビューワには、チェックされている画像がグレースケースで表示されている。チェックを外すと非表示になる。４つの画像を一つ一つ表示して、それぞれの波長のデータを確認する。

- 6 -衛星・GIS データ処理の基礎② 画像ヒストグラムの確認

画像のプロパティを確認する。レイヤパネルで画像を一つ選択して右クリックし（Mac の場合、control+クリック）、サブメニューから[プロパティ]を選択。
左の項目を選択すると画像のプロパティを確認・変更できる。ここでは、[ヒストグラム]を選択して、画像ヒストグラムを確認する。青、緑、赤、近赤外のヒストグラムを比較する。

- 7 -衛星・GIS データ処理の基礎

このヒストグラムは画像中の各画素（ピクセル）のデジタル値（DN）の頻度分布である。
③ 画像ストレッチ

Layer panel で画像を選択して右クリックし（Mac の場合、control+クリック）、サブメニューから以下を選択する。
(英語) [Properties]̶[Symbology]
ここでは表示画像のコントラストを変更できる。表示されている画像は、このダイアログに表示されている最小と最大のデジタル値をグレースケールの濃淡にスケーリングされて表示されている。デフォルトでは、最小値と最大値で画層の濃淡を表現している。最小と最大の値を変えることで画像のコ- 8 -衛星・GIS データ処理の基礎

ントラストを変更できる。
④ 表示色の変更

表示画像はデフォルトでは、グレースケール（白黒濃淡）で表示される。
[Symbology]のレンダータイプを変更することで、表示色を変更できる。レ
ンダータイプのドロップタブを単バンド疑似カラー(Singlebandpseudocolor)に変更する。[新規カラーマップを作成]から表示色を変更する。
下の図の例では、Spectral を選択している。選択後、[Apply]ボタンをクリックする。

表示色が正しく変更された場合、以下のように画像が虹色で表示される。

- 9 -衛星・GIS データ処理の基礎- 10 -

衛星・GIS データ処理の基礎

2.４トゥルーカラー画像の作成

衛星画像を写真の様にカラー表示させる。QGIS では独立した画像レイヤをカラー表示することはできない。そのため、２．３節で読み込んだ４つの画像をひとまとめにしたマルチバンドデータを作成する。
①

マージツールの利用

ツールバーから（日本語）[ラスタ]-[その他]-[結合]（英語）[Raster]-[Miscellaneous]-[Merge]を選択。Landsat-8 では、青、緑、赤、近赤外の反射率データは 2.3 節で説明したように band2, band3, band4, band5 に対応する。入力ファイ
ルで青、緑、赤、近赤外の４つの波長に対応する入力ファイルを Input layerのとこで選択する。データは選択した順番で保存される。 さらに、[Please
each input file into a separate band] (別々のバンドに各入力ファイルを配置する)のチェックボックスをチェックする。Output data type はInt16(符号付き 16 ビット整数)を選択する。
出力先のファイル名は Advanced parameters の Merged のところで入力する。

- 11 -

衛星・GIS データ処理の基礎

マルチバンド画像はウインドウ内の GDAL/OGR console call に書かれている gdal スクリプトで作成される。Run をクリックすると、マルチバンド画像が作成されレイヤパネルに追加される。
②

カラーコンポジット画像の表示

レイヤパネル内のマルチバンド画像をチェックすると下のようなカラー画像が表示される。２．３の作業と同様にレイヤパネルでマルチバンド画像を選択し、右クリックでポップアップメニューから[Properties]を選択し、[Symbology]タブを選択する。マルチバンド画像の場合、Red band,Green band, Blue band の三色にそれぞれ、Band 1(青)、Band 2(緑)、Band 3(赤)が割り当てられている。ここで、Red band に Band 3(赤)、Blue band に Band 1(青)を指定し、それぞれのバンドの Min/max の値を調整した後 Apply をクリックすると、トゥルーカラー画像が表示される。

- 12 -

衛星・GIS データ処理の基礎

さらに、トゥルーカラー画像とは異なるデータを Red band, Green band,Blue band に割り当てて、フォルスカラー画像を作成してみる。Red bandに Band４(近赤外)、Green band に Band 3(赤)、Blue band に Band 2(緑)を割り当ててみる。フォルスカラー画像で人間の目にみえない波長の画像を可視化することで、トゥルーカラー画像では見えない地表面の特性を把握することができる。

- 13 -

衛星・GIS データ処理の基礎

2.５基盤地図情報の取得（WMS）

QGIS は、Web Map Service（WMS）に対応している。WMS はインターネットを利用して地図情報を取得するサービスである。ここでは、国土地理院が作成し、農研機構が配信サービスを提供する２５０００分の１地図を取得して衛星データに重ねて表示する。
ツールバーから[Layer]-[Add WMS/WMTS layer]を選択し、Data Source
Manager ダイアログボックスの Layers タブの上部にある NEW を選択する。
Name の欄に「Kiban２５０００」などわかりやすい名前を付ける。そして、URL の欄に以下の URL を入力する。して OK をクリックする。
http://www.finds.jp/ws/wms.php
Data Source Manager ダイアログボックスの Layers タブの上部にあるConnect をクリックすると基盤地図２５０００情報がダイアログボックスに表示される。ここで、ID１の FGD25000 を選択し、[add]をクリックすると、Layerpanel に FGD25000 が追加される。
Layer panel で FGD25000 を選択し、右クリックでサブメニューから[Properties]-[transparency]を選択し、透過率を５０％にすると、以下のようにフォールス画像に基盤地図が重なって表示される。

- 14 -

衛星・GIS データ処理の基礎- 15 -

衛星・GIS データ処理の基礎

３．
３．１

反射率と植生指数の計算
反射率の算出

！本編では地表面反射率のプロダクトの利用を前提としている。その場合は、３．１節は読み飛ばしてもよい。大気上端輝度を利用する場合には以下の処理が必要となる。また、２０１６年４月現在、Landsat の地表面反射率データの公開が始まっており、それを利用する場合には、以下の計算を行う必要はない。
２章で表示した画像は、地表の色を反射の強度に比例したデジタル値で表現していた。画像の定性的な確認であれば、デジタル値の画像でも問題ない。しかし、より定量的な解析を行うためには、各波長のデータを反射率などの物理量に変換する必要がある。デジタル値は、以下の式で大気上端の反射率に変換できる。

𝑅" = 𝜋𝐿" =%& ' ()*+,-./01345 6,789∗;8<'=>'?@

Rλ:反射率
DN：デジタル値
Cal:校正係数
Lλ:大気上端輝度（Wm-2 um str）
Bandwidth:センサの観測波長幅d:地球̶太陽の平均距離（天文単位 AU）
ESUN:観測波長における分光太陽放射照度
θs:太陽天頂角- 16 -（１）（２）

衛星・GIS データ処理の基礎

３．２

植生指数の計算

正規化植生指数（Normalized Difference Vegetation Index, NDVI）を計算し、その空間分布を把握する。NDVI は、以下の式で算出できる。
NDVI＝（近赤外反射率-赤反射率）/(近赤外反射率+赤反射率)(3)

画像演算には、Raster calculator を利用する。ツールバーから[Raster]-[Rastercalculator]をクリックすると、Raster calculator ダイアログが現れる。ここで、赤反射率と近赤外反射率のデータを選択して、式（３）を記述する。正しい計算式が入力できると、expression valid と表示される。ここで、output layer
名を入力し OK を押すと、計算結果が layer panel に追加される。
・ NDVI のレイヤを Single pseudocolor で表示する。
・ NDVI の空間分布を確認する。
・ ヒストグラムで NDVI の頻度分布を確認する。

- 17 -

衛星・GIS データ処理の基礎

４．

地図投影変換と標高データ（DEM）の利用

デジタル標高モデル（DEM）を読み込んで表示する。また、地形解析ツールを使って、傾斜角や斜面方位を計算する。DEM データはラスタデータであり、土地の標高が二次元のグリッドデータとして格納されている。今回は、Terra 衛
星の ASTER（アスター）センサのステレオ視観測で作成されたグローバル DEMデータ、ASTER GDEM version 2 を利用する。ASTER GDEM は以下の NASA
EARTH DATA サイトから無償でダウンロードできる。
https://search.earthdata.nasa.gov/
４．１ASTER GDEM データの読み込みツールバーから[layer]-[Add raster layer…]を選択し、DEM データ
Kanto̲DEM.tif を選択。
・ Layer panel で Kanto̲DEM をチェックしてデータを表示
・ NDVI や各波長のデータ、基盤地図データと重ねあわせて比較する。
４．2

地図投影変換現在公開されている ASTER GDEM データが採用している空間参照システ
ムは WGS 84 緯度経度座標(EPSG: 4326)である。データの空間参照システム
は、レイヤパネルの Kanto̲DEM データを右クリックして[Properties]を選択し、[Source] を選択することで確認できる。
一方で、Landsat 画像の空間参照システムを同様の方法で確認すると UTM
ゾーン 54N/WGS84 (EPSG:32654)と表示されている。QGIS では異なる空間参照システムのデータ同士を直接解析することはできない。そこで、ASTER
GDEM データ（Kanto̲DEM.tif）を Landsat 画像と同じ UTM ゾーン54N/WGS84 (EPSG:32654)に座標変換する。
[Raster]-[Projection]-[Warp(reproject)]を選択肢、入力ファイルとして- 18 -

衛星・GIS データ処理の基礎

Kanto̲DEM.tif を選択する。また出力ファイル名を”Kanto̲DEM̲UTM.tif”とする。Source SRS には、Kanto̲DEM.tif の空間参照システムである
EPSG:4326 が表示される。Target SRS については Target SRS の右側のSelect ボタンをクリックし、UTM ゾーン 54N/WGS84 (EPSG:32654)を選択する。OK をクリックすると、投影変換された Kanto̲DEM̲UTM がレイヤパネ
ルに表示される。

４．3

地形解析

ツールバーから[Raster]-[Analysis]-[Aspect]（傾斜方位画像）を選択。Input file- 19 -

衛星・GIS データ処理の基礎

で Kanto̲DEM̲UTM.tif を選択する。地形解析ツールには、Aspect のほかにも陰影画像(Hillshade)、傾斜角画像(Slope)がある。それぞれの画像を作成する。
・陰影画像、傾斜角画像、傾斜方位画像を作成し表示する。

陰影画像の表示結果- 20 -

衛星・GIS データ処理の基礎

５．

ベクタデータを用いた解析

実際の衛星データ解析では、植生の分布や量を河川、標高などと対応させて調べたり、行政区分や動植物の生息域ごとに解析したりする場合がある。ここでは、このような解析の一例として、ベクタデータを活用して特定の領域や河川周辺の植生指数を調べてみる。QGIS ではベクタデータは Shape ファイルとして保存される。
５．１

新規ベクタレイヤの作成①

新規ベクタレイヤの追加[Layer]-[Create Layer]-[New shapefile layer…]を選択。下のダイアログボッ
クスの様に、ベクタの Geometry type を Polygon に設定、空間参照系はLANDSAT 画像に合わせて EPSG 32654 WGS 84/UTM zone 54N とする。
今回は既存の Attributes list「id」を使うため、New attribute は空欄のままで良い。OK をクリックし、保存先ファイル名を入力する。これで layer panel に新しいベクタレイヤが追加される。

③

データの追加- 21 -

衛星・GIS データ処理の基礎

layer panel で追加したベクタレイヤを右クリックしポップアップダイアログから[Toggle Editing]を選択。これでベクタレイヤが編集可能な状態になる。次にツールバーから[edit]-[add polygon feature]で画像上にポリゴンを描画。下の図は３つのポリゴンを追加した状態。ポリゴン描画が終わったら、再度ベクタレイヤを右クリックしポップアップダイアログから[Toggle Editing]を選択し、編集を終了する。

- 22 -

衛星・GIS データ処理の基礎

５．２

ポリゴン内の平均 NDVI と平均標高の計算

Processing toolbox の検索フィールドに Zonal statistics とタイプして Zonalstatistics を検索する。Zonal statistics をクリックしたら Raster レイヤにNDVI を選択し、ポリゴンデータとして作成したデータ名を指定する。Output
column prefix に NDVI と入力し OK を押すと、各ポリゴン内の平均値が算出される。
結果を確認するためには、ベクタレイヤを右クリックしポップアップダイアロ
グから[Open attribute table]を選択する。すると各ポリゴンの統計値（ポリゴン内のピクセル数、合計値、平均値など）が表示される。同様に標高や他のデータも Zonal statistics を使用することで、統計値を算出しそれぞれのポリゴン領域間の比較を行うことができる。
・ NDVI の空間分布と地形の影響を比較する- 23 -

衛星・GIS データ処理の基礎

５．3

河川周辺の植生解析（バッファ解析）

河川周辺の植生解析を行う。今回利用する Landsat-8 画像には利根川が含まれている。この利根川をマウスでなぞってベクタデータの要素の一つであるラインデータを作成する。
①

河川データの作成

まず、５．１節「新規ベクタレイヤの追加」と同じ要領で、新規ベクタレイヤを作成する。その際に、ベクタの type を Line に設定する。
③

河川データの表示

次に５．１節「ベクタデータの追加」と同じ要領で河川のラインデータを追加する。まず、layer panel で追加したベクタレイヤを右クリックしポップアップ
ダイアログから[Toggle Editing]を選択。これでベクタレイヤが編集可能な状態になる。次にツールバーから[edit]-[add line feature]で画像上にラインを描画。
河川データを作成する際に、カラー画像で河川が判別しにくい時には背景画像として近赤外画像を利用すると川の流れが判別しやすい。これは、近赤外では植生の反射は高く、水面の反射が低い特性があるためである。下の画像は、河川のラインデータを作成した様子である。

- 24 -

衛星・GIS データ処理の基礎

ポリゴン描画が終わったら、再度ベクタレイヤを右クリックしポップアップダ
イアログから[Toggle Editing]を選択し、編集を終了する。
③

バッファ解析

作成した河川データから「河川の周辺２０００ｍの領域」を抽出する。
ツールバーから[Vector]-[Geoprocessing Tools]-[Fixed Distance Buffer]を選択。Input vector layer として作成した River を指定し、Distance を 2000[m]とする。Buffer の部分に出力ファイルの名前を入力し（ここでは Buffer とした）OK を選択すると河川から 20００ｍの領域がベクタデータとして作成される。

- 25 -

衛星・GIS データ処理の基礎- 26 -

衛星・GIS データ処理の基礎

・ ５．２節のポリゴン内の平均 NDVI 計算と同様のやり方で、Zonal statisticsを使って河川周辺２０００ｍの平均 NDVI を計算する。
なお、ベクタデータには、国や県市町村境界ポリゴン、道路網、植生分布、など実に様々なものがあり、政府や自治体などからデータが公開されている。
今回の演習では、自前のベクタデータ（ポリゴン・ライン）を作成したが、より広域な衛星データ解析を行う際には、公開データをうまく活用すると良い。
日本国内の GIS データの例
・ 国土地理院基盤地図情報ダウンロードサービス
Ø

http://fgd.gsi.go.jp/download/

・ 国土交通省国土政策局 国土数値情報Øhttp://nlftp.mlit.go.jp/ksj/

・ 環境省自然環境局
Ø生物多様性センター

http://www.biodic.go.jp/trialSystem/top.html- 27 -

衛星・GIS データ処理の基礎

以上- 28 -