【実習編】～Landsat8衛星画像で植生を映えさせた地図を描こう～
I. はじめに
衛星画像というと、最近は機械学習や深層学習のためのおもちゃデータと化し、データの入手から分析までの一連の流れも一気に変わったように感じます。
現在の分析環境の主流は2つ・・・
・pythonやRを用いてJupyter notebook環境（Google ColaboratoryやTellusを含む）で処理を回す
・Google Earth Engine（JavaScriptを使ってブラウザ上で操作）を使って処理を回す
でしょうか。高速で進捗を得るために、分析環境が日々進化しています。
その一方で、GISソフトを使った衛星画像の処理は、あっという間に世界の流れから置き去りにされてしまったような感があります。GISは、衛星画像解析の機能を豊富に抱えているので、外部からそれら機能が呼び出されているという点でかろうじて貢献している感じですね。
とはいえ、上記の環境はデータサイエンスの最先端の話。実務や教育の現場ではGISによる衛星画像処理もまだまだ使われています。衛星画像の扱いに慣れるという点では、GISを用いた対話型の操作の方が理解しやすいかもしれません。
今回のnoteの内容は、QGISで衛星画像を扱う方法です。上記の通り、この方法は最先端のトピックではありませんが、このnoteのコンセプトである「非専門家のための」方法として、まずはここから衛星画像の扱いに慣れていきましょう。
II. 衛星画像と空中写真の違い
空から地表を撮影したデータには衛星画像と空中写真（航空写真）があります。この両者の違いを簡単に確認しておきましょう。
衛星画像は、人工衛星に搭載された光学センサーを用いて取得されたデジタルデータです。光学センサーは光をある範囲の波長ごとに感知し、その値をデータ化したものが衛星画像です。つまり、赤色の波長の画像、緑色の波長の画像というように波長ごとの画像があります。また、光学センサーの中には人間の目には見えない近赤外線を捉えるものもあります。人工衛星によって搭載している光学センサーが異なり、得られる画像の波長域も異なります。人工衛星は一定の速度で地球の周りをまわっているので、ある地域を定期的に撮影した時系列データが得られます。
空中写真（航空写真）は、航空機に搭載したカメラを使って地上を撮影した写真です。人工衛星に比べれば低い高度から撮影するため、高解像度の画像が得られます（※1）。一般に、撮影された空中写真はRGBで記録され、衛星画像のように波長ごとに分けることはできません。航空写真は、航空機を飛ばせばすぐに撮影できるので、災害時対応などに向いていますが、コストがかかるので頻繁に撮影することは現実的ではありません。
端的に言えば、衛星画像は分析に適したデータです。一方、空中写真は目視判読に適したデータで、GISでは背景図としても利用されるデータになります。
※1 人工衛星の高度は400km～800km程度。航空機の飛行高度は300m～4000m程度。単位が違います。超低高度衛星と呼ばれた「つばめ」（2019年運用終了）でさえ高度は約200kmなので飛行域が全然異なります。
III. 衛星画像の入手
かつて衛星画像は軍事的懸念から公開が制限されたりわざと解像度を落として公開されたりしていました。しかし、最近は衛星画像の積極的な活用の機運が高まっており、特に解像度が10m以上のものはオープンデータ化が進められています。現在、無料で利用できる衛星画像の有名どころといえば、「Landsat8」と「Sentinel-2」でしょう（※2）。今回はLandsat8を利用してみましょう。
Landsat8は、米国地質調査所（USGS）が運営する人工衛星です。空間分解能（解像度）は30mで高くはないのですが、11バンドの画像を利用でき、そのバンド区分が分かりやすいので最初に利用する衛星画像としてはおすすめです。
※2 Landsat8は「ランドサット8号」と読みます。Sentinel-2は「センチネル2号」です。
日本上空を撮影したものであれば、産業総合研究所（産総研）が運営するLandBrowserというサービスから、衛星画像をダウンロードできます。下記URLからLandBrowserにアクセスしましょう。
LandBrowserを開くと、下図のような地図をベースとしたページが開きます。画面に表示された地域から画像が検索されます。まずは、衛星画像をダウンロードしたい地域に地図を移動させましょう。
画面の右上にあるツールバーやメニューボックスから、衛星画像の検索や保存などの操作ができます。よく使うものをピックアップして紹介します。
Data Selector：
衛星画像の撮影日時を選択します。バーの右端が最新で、左にいくほど過去のデータになります。
Satellite：
検索対象とする人工衛星を選択できます。
Cloud：
衛星画像の中には、地表が雲で覆われ何も見えないものもあります。このバーの数値は雲に覆われている割合を示しており、例えば「50」にすると、雲に覆われている範囲が画像全域の50%以下である衛星画像のみが検索対象となります。
Save：
表示している衛星画像を保存します。ボタンをクリックするとダウンロード設定ウィンドウが別に開きます。このウィンドウは、撮影した人工衛星によって項目が異なります。
産総研ウェブサイトに操作マニュアルがPDFで公開されています。
ここでは、例として新潟県を撮影したLandsat8データをダウンロードします。まずは、欲しい衛星画像を地図上に表示させましょう。次のような手順で操作します。
1. 地図を新潟県が中心に来るように移動する。
2.「Satellite」を「LANDSAT8」にして、検索対象をLandsat8のみにする。
3. 「Data Selector」のバーを左に動かし、2019年5月20日撮影の衛星画像を探す。
地図上に衛星画像が表示されたら、[ Save ] ボタンをクリックしましょう。すると、ダウンロード設定ウィンドウが開きます。
まず、EPSGを4326、画像の形式をTIFFに設定しましょう（※3）。その状態で、band2、band3、band4、band5の4種類をダウンロードします。1つの画像につき、約60MBあります。
ここで、Landsat8の各バンドが何を表すのか確認しておきましょう。前述のように、Landsat8では電磁波を11のバンド帯で測定しています（※4）。そのバンド域は下の表のとおりです。
植物や土、水面、雲などはそれぞれ電磁波に対する固有の反射特性、放射特性を持っています。そこで、バンド帯を組み合わせて反射特性から地物の判別が行えるように、バンドが区分されています（※5）。
今回ダウンロードしたband2、band3、band4、band5はそれぞれ、青色、緑色、赤色、近赤外の波長に対応した画像ということになります。
※3 EPSG：3857はWebメルカトル図法と呼ばれる図法で、GoogleMapsなどウェブ上で地図を表示する際によく利用される図法です。EPSG：4326はWGS84の地理座標系です。EPSG：4326は東西方向に間延びするので、描いた地図をアウトプットする際には投影座標系に変換しましょう。
画像の形式には他にJPEGとPNGが選択できます。しかし、位置情報を持つことができる形式はTIFFのみなので、GISで使うならばTIFF一択です。
※4 Landsat8は2013年に打ち上げられた人工衛星で、それ以前にはLandsat1から7までありました。Landsat1からLandsat8まで、搭載するセンサーは同じではありません。そのため、バンド数や観測するバンド帯が過去のLandsatと現在のLandsat8とは少し異なるので要注意です。
※5 Landsat8はマルチスペクトルセンサーによって、可視光から近赤外までの電磁波を測定しています。このセンサーは広域のスペクトルを区別できますが、代わりに空間分解能（解像度）が低いという欠点があります。Landsat8のバンド8はパンクロマチック（パン、パンクロ）と呼ばれる、特定の1つの波長帯のみを観測する代わりに空間分解能が高い画像になっています。このパン画像とマルチ画像を融合し、マルチ画像でありながら高分解能な画像を作り出す、パンシャープンと呼ばれる技術があります。
IV. QGISでランドサット画像を表示する
それではダウンロードした4つの画像をQGISで読み込んでみましょう。読み込み方法は、データソースマネージャからでも、ドラッグ&ドロップでも構いません。
QGISに表示された画像は、新潟県を描いていますが、どれも白黒です。LandBrowserで表示されていたプレビュー画像と異なりますね。
衛星画像が表すデータは、物体から反射された光の強さです。例えばband2の画像であれば、青色の光の反射強度がどの程度か、band3の画像ならば緑色の光の反射強度がどの程度かを数値で保存しています。現在表示されている画像は、その数値を白黒のグラデーションで描いたものになります（※6）。
一方、実際に見える色は、青色・緑色・赤色それぞれの反射強度が組み合わさって決定されます。したがって、各バンドのファイルが別々になっていると、普段目にする地表のカラー画像を描くことができません。
※6 プロパティでレンダリングタイプを「単バンド疑似カラー」にすれば、値の大小をカラーグラデーションで描くことはできます。しかし、どう頑張っても単バンド画像では空中写真のような自然な色合いを描くことはできません。
そこで、各バンドに分かれた画像を1つにまとめる操作を行います。
メニューバーの「ラスタ」から「その他」→「仮想ラスタの構築」を選択しましょう。
Inputlayersでは、QGISに読み込んだ4つのラスタにチェックを入れます。この時、ファイルの並びを上からバンド2、バンド3と順に並べておきます。並びの入れ替えはドラッグで行います。
保存ファイルは「.vrt」という形式のテキストデータです。.vrtファイルには仮想ラスタを構築元となったラスタファイルのパスなどが保存されています。したがって、.vrtファイルを作成した後にLandsat画像を移動すると、仮想ラスタを正しく読み込めなくなる場合があるので注意しましょう。処理を実行するとレイヤに仮想ラスタが追加されます。
「仮想ラスタの構築」は、複数のラスタデータをレイヤー上で1つのファイルとして振舞うように設定する機能です。ここで重要なのは、仮想レイヤにすることで各ラスタをマルチバンドカラーに割り当てられるようになることです。
実際に操作してみましょう。仮想レイヤのプロパティを開き、シンボロジを選択します。レンダリングタイプを「マルチバンドカラー」にすると、「赤のバンド」、「緑のバンド」、「青のバンド」にそれぞれバンド1～3が割り当てられており、ドロップダウンから割り当てを変更できます。ここでバンド1～4は仮想ラスタを作成したときの並び順なので、バンド2から順に並べていれば、バンド1が青色、バンド2が緑色、バンド3が赤色、バンド4が近赤外に該当しています。
したがって、次のようにバンドを割り当てれば、自然な色合いで衛星画像が描画されるはずです（※7）。
※7 複数のバンドをもつラスタのことをマルチバンドラスタと呼びます。また、各バンドに赤・緑・青の色を割り当て、カラー合成することをRGBコンポジットと呼びます。
「赤のバンド」・・・バンド3（赤の可視光）
「緑のバンド」・・・バンド2（緑の可視光）
「青のバンド」・・・バンド1（青の可視光）
バンドの割り当てを変更した結果が下の画像です。海は青く、山は緑、都市部は灰色っぽくなりました・・・。色合いは良さそうですが、全体的に暗く鮮やかさが足りないですね。
仮想ラスタのプロパティを開き、「ヒストグラム」を表示しましょう。ヒストラムが表示されない場合は、[ ヒストグラムの計算 ] ボタンをクリックします。
ヒストグラムでは、各バンドが取る値の分布を確認できます。Landsat8では、値（反射強度）を16bitで表現しているので、取り得る値は0～65535まであります。
しかし、実際には各バンドとも値は5000～10000に集中しており、その他の値はほとんど無いことが分かります。そのため、全体に暗く色彩に乏しい表現となっています。
そこで、有効な値の幅を変更し、使える色の幅を広げてあげましょう。
ヒストグラムを確認すると、バンド1～3は5000～15000の幅があれば良さそうです。バンド4は5000～30000の幅があれば良いでしょう。
幅の設定は、ヒストグラムの下にある項目か、シンボロジから行うことができます。ヒストグラムを見ながら詳細に設定をする場合はヒストグラムウィンドウから、大雑把な設定でよければシンボロジから設定しましょう（※8）。下図では、バンド1～3の値の有効範囲をシンボロジからバッサリと指定しています。
最小、最大の幅変更を適用したものが下の画像です。全体的に鮮やかになり、GoogleMapsの航空写真などと同じような画像になりました。
※8 取り得る値の幅を狭めることで、変更後のヒストグラムは横に伸びます。このように、値の分布を均す操作のことをヒストグラム平坦化といいます。
V. 植生を強調した画像表現
衛星画像の赤の波長域画像をQGISで赤色に、緑の波長域画像を緑色に、青の波長域画像を青色に割り当てれば自然な色合いの画像が描ける。これは当たり前のことであるとともに、これだけでは何も面白くありません。むしろ自然な色合いの画像が欲しければ、より解像度が高く雲に覆われることもない空中写真の方が適しているでしょう。
衛星画像の優れたところは、波長域と色の割り当てを自由に決められるところにあります。特に、近赤外線という人の目には見えない光の値に色を割り当てられる点が重要です。
次の図は光の波長帯と反射・放射特性を表したものです。植物は、緑色の光をよく反射するので私たちには緑色に見えています。しかし、実は緑色以上に近赤外の光を反射しています。ということは、近赤外の画像を緑色に割り当てれば植生を強調した画像になるのではないでしょうか。
JAXAウェブサイトより：http://www.sapc.jaxa.jp/use/data_view/
実際に設定してみましょう。仮想ラスタのシンボロジで、「緑のバンド」をバンド4（近赤外）に変更します。また、最小・最大は先ほど確認したとおり、5000～30000に設定します。「青のバンド」にはバンド2（緑の可視光）を割り当てます。これで設定を適用すると次のような図になります。
森林部分の緑が輝きを増したような画像になりました。また、都市部でも川に沿って緑が強く出ていることが分かります。
続いて、次のようなバンドの割り当て方をするとどうなるでしょうか。
「赤のバンド」・・・バンド4（近赤外）
「緑のバンド」・・・バンド3（赤の可視光）
「青のバンド」・・・バンド2（緑の可視光）
すると、植生が分布するところが赤色で表現されます。リモートセンシングの分野では、植生を強調した図を描く場合には、このような植生を赤で表す表現方法がよく用いられます。
ここまで、作成した3つの画像には、それぞれ名前が付いています。各画像の名前とその説明は以下の通りです。
トゥルーカラー画像：
可視光の赤・緑・青の波長域と、色の赤・緑・青の割り当てを対応させて作成した画像。普段見るような、自然な配色の画像になります。
ナチュラルカラー画像：
緑に近赤外の波長域を割り当てて作成した画像。植物は近赤外を強く反射するので、植生部分が強調された画像になります。
フォールスカラー画像：
近赤外の波長域を赤に割り当て作成した画像。人は赤色を他の色より強調して認識する性質があります。そこで、近赤外の波長域を赤色に割り当てることで、ナチュラルカラー画像よりも植生部分を強調して見せることができます。
VI. 正規化植生指標（NDVI）を算出
前章では描画に使用するバンドや割り当てる色を工夫することで植生の分布を強調した画像を作成しました。次は、バンド演算という方法を用いて、植生の状況を可視化する方法を紹介します（※9）。
※9 バンド演算とは、マルチバンドラスタのバンド間で計算処理を行うことです。特別な名前が付いていますが、やることはラスタファイル同士で計算処理を行うラスタ演算と同じです。
ここでは正規化植生指標（NDVI）を算出してみましょう。植生指標とは、植生の活性状況を表す指標です。元気な植物は活発に光合成を行うため、光合成に使う赤の波長域の光をより多く吸収しているはずです。そこで、次のような式を考えます。
バンド4 - バンド3
これは近赤外から赤の可視光の値を引き算する計算になります。光合成を活発に行っている植物が多い場所では、赤の波長域の反射強度を表すバンド3の値が小さいはずなので、この計算結果の値は大きくなります。一方、近赤外を強く反射するけど赤の波長域の光も反射するようなものは（※10）、値が小さくなるので活性の高い植物と区別できます。
※10 光合成が十分に行われていない病気の植物など。
実際には、上記の式を次のように修正（正規化）します。こうすると、算出される値は-1～1の値を取るようになり、指標として扱いやすくなります。
(バンド4 - バンド3) / (バンド4 + バンド3)バンド3やバンド4は今作業している仮想ラスタでの表現なので、一般化して表現すれば次のようになります。ここでIRは近赤外、Rは赤の可視光を意味します。
NDVI = ( IR - R ) / ( IR + R )このNDVIの値が1に近いほど植物活性が高いことを意味します。一般にはこの式で求められる値をNDVIと呼びますが、マイナスの値が出ることを嫌う場合は、さらに修正項を追加して値が非負となるようにすることもあります。
それでは実際に計算してみましょう。メニューバーの「ラスタ」から「ラスタ計算機」を選択します。新しく開いたラスタ計算機ウィンドウで、下図のように計算式を入力します。[ OK ] ボタンをクリックすると処理が開始され、処理が終わるとレイヤにNDVIのラスタが追加されます。
NDVIラスタのシンボルを設定しましょう。プロパティのシンボロジでレンダリングタイプを「単バンド疑似カラー」に変更し、グラデーションを選択します。ここではグラデーションに「Viridis」を選択しました。
作成したNDVIを眺めてみましょう。ここで、背景に地理院地図を描画し、NDVIのカラーレンダリングの混合モードを「乗算」にします。すると、NDVIと地図が重なって表示され、植生の活性が高い場所とその土地の情報が分かりやすくなります。
新潟県は米の産地として有名ですが、NDVIを見る限り水田に植物活性がみられません。新潟県の田植えは5月下旬から6月上旬に行われるのが一般的のようです。この衛星画像が撮影された5月20日は、まだ田植えが行われてなさそうです。
一方、畑や果樹園の地図記号が見られる場所では、植物活性が高い傾向にあります。また、河川敷や荒れ地も植物活性が高くなっています。雑草繁茂の勢いが見て取れますね。このように、土地利用の状況を植物活性から確認できます。
VII. おわりに
今回はLandsat8が撮影した衛星画像を入手する方法と、その衛星画像を用いて植生を強調した画像の作成方法を紹介しました。バンドの割り当てやバンド演算を使って、植生の特徴引き出す過程は機械学習の基本的な操作とも共通する部分があります。
Landsat8は解像度が30m程度ですが、Sentinel-2やALOSのようにLandsat8より高解像度でありながら無料で利用できるデータも続々と公開されています。それぞれ搭載しているセンサーも異なり、Landsat8では捉えられない波長のデータを持っているものもあります。
衛星画像の高解像度化やセンサーの種類の増加に伴って、衛星画像の利用方法も物体検知や山火事の被害状況把握、サンゴ礁海域の水深・海底地形把握など新しい応用が日々登場しています。
他の衛星データを使った、他の衛星画像の利用事例も今後のnoteで紹介したいと思います。
記事の内容でもそれ以外でも、地理やGISに関して疑問な点があれば、可能な限りお答えしたいと思います。