GDALの使い方

この章で学習すること

この章では、代表的な地球観測衛星「LANDSAT-8」が取得した観測データを用いて、GDALライブラリによる衛星データの加工方法を学習します。

  • gdal.Openmatplotlibを用いた画像の表示方法
  • gdal.Translateを用いた衛星画像の切り出し方法、カラー合成方法
  • gdal.Warpを用いた衛星画像の座標系変換方法
  • gdal_pansharpenを用いたパンシャープン画像の作成方法
  • gdal.Translateを用いたフォーマット変換方法

色合いを補正して好みの場所を切り出し、PNGフォーマットへの変換を行うことで企画書の資料に衛星データを貼り付けることも可能です。

企画書の資料に衛星画像を貼り付けて、プレゼン力を高めましょう!!

背景

画像解析に利用するライブラリの多くは、写真や絵などの画像を対象としていますが、衛星画像のような地理空間情報を持つ画像の解析には、専用のライブラリを用いる必要があります。

GDALとは、Geospatial Data Abstraction Libraryの略で、衛星データのような 地理空間情報を含む画像 を扱う際に、地理空間情報を考慮した加工が可能なライブラリのことです。

使用するデータ

衛星概要

この章では、LANDSAT-8衛星のデータを利用します。

LANDSATシリーズは、NASAが打ち上げた陸域観測衛星で、1972年に打ち上げられたLANDSAT-1から、およそ40年間にわたり、継続して地球の状態が観測されています。

継続された観測と並んで、もう一つ大きな特徴として、データの入手が無料です。 費用をかけず、解析が行えるため、スタートアップに適しているデータといえます。

現在は、LANDAST-7(1999年4月打ち上げ)とLANDAST-8(2013年2月打ち上げ)が運用されています。

LANDAST-8に搭載されているセンサー「Operational Land Imager (OLI)」は、以下の通り9つのバンドを持っています。

[バンド番号]: [波長域]  [特徴]

  • バンド1:0.43-0.45μm 沿岸/エアロゾル
  • バンド2:0.45-0.51μm 青色
  • バンド3:0.53-0.60μm 緑色
  • バンド4:0.63-0.68μm 赤色
  • バンド5:0.85-0.88μm 近赤外域
  • バンド6:1.56-1.66μm 中間赤外1
  • バンド7:2.10-2.30μm 中間赤外2
  • バンド8:0.50-0.68μm パンクロマティック
  • バンド9:1.36-1.39μm 中間赤外3

wavelength.png

私たちの目に見えている色は、光の反射を見ています。光は波の性質を持っており、波の谷から谷までの距離を波長と呼びます。 その中でも人間の目に見えている波長の帯域は、「可視域」と呼ばれています。

可視域の中で短い波長は青く見え、長い波長は赤く見えます。赤色よりさらに長い波長には、赤外域があります。

また、物質に光を当てると、その物質によって波長の反射の強さは異なるという特性があります。

OLIは9つのバンドを持っていますが、それぞれのバンド毎に捉える波長が異なっており、人間の目で見えない波長域も捉えることで、多くの情報を得ています。詳細については、リモート・センシング技術センターのサイトを参照してください。

画像のダウンロード方法


それでは、LANDAST-8-OLIのデータをダウンロードしてみましょう。

今回は、AWS(Amazon Web Service)のサイトから、LANDAST-8の画像を呼び出す方法について紹介します。

AWSに格納されているLANDSAT-8のデータを取得するには、まず、入手したいデータの「パス番号」「フレーム番号」「プロダクト名」を指定する必要があります。

この情報を得るためには、LANDSAT-8のデータ検索サイトであるEarth Explorerで検索する必要があります。

Earth Explorerでは、以下の3要素を指定し、「Results」ボタンを押すと、画像を検索することができます。(図1、図2)

  • 観測時期
  • 観測範囲
  • 衛星データ

図1は、サンプルとして、「2020年3月1日~3月31日」の期間の「沖縄」の領域を指定することにします。 左下の赤矩形部分に、期間を入力します。 次に、画面中央のMAP上で、検索したい領域範囲をクリックすることで指定できます。

期間と領域が指定出来たら、左下のData Setsをクリックします。

In [ ]:

usgs_A.png

図1. USGS EarthExoplorer 説明(1/2)

衛星を選択する画面に移るので、今回は、LANDSATのCollection1のLevel-2というデータを選択します。

Resultsボタンをクリックすることで 検索結果が表れます。

usgs_B.png

図2. USGS EarthExoplorer 説明(2/2)

検索結果が得られたら、

1.{プロダクトID}、

2.{Path番号}、

3.{Row番号}

を取得し、以下のURLを指定することで AWSに格納されているLANDSAT-8のデータを取得することができます。

"/vsicurl/http://landsat-pds.s3.amazonaws.com/c1/L8/{パス番号}/{フレーム番号}/{プロダクト名}/{プロダクト名}_{バンド番号}.TIF"

GDALのインストール方法

続いて、今回利用するGDALをインストールしましょう。

ここでは、Window環境に対応するGDALのインストール方法についてご紹介します。

これから、Python環境でGDALを扱う方には、Anacondaがおすすめです。

Anacondaのインストールについては、このサイトで丁寧に説明されているので参考にしてください。

Anacondaのインストールが完了したら、以下の手順でGDALを設定します。

①「Windowsボタン」→「Anaconda Prompt」を右クリックで起動し、管理者権限で実行する。

②ターミナルが出るので、以下のコマンドを実行します。

conda update --all

実行すると、関連パッケージが最新版に更新されます。

③続いて、以下のコマンドを実行します。

conda install -c conda-forge gdal==3.2.0

※3.2.0はバージョン番号であり、2020年12月時点での最新バージョンです。アップデート状況により、適宜変更してください。

これでGDALの設定は終了です。それでは、GDALを使った画像加工処理を始めましょう。

衛星データの表示

まずは、画像を読み込みましょう。 画像を読み込むには「gdal.Open()」を利用します。

AWSのサイトから、LANDSAT-8のバンド5(近赤外)のデータを取得することとします。

今回は、雲の少ない観測データを例として、「2020年4月29日」に観測された関東近辺のデータを取得します。

In [ ]:
#LANDSAT-8 ファイル指定
fpath5=  "/vsicurl/http://landsat-pds.s3.amazonaws.com/c1/L8/107/035/LC08_L1TP_107035_20200429_20200509_01_T1/LC08_L1TP_107035_20200429_20200509_01_T1_B5.TIF"

次に、gdal.Openを利用して、画像の読み込み作業を行います。

pythonでは、以下のように書きます。

{変数名}=gdal.Open({ファイル名})

先ほどの、LANDAST-8のデータを読み込んでみましょう。

In [ ]:
import gdal

#LANDSAT-8 gdal.Openを利用して画像の読み込み作業
band5_image=gdal.Open(fpath5)

gdal.Openを使うことで、{band5_image}という変数に画像データを読み込みました。

読み込んだ画像を表示してみましょう。

ここでは、画像を配列情報[Array]に変換し、表示する方法を紹介します。 画像情報を配列に変換する方法は、以下になります。

{変数名}={画像情報}.ReadAsArray()

また、配列情報を、画像表示するためには、plt.imshowを利用します。

In [ ]:
import numpy as np
import matplotlib.pyplot as plt

#IR_Band_arrayという変数に、近赤外バンドの配列を入れます
NIR_Band_array   = band5_image.ReadAsArray()

#plt.imshowで、配列を画像に表示します。
plt.imshow(NIR_Band_array)
Out[ ]:
<matplotlib.image.AxesImage at 0x7efcecd8b780>

上のような画像が出てきました。

plt.imshowで、そのまま表示すると、暗めの画像が表示されてしまいます。

今回の画像は、ひとつひとつのピクセルの持つ値を「0」から「65,545」の間で持つように定められているのですが、実際は、「30,000」以下の値を持つピクセルが多いです。

何も指定しないと、plt.imshowでは、「0」から「65,545」の尺度で表示しようとするので暗く見えますが、表示する値の最小値と、最大値の指定を行うことで、明るさを変えることができます。

そのためには、それぞれのピクセルがどのような値を持っているのか、「ヒストグラム」を作成して確認します。

ヒストグラムを作成するためには、「plt.hist」を使用します。

In [ ]:
import matplotlib.pyplot as plt
plt.figure(figsize=(20,3))
#先ほど配列に変換したデータの輝度値について、300本の柱で表示することにします
plt.hist(NIR_Band_array.flatten(),bins=300,range=(1, 65545)) 
# 図表の表示
plt.show()

ヒストグラムを確認すると、今回使用する画像は、画素値が6,000から25,000のあたりに集まっているようです。 反対に、30,000以上の値を持つデータはほとんどありません。

上記より、今回の画像では、最小値と最大値を、6,000と25,000に定めたいと思います。

最小値と最大値の指定は、以下の形で指定します。

vmin={最小値}, vmax={最大値}

また、値が低いところは黒、値が高いところは白とすることにしましょう。

色の指定は"cmap"で行います。

cmap={カラーバー指定}

白から黒のカラーは、'gray'とすることで設定できます。

合わせて、以下のように、タイトルと、カラーバーも入れてみましょう。

In [ ]:
#図のサイズを設定します。
plt.figure(figsize=(10,10))

#"NIR_band_array"を、6,000から25,000の範囲で表示します。
plt.imshow(NIR_Band_array,vmin=6000,vmax=25000,cmap='gray')

#タイトルを設定します
plt.title("NIR_Band",fontsize=18)

#カラーバーを設定します
plt.colorbar()

plt.show
Out[ ]:
<function matplotlib.pyplot.show>

LANDSAT-8画像の近赤外バンドを表示することができました。

同様に、青、緑、赤の3つのバンドも、それぞれヒストグラムををそれぞれ表示してみます。

In [ ]:
fpath2 = "/vsicurl/http://landsat-pds.s3.amazonaws.com/c1/L8/107/035/LC08_L1TP_107035_20200429_20200509_01_T1/LC08_L1TP_107035_20200429_20200509_01_T1_B2.TIF"
fpath3 = "/vsicurl/http://landsat-pds.s3.amazonaws.com/c1/L8/107/035/LC08_L1TP_107035_20200429_20200509_01_T1/LC08_L1TP_107035_20200429_20200509_01_T1_B3.TIF"
fpath4 = "/vsicurl/http://landsat-pds.s3.amazonaws.com/c1/L8/107/035/LC08_L1TP_107035_20200429_20200509_01_T1/LC08_L1TP_107035_20200429_20200509_01_T1_B4.TIF"

band2_image=gdal.Open(fpath2)
band3_image=gdal.Open(fpath3)
band4_image=gdal.Open(fpath4)

BlueBand_array  = band2_image.ReadAsArray()
GreenBand_array = band3_image.ReadAsArray()
RedBand_array   = band4_image.ReadAsArray()
In [ ]:
plt.figure(figsize=(20,9))

plt.subplot(3,1,1)
plt.hist(BlueBand_array.flatten(),bins=300,range=(1, 65545)) #青バンドの配列のヒストグラムを表示
plt.title("BlueBand",fontsize=10)

plt.subplot(3,1,2)
plt.hist(GreenBand_array.flatten(),bins=300,range=(1, 65545)) #緑バンドの配列のヒストグラムを表示
plt.title("GreenBand",fontsize=10)

plt.subplot(3,1,3)
plt.hist(RedBand_array.flatten(),bins=300,range=(1, 65545)) #赤バンドの配列のヒストグラムを表示
plt.title("RedBand",fontsize=10)

# 図表の表示

plt.show()

ヒストグラムを確認すると、青バンドは【8,000から15,000】のあたりに値が集まっています。

緑バンドも、少し低くなり、およそ【7,000から13,000】のあたりに値が集まっています。

赤バンドは、【6,000から12,000】あたりで表示したいと思います。

どのバンドの画像かわかりやすくするため、今回は疑似的に『バンドの色』を使って濃淡を表示することにします。

In [ ]:
plt.figure(figsize=(30,30))

#Blueバンドの表示
plt.subplot(1,3,1)
plt.imshow(BlueBand_array,vmin=8000,vmax=15000,cmap='Blues')
plt.title("BlueBand")

#Greenバンドの表示
plt.subplot(1,3,2)
plt.imshow(GreenBand_array,vmin=7000,vmax=13000,cmap='Greens')
plt.title("GreenBand")

#Redバンドの表示
plt.subplot(1,3,3)
plt.imshow(RedBand_array,vmin=6000,vmax=12000,cmap='Reds')
plt.title("RedBand")
plt.show
Out[ ]:
<function matplotlib.pyplot.show>