GDALの使い方¶
この章で学習すること¶
この章では、代表的な地球観測衛星「LANDSAT-8」が取得した観測データを用いて、GDALライブラリによる衛星データの加工方法を学習します。
gdal.Openとmatplotlibを用いた画像の表示方法gdal.Translateを用いた衛星画像の切り出し方法、カラー合成方法gdal.Warpを用いた衛星画像の座標系変換方法gdal_pansharpenを用いたパンシャープン画像の作成方法gdal.Translateを用いたフォーマット変換方法
色合いを補正して好みの場所を切り出し、PNGフォーマットへの変換を行うことで企画書の資料に衛星データを貼り付けることも可能です。
企画書の資料に衛星画像を貼り付けて、プレゼン力を高めましょう！！
背景¶
画像解析に利用するライブラリの多くは、写真や絵などの画像を対象としていますが、衛星画像のような地理空間情報を持つ画像の解析には、専用のライブラリを用いる必要があります。
GDALとは、Geospatial Data Abstraction Libraryの略で、衛星データのような 地理空間情報を含む画像 を扱う際に、地理空間情報を考慮した加工が可能なライブラリのことです。
現在は、LANDAST-7（1999年4月打ち上げ）とLANDAST-8（2013年2月打ち上げ）が運用されています。
LANDAST-8に搭載されているセンサー「Operational Land Imager (OLI)」は、以下の通り９つのバンドを持っています。
[バンド番号]： [波長域] [特徴]
- バンド１：0.43-0.45μm 沿岸/エアロゾル
- バンド２：0.45-0.51μm 青色
- バンド３：0.53-0.60μm 緑色
- バンド４：0.63-0.68μm 赤色
- バンド５：0.85-0.88μm 近赤外域
- バンド６：1.56-1.66μm 中間赤外１
- バンド７：2.10-2.30μm 中間赤外２
- バンド８：0.50-0.68μm パンクロマティック
- バンド９：1.36-1.39μm 中間赤外３
私たちの目に見えている色は、光の反射を見ています。光は波の性質を持っており、波の谷から谷までの距離を波長と呼びます。 その中でも人間の目に見えている波長の帯域は、「可視域」と呼ばれています。
可視域の中で短い波長は青く見え、長い波長は赤く見えます。赤色よりさらに長い波長には、赤外域があります。
また、物質に光を当てると、その物質によって波長の反射の強さは異なるという特性があります。
OLIは９つのバンドを持っていますが、それぞれのバンド毎に捉える波長が異なっており、人間の目で見えない波長域も捉えることで、多くの情報を得ています。詳細については、リモート・センシング技術センターのサイトを参照してください。
画像のダウンロード方法
それでは、LANDAST-8-OLIのデータをダウンロードしてみましょう。
今回は、AWS（Amazon Web Service)のサイトから、LANDAST-8の画像を呼び出す方法について紹介します。
AWSに格納されているLANDSAT-8のデータを取得するには、まず、入手したいデータの「パス番号」「フレーム番号」「プロダクト名」を指定する必要があります。
この情報を得るためには、LANDSAT-8のデータ検索サイトであるEarth Explorerで検索する必要があります。
Earth Explorerでは、以下の3要素を指定し、「Results」ボタンを押すと、画像を検索することができます。（図１、図２）
- 観測時期
- 観測範囲
- 衛星データ
図１は、サンプルとして、「2020年3月1日～3月31日」の期間の「沖縄」の領域を指定することにします。 左下の赤矩形部分に、期間を入力します。 次に、画面中央のMAP上で、検索したい領域範囲をクリックすることで指定できます。
期間と領域が指定出来たら、左下のData Setsをクリックします。
図１． USGS EarthExoplorer 説明(1/2)
衛星を選択する画面に移るので、今回は、LANDSATのCollection1のLevel-2というデータを選択します。
Resultsボタンをクリックすることで 検索結果が表れます。
図２． USGS EarthExoplorer 説明(2/2)
検索結果が得られたら、
１．{プロダクトID}、
２．{Path番号}、
３．{Row番号}
を取得し、以下のURLを指定することで AWSに格納されているLANDSAT-8のデータを取得することができます。
"/vsicurl/http://landsat-pds.s3.amazonaws.com/c1/L8/{パス番号}/{フレーム番号}/{プロダクト名}/{プロダクト名}_{バンド番号}.TIF"
GDALのインストール方法¶
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のバンド５（近赤外）のデータを取得することとします。
今回は、雲の少ない観測データを例として、「2020年4月29日」に観測された関東近辺のデータを取得します。
#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のデータを読み込んでみましょう。
import gdal
#LANDSAT-8 gdal.Openを利用して画像の読み込み作業
band5_image=gdal.Open(fpath5)
gdal.Openを使うことで、{band5_image}という変数に画像データを読み込みました。
読み込んだ画像を表示してみましょう。
ここでは、画像を配列情報[Array]に変換し、表示する方法を紹介します。 画像情報を配列に変換する方法は、以下になります。
{変数名}={画像情報}.ReadAsArray()また、配列情報を、画像表示するためには、plt.imshowを利用します。
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)
上のような画像が出てきました。
plt.imshowで、そのまま表示すると、暗めの画像が表示されてしまいます。
今回の画像は、ひとつひとつのピクセルの持つ値を「0」から「65,545」の間で持つように定められているのですが、実際は、「30,000」以下の値を持つピクセルが多いです。
何も指定しないと、plt.imshowでは、「0」から「65,545」の尺度で表示しようとするので暗く見えますが、表示する値の最小値と、最大値の指定を行うことで、明るさを変えることができます。
そのためには、それぞれのピクセルがどのような値を持っているのか、「ヒストグラム」を作成して確認します。
ヒストグラムを作成するためには、「plt.hist」を使用します。
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'とすることで設定できます。
合わせて、以下のように、タイトルと、カラーバーも入れてみましょう。
#図のサイズを設定します。
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
LANDSAT-8画像の近赤外バンドを表示することができました。
同様に、青、緑、赤の３つのバンドも、それぞれヒストグラムををそれぞれ表示してみます。
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()
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】あたりで表示したいと思います。
どのバンドの画像かわかりやすくするため、今回は疑似的に『バンドの色』を使って濃淡を表示することにします。
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
LANDASAT-8の画像表示の方法を学びました。
続いて、画像の切り出しについて学びましょう。
画像の切り出し¶
衛星画像を用いた解析処理において、衛星データによっては、観測範囲が広いことから、後段の処理が重くなることを避けるために、関心領域のみを切り出して作業を行うことがあります。
「gdal.Translate」を用いることで、画像を切り出すことができます。
今回、切り出しを行う方法について、以下の2通りの方法を示します。
【方法① 切り出し開始位置と幅を設定して切り出す方法】
【方法② 緯度/経度を指定して切り出す方法】
方法① 開始位置と幅を指定した切り出し¶
全体の画像から、『X方向開始位置』、『Y方向開始位置』、『X方向の切り出し幅』、そして『Y方向の切り出し幅』を設定することで、切り出すことができます。
gdal.Translateの"srcWin"オプションを利用します。
gdal.Translate(｛出力画像名｝,{入力画像名},srcWin=[minX,minY,deltaX,deltaY])
画像全体から、東京の都市部を切り出してみようと思います。
全体の画像を見ながら、例えば、以下のように切り出し位置を設定しました。
X方向開始位置：2700
Y方向開始位置：5100
X方向の切り出し幅：300
Y方向の切り出し幅：300
#全体画像から、東京都の切り出し位置を設定
minX=2700
minY=5100
deltaX=300
deltaY=300
#出力画像の名前を指定
cut_NIR_path="TOKYO_NIR_.tif"
#gdal.Translate(｛出力画像名｝,{入力画像名},srcWin=[minX,minY,deltaX,deltaY])
ds=gdal.Translate(cut_NIR_path, fpath5, srcWin=[minX,minY,deltaX,deltaY])
ds=None
切り出し処理が完了したかどうか、画像表示で確認しましょう。
#gdal.Openで画像を読み込みます
cut_NIR_img=gdal.Open(cut_NIR_path)
#画像情報を配列に変換します
cut_NIR_array = cut_NIR_img.ReadAsArray()
#表示パラメータを決めて、画像を表示します
plt.figure(figsize=(10,10))
plt.imshow(cut_NIR_array,vmin=6000,vmax=25000,cmap='Oranges')
plt.title("NIRBand(TOKYO)")
plt.colorbar()
plt.show
皇居付近の”近赤外域のバンド”を切り出して表示することができました。
近赤外域の波長は、植物が持つクロロフィルの活性度と比例関係があるため、植生が多いところで、値は高くなります。 皇居の周りのビル群では、近赤外の値は低く、皇居の内側は近赤外の値が高いことが確認できます。
方法② 4隅の緯度経度を指定した切り出し¶
続いて、4隅の緯度経度を指定して、切り出す方法を示します。
方法①で学習した、『開始位置と幅を指定した切り出し』は、単画像の切り出しの際には有効ですが、複数枚の画像で同じ場所を切り出したい場合は、画像ごとに開始位置などの条件が異なってくるため、それぞれ開始位置と幅を計算する必要があります。
条件が異なる複数枚の画像から、同じ場所を切り出すためには、『緯度経度を指定した切り出し』が有効です。
そこで、解像度の異なる、近赤外バンド（30m解像度）とパンクロマティックバンド（15m解像度）を用いて、同じ場所を切り出してみましょう。
※緯度経度を指定して切り出す方法は、”災害分野”の演習の中でも同じ方法を用いていますので、そちらも参考にしてください。
まず、使用するデータの座標系を確認する必要があります。
データの座標系情報を知るためには、GetProjection()を利用します。
{変数}＝{画像情報}.GetProjection()band5_image=gdal.Open(fpath5)
src_prj=band5_image.GetProjection()
print ("座標系などの情報：")
print (src_prj)
出力された項目を確認することで、今回使用するデータの持つ情報は、測地系が"WGS84"、座標系は、"UTM zone 54N(EPSG:32654)"であることが確認できました。
UTM図法では、位置情報を扱う際の単位はメートルです。
緯度/経度を用いて切り出しを行うためには、緯度経度情報をもつ座標系に変換する必要があります。
画像の持つ座標系を変換するためには、「gdal.Warp」を使います。
今回、EPSG：32654から、EPSG：4326に変換することで、緯度経度情報を持つデータとして扱うことが可能です。
※EPSGコードの詳細については、空間情報クラブのサイトを参照してください。
latlon_band5_path="latlon_band5_img.tif"
ds=gdal.Warp(latlon_band5_path,fpath5,srcSRS="EPSG:32654",dstSRS="EPSG:4326")
ds=None
座標系が変換されているかどうか、もう一度GetProjectionで確認します。
latlon_band5_img=gdal.Open(latlon_band5_path)
src_prj=latlon_band5_img.GetProjection()
print ("座標系などの情報：")
print (src_prj)
座標系が"EPSG:4326"に、変更されていることが確認できました。
緯度経度が参照できる座標系に変換することができたので、緯度経度の情報を使って、画像の切り出しを行います。
事前に、GoogleEarthなどから、関心領域の範囲の4隅の緯度経度を確認しておきます。
今回は、国立競技場近辺を切り出してみましょう。以下の範囲で切り出します。
西の端（minX)：139.69 (deg)
南の端（minY)：35.66 (deg)
東の端（maxX)：139.74 (deg)
北の端（maxY)：35.70 (deg)
緯度経度を用いた切り出しは、以下のように書きます。
方法①では、「srcWin」を用いましたが、今回は、「projWin」を使います。
gdal.Translate(｛出力画像名｝,{入力画像名},prjWin=[minX,maxY,maxX,minY])赤色バンド(バンド４：30m解像度)とパンクロマティックバンド(バンド８：15m解像度)について、同じ緯度経度で切り出してみましょう。
まずは、２つのバンドに対して、それぞれ座標の変換を行います。
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"
fpath8 = "/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_B8.TIF"
#緯度経度の情報を持ったバンド8のファイル名を定義
latlon_band4_path="latlon_band4_img.tif"
latlon_band8_path="latlon_band8_img.tif"
#gdal.Warpで座標系を変換
ds=gdal.Warp(latlon_band4_path,fpath4,srcSRS="EPSG:32654",dstSRS="EPSG:4326")
ds=gdal.Warp(latlon_band8_path,fpath8,srcSRS="EPSG:32654",dstSRS="EPSG:4326")
ds=None
続けて、変換を行った画像に対して、gdal.Warpで切り出しを行います。
minX=139.69
minY=35.66
maxX=139.74
maxY=35.70
#切り出し後のファイルのファイル名を定義
cut_latlon_band4_path="cut_latlon_band4_img.tif"
cut_latlon_band8_path="cut_latlon_band8_img.tif"
#gdal.Translate で切り出しを実行
ds=gdal.Translate(cut_latlon_band4_path,latlon_band4_path, projWin=[minX,maxY,maxX,minY])
ds=gdal.Translate(cut_latlon_band8_path,latlon_band8_path, projWin=[minX,maxY,maxX,minY])
ds=None
これで、画像の切り出しが完了しました。
切り出された範囲を確認するため表示してみましょう。
#それぞれgdal.Openで画像を読み込み
cut_latlon_band4_img=gdal.Open(cut_latlon_band4_path)
cut_latlon_band8_img=gdal.Open(cut_latlon_band8_path)
#それぞれ配列情報を読み込み
cutRedBand_array = cut_latlon_band4_img.ReadAsArray()
cutPanBand_array = cut_latlon_band8_img.ReadAsArray()
#画像を表示
plt.figure(figsize=(20,20))
plt.subplot(1,2,1)
plt.imshow(cutRedBand_array,vmin=6000,vmax=12000,cmap='gray')
plt.title("RedBand")
plt.subplot(1,2,2)
plt.imshow(cutPanBand_array,vmin=5000,vmax=15000,cmap='gray')
plt.title("PanchroBand")
plt.show
緯度経度の情報を用いて、国立新競技場の周辺を切り出すことができました。 赤バンドより、パンクロマティックバンドの方が、より細かく見えていることが確認できます。
続いて、衛星データを用いたカラー合成について学んでいきましょう。
コラム１ gdal.Translateのオプション紹介¶
この演習で利用している"gdal.Translate"には、今回の事例以外にも様々なオプションがあります。その中から一部、よく使うオプションを紹介します。
Pythonにおける gdal.Translate の使い方
gdal.Translate ({出力ファイルパス},{入力ファイルパス},{オプション１},{オプション２},・・・)| オプション名 | 説明 | 使用例 |
|---|---|---|
| format | 出力ファイルの形式を指定する（"GTiff", "JPEG",etc….） | format="Gtiff" |
| outputType | 出力ファイルのデータ型を指定する（gdal.GDT_Byte, gdal.GDT_Int16,etc….) | outputType="gdal.GDT_Byte" |
| width | 出力ファイルのX軸方向のピクセル数を指定する | width=10000 |
| height | 出力ファイルのY軸方向のピクセル数を指定する | height=10000 |
| widthPct | 出力ファイルのX軸方向の縮尺率を指定する（例：200%） | widthPct=200 |
| heightPct | 出力ファイルのY軸方向の縮尺率を指定する（例：200%） | heighPct=200 |
| xRes | 出力ファイルのピクセルサイズ（X方向）を指定する | xRes=1.5 |
| yRes | 出力ファイルのピクセルサイズ（Y方向）を指定する | yRes=1.5 |
| srcWin | 開始位置と幅を指定することで領域を切り出す（left_x,top_y,width,height) | srcWin=(300,450,100,100) |
| projWin | 左上の位置情報と、右下の位置情報を指定することで領域を切り出す | projWin=(36.035,140.058,36.025,140.068) |
| resampleAlg | 出力ファイル生成時の、リサンプリング方法を指定する（例：NN、CC、etc…） | resampleAlg="NN" |
| noData | 出力ファイルのNoDataの値を指定する | noData=-9999 |
| scaleParams | 出力ファイルの表示範囲を指定する。[DN最小値、DN最大値] | scaleParams=[[3000,16000]] |
カラー合成¶
LANDASTデータは、バンドごとに分かれているので、色々なバンドを組み合わせてカラー合成画像を作成してみましょう。
ここでは、「トゥルーカラー（True Color）」「フォールスカラー（False Color）」「ナチュラルカラー（Natural Color）」の画像を作成します。
トゥルーカラー画像とは、人間の目で見える色と同じ配合で作成される画像のことです。
衛星画像の赤の波長を赤バンドに、緑波長を緑バンドに、青波長を青バンドに割り当てた画像になります。
フォールスカラー画像とは、植生部分を赤く表示させた画像のことです。
人間の目は「赤色が強調されて見える」といった特性を利用して、植生の様子を詳細に把握することができます。
近赤外の波長を赤バンドに、赤の波長は緑バンドに、緑の波長を青バンドに割り当てた画像をフォールスカラーと呼んでいます。
ナチュラルカラー画像とは、植生を自然な色合いで強調した画像のことです。
赤の波長を赤バンドに、近赤外の波長を緑バンドに、緑の波長を青バンドに割り当てます。
それでは、カラー画像の作成方法を以下に示します。
まず最初に、LANDSAT-8のデータは、16biのデータ形式ですが、カラー表示を行うため、8bitの形式に変換します。
変換には、gdal.Translateを用います。
gdal.Translate({出力画像名},{入力画像名},outputType={データ型},scaleParams=[[min,max]])さらに、Section2で行った、切り出し作業も、同じgdal.Translateを用いているので、つなげて書くことが可能です。
gdal.Translate({出力画像名},{入力画像名},outputType={データ型},scaleParams=[[min,max]],srcWin=[minX,minY,deltaX,deltaY])LANDAST-8のバンド２、３、４、５を読み込んで、それぞれを8bit画像に変換し、かつ都心部に切り出すことを以下のスクリプトで実行します。
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"
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"
#8bit出力データ名を指定しておく
band5_8bit_path="Band5_8bit.tif"
band4_8bit_path="Band4_8bit.tif"
band3_8bit_path="Band3_8bit.tif"
band2_8bit_path="Band2_8bit.tif"
#切り出しの詳細
minX=2400
minY=4500
deltaX=1500
deltaY=1500
#各バンドのファイルを、それぞれ、関心領域のみ切り出す。出力は8bitのgeotifとする
#gdal.Translate({出力画像名}, {入力画像名}, outputType={データ形式設定} , scaleParams=[[min,max]])
#
gdal.Translate(band2_8bit_path,fpath2,outputType=gdal.GDT_Byte,scaleParams=[[9700,13400]],srcWin=[minX,minY,deltaX,deltaY] )
gdal.Translate(band3_8bit_path,fpath3,outputType=gdal.GDT_Byte,scaleParams=[[8100,12800]], srcWin=[minX,minY,deltaX,deltaY] )
gdal.Translate(band4_8bit_path,fpath4,outputType=gdal.GDT_Byte,scaleParams=[[7000,13000]],srcWin=[minX,minY,deltaX,deltaY] )
gdal.Translate(band5_8bit_path,fpath5,outputType=gdal.GDT_Byte,scaleParams=[[6000,18300]],srcWin=[minX,minY,deltaX,deltaY] )
#作成した8bitの切り出し画像を読み込む
b2_image=gdal.Open(band2_8bit_path)
b3_image=gdal.Open(band3_8bit_path)
b4_image=gdal.Open(band4_8bit_path)
b5_image=gdal.Open(band5_8bit_path)
#読み込んだ画像を配列に変換する
BlueBand_array = b2_image.ReadAsArray()
GreenBand_array = b3_image.ReadAsArray()
RedBand_array = b4_image.ReadAsArray()
NIRBand_array = b5_image.ReadAsArray()
続けて、出力画像の設定、および書き込みを行います。
#出力ファイルの設定のために、入力ファイルのX方向のピクセル数、Y方向のピクセル数を読み出す
Xsize=b2_image.RasterXSize #band2の画像のX方向ピクセル数
Ysize=b2_image.RasterYSize #band2の画像のY方向ピクセル数
dtype=gdal.GDT_Byte
band=3
#出力ファイルの設定を行う(True Color)
out_True_path ="TrueColor_TOKYO.tif" #出力ファイル名
#空の出力ファイルを作成する
out1= gdal.GetDriverByName('GTiff').Create(out_True_path, Xsize, Ysize, band, dtype)#({出力ファイル名}, {X方向のピクセル数},{Y方向のピクセル数},{バンド数},{データ形式})
#出力ファイルの座標系を設定する
out1.SetProjection(b2_image.GetProjection()) #{出力変数}.SetProjection(座標系情報)
out1.SetGeoTransform(b2_image.GetGeoTransform()) #{出力変数}.SetGeoTransform(座標に関する６つの数字)
#Red、Green、Blueバンドの配列を、WriteArrayを用いて出力ファイルの3バンドに書き込む
out1.GetRasterBand(1).WriteArray(RedBand_array) #赤の配列を赤バンドに書き込む
out1.GetRasterBand(2).WriteArray(GreenBand_array) #緑の配列を緑バンドに書き込む
out1.GetRasterBand(3).WriteArray(BlueBand_array) #青の配列を青バンドに書き込む
out1.FlushCache()
上記で作成した、”TrueColor_TOKYO.tif”を、表示してみましょう。
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
plt.figure(figsize=(10,10))
image1 = mpimg.imread(out_True_path)
plt.imshow(image1)
plt.title("True Color TOKYO")
plt.show
綺麗な都心部の画像を得ることができました。
続けて、フォールスカラー、ナチュラルカラーの画像も作成してみましょう。
作成方法はトゥルーカラーと同じです。最後の、配列の書き込み部分が異なってきます。
####出力ファイルの設定を行う(False Color)
#出力ファイルの設定
out_False_path="FalseColor_TOKYO.tif"
out2= gdal.GetDriverByName('GTiff').Create(out_False_path, Xsize, Ysize, band, dtype)
out2.SetProjection(b2_image.GetProjection())
out2.SetGeoTransform(b2_image.GetGeoTransform())
#Red、Green、Blueバンドの配列を、出力ファイルの3バンドに書き出す
out2.GetRasterBand(1).WriteArray(NIRBand_array) #近赤外の配列を赤バンドに書き込む
out2.GetRasterBand(2).WriteArray(RedBand_array) # 赤 の配列を緑バンドに書き込む
out2.GetRasterBand(3).WriteArray(GreenBand_array)# 緑 の配列を青バンドに書き込む
out2.FlushCache()
####出力ファイルの設定を行う(Natural Color)
#出力ファイルの設定
out_Natural_path="NaturalColor_TOKYO.tif"
out3= gdal.GetDriverByName('GTiff').Create(out_Natural_path, Xsize, Ysize, band, dtype)
out3.SetProjection(b2_image.GetProjection())
out3.SetGeoTransform(b2_image.GetGeoTransform())
#Red、Green、Blueバンドの配列を、出力ファイルの3バンドに書き出す
out3.GetRasterBand(1).WriteArray(RedBand_array) # 赤 の配列を赤バンドに書き込む
out3.GetRasterBand(2).WriteArray(NIRBand_array) #近赤外の配列を緑バンドに書き込む
out3.GetRasterBand(3).WriteArray(GreenBand_array)# 緑 の配列を青バンドに書き込む
out3.FlushCache()
3つの画像を、並べて表示してみましょう。
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
plt.figure(figsize=(20,20))
plt.subplot(1,3,1)
image1 = mpimg.imread(out_True_path)
plt.imshow(image1)
plt.title("True Color TOKYO")
plt.subplot(1,3,2)
image2 = mpimg.imread(out_False_path)
plt.imshow(image2)
plt.title("False Color TOKYO")
plt.subplot(1,3,3)
image3 = mpimg.imread(out_Natural_path)
plt.imshow(image3)
plt.title("Natural Color TOKYO")
plt.show
東京湾近辺のカラー画像を確認することができました。 フォールスカラーや、ナチュラルカラーの画像を確認すると、都心部の中でも、皇居や新宿御苑等、まとまった緑があることが良くわかります。
パンシャープン画像の作成¶
LANDAST-8のOLIが持つバンド８（パンクロマティック）は高精度（解像度15m）のデータを取得することができますが、モノクロの画像です。
一方、その他のバンドは、可視光から近赤外線まで複数の波長帯の電磁波を記録しているため多くのスペクトル情報を持っています（＝マルチスペクトル）が、低解像度（解像度30m）です。
パンクロマティック画像（解像度15m）と、マルチスペクトル画像（解像度30m）を用いて、パンシャープン処理を行うことで、解像度15mのカラー画像を作成することができます。
GDALでは、gdal_pansharpenという機能で、パンシャープン画像を作成することができるので、ここではパンシャープン画像作成について紹介します。
gdal_pansharpen.pyに関する詳細は、以下のサイトをご確認ください。
gdal_pansharpen.pyは、Python内のGDALライブラリの標準では無いため、Python内で使用するためには、機能を呼び出す必要があります。
まずは、自身の環境から、gdal_pansharpenがどこにあるか以下のように確認します。
!which gdal_pansharpen.py
"/usr/bin"の下にあることが分かったので、以下のようにimportを行います。確認ができたら、次のようにgdal_pansharpenをインポートします。
import sys
sys.path.append('/usr/bin/')
from gdal_pansharpen import gdal_pansharpen
これで、gdal_pansharpenが使えるようになったので、画像をパンシャープン画像にしてみましょう。
先ほど作成した、カラー合成と同じ範囲のパンクロマティックデータを切り出してパンシャープンを行います。
fpath8 = "/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_B8.TIF"
#8bit出力データ名を指定しておく
band8_8bit_path="Band8_8bit.tif"
#切り出しの詳細（ここでは、マルチカラーの2倍を指定することで同じ場所を切り出すことが可能）
minX=4800
minY=9000
deltaX=3000
deltaY=3000
ds=gdal.Translate(band8_8bit_path,fpath8,outputType=gdal.GDT_Byte,scaleParams=[[5000,15000]],srcWin=[minX,minY,deltaX,deltaY] )
ds=None
パンシャープン処理を行うコマンドは、以下になります。
gdal_pansharpen(["空白",1番目のバンドのファイル名,2番目のバンドのファイル名,3番目のバンドのファイル名,パンクロマティックバンドのファイル名])もしくは、すでにカラーバンドを作成している場合は、以下でも可能です。
gdal_pansharpen(["空白",RGB画像のファイル名,パンクロマティックバンドのファイル名])今回は、順番にバンドを入れる方法を紹介します。 1番目に赤のバンド、2番目に緑のバンド、3番目に青のバンドを入れることで、トゥルー画像のパンシャープンを作成することができます。
#パンシャープン画像のファイル名を指定します。
pansharpen_path="PANSHARPEN.tif"
#ds=gdal_pansharpen(["",band8_8bit_path,out_Natural_path,pansharpen_path])
ds=gdal_pansharpen(["",band8_8bit_path,band4_8bit_path,band3_8bit_path,band2_8bit_path ,pansharpen_path])
ds=None
このような操作だけで、パンシャープン画像を作成することが可能です。
ではここで、パンシャープンする前のトゥルーカラー画像と比較してみましょう。
今回作成した画像では、マルチ画像との違いが分かりづらいため、両方の画像について任意の場所を切り出して表示します。
#パンシャープン用切り出しの詳細
minX=1500
minY=1800
deltaX=300
deltaY=200
cut_pans="cut_pans.tif"
ds=gdal.Translate(cut_pans,pansharpen_path,srcWin=[minX,minY,deltaX,deltaY] )
ds=None
#マルチ画像用切り出しの詳細
minX2=750
minY2=900
deltaX2=150
deltaY2=100
cut_RGB="cut_RGB.tif"
#ds=gdal.Translate(cut_RGB,out_True_path,srcWin=[minX2,minY2,deltaX2,deltaY2] )
ds=gdal.Translate(cut_RGB,"TrueColor_TOKYO.tif",srcWin=[minX2,minY2,deltaX2,deltaY2])
ds=None
#画像の表示
plt.figure(figsize=(20,20))
plt.subplot(1,2,1)
image1 = mpimg.imread(cut_RGB)
plt.imshow(image1)
plt.title("multi image TDL")
plt.subplot(1,2,2)
image2 = mpimg.imread(cut_pans)
plt.imshow(image2)
plt.title("pansharpen image TDL")
plt.show
マルチ画像と比較して、パンシャープン画像は解像度が高くシャープに映っていることが確認できます。
gdal_pansharpenの機能を用いて、パンシャープン画像の作成方法を学ぶことができました。
フォーマットの変換¶
最後に、出力したカラー画像（GeoTIFフォーマット）をPNGフォーマットに変換する方法を覚えましょう。
資料に使うために、ファイルサイズを小さくする方法も覚えましょう。
gdal.Translateを用いて、以下のように書きます。
gdal.Translate({出力画像名},{入力画像名},format={フォーマット名},widthPct={X方向の拡大/縮小パーセンテージ},heightPct={Y方向の拡大/縮小パーセンテージ})試しに、先ほど作成したパンシャープンの画像をPNGフォーマットに変換してみましょう。
PNG_img_path="Pansharpen.png" #出力のファイル名を定義
gdal.Translate(
PNG_img_path, #出力画像名
cut_pans, #入力画像名
format='PNG'#出力ファイルのフォーマット指定
)
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
plt.figure(figsize=(10,10))
PNGimg= mpimg.imread(PNG_img_path)
plt.imshow(PNGimg)
plt.title("TDL pansharpen PNG")
plt.show
GeoTIFF画像を"Pansharpen.png"というファイル名で、PNGフォーマットに変換することができました。
その他、JPEGやENVIフォーマット等、色々なフォーマットに変換することが可能です。
コラム2：gdal_mergeによるカラー合成¶
衛星画像を解析する際、良く使うGDALの機能として"gdal_merge"を紹介します。
例えば、複数の位置情報を持った画像データを1枚の画像として結合したいとき、gdal_mergeが有効です。 また、演習と紹介したカラー画像作成のように、同じ観測データの異なるバンドを結合してカラー画像を作成することができます。
gdal_merge.pyの詳しい説明は、以下サイトをご確認ください。
https://gdal.org/programs/gdal_merge.html
ここでは、gdal_mergeを使ったカラー画像の作成方法をご紹介します。
gdal_merge.pyはGDALライブラリの標準では無いため、Python内で使用するためには、機能を呼び出す必要があります。
まずは、自身の環境から、gdal_mergeがどこにあるか以下のように確認します。
端末の中に、gdal_merge.pyがどこにあるか、"which"コマンドで確認しましょう。
!which gdal_merge.py
"/usr/bin"の下にあることが分かったので、以下のようにimportを行います。確認ができたら、次のようにgdal_mergeをインポートします。
import sys
sys.path.append('/usr/bin/')
import gdal
import gdal_merge
まず、本編と同様に青、緑、赤バンドを用意し、8bitの画像に変換します。 ここでは、画像全体をRGB合成してみましょう。
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"
#出力画像の名前を指定
Blue_path_8bit ="8bitTOKYO_Blue.tif"
Green_path_8bit="8bitTOKYO_Green.tif"
Red_path_8bit ="8bitTOKYO_Red.tif"
#gdal.Translate(｛出力画像名｝,{入力画像名},srcWin=[minX,minY,deltaX,deltaY])
ds=gdal.Translate(Blue_path_8bit , fpath2,outputType=gdal.GDT_Byte,scaleParams=[[9700,13400]])
ds=gdal.Translate(Green_path_8bit, fpath3,outputType=gdal.GDT_Byte,scaleParams=[[8100,12800]])
ds=gdal.Translate(Red_path_8bit , fpath4,outputType=gdal.GDT_Byte,scaleParams=[[7000,13000]])
ds=None
gdal_mergeのseparateオプションを用いて、以下のように書きます。
gdal_merge.main(["","-o",{出力画像},"-separate",{入力画像１},{入力画像２},{入力画像３}])早速、実行してみましょう。
gdal_merge.main(["","-o","RGB_KANTO2.tif","-separate",Red_path_8bit,Green_path_8bit,Blue_path_8bit])
これで、カラー画像の作成ができました。画像を見てみましょう。
plt.figure(figsize=(10,10))
img= mpimg.imread("RGB_KANTO2.tif")
plt.imshow(img)
plt.title("KANTO")
plt.show
まとめ¶
この章では、主にPythonのGDALライブラリを用いた、LANDAST画像の加工方法について学習しました。 皆さんも、ご自身の興味のある場所について、自由に衛星画像を加工してみてください。特に、LANDASATは40年前からのデータが豊富にあるため、新しい発見に出会えると思います。