(更新)合成画像の作成方法を更新．
概要
これまで，Google Earth Engineを用いて欧州の人工衛星 Sentinelシリーズの衛星画像の取得，およびその解析方法について紹介してました．今回は，Sentinelシリーズでも皆さんに馴染みのある光学観測衛星であるSentinel-2の衛星画像の取得，および解析方法についてご紹介します．
今回は，海外ニュースで報道されているアメリカ西海岸の森林火災を見てみます．
こちらは，アメリカ西海岸のポートランド近郊の森林火災現場の衛星画像（2020年9月9日撮像）になります．
海外のニュースサイトでは以下のように報道されています．こちらのサイトでは，森林火災の履歴があり，どこで火災があったのか把握することができます．
Fire Maps: Glass and Zogg Wildfire Tracker
現場の様子は以下のサイトがリアルでせまってくるものがあります． New York Timesは興味深い記事が多く毎日みています．
Historic Wild Fires Range in Western States
ここでは，Google Earth EngineでのSentinel-2の観測画像の取得方法や，その画像データをカラー画像にする方法例について紹介します．
また，記事で紹介したコードは Githubに置きましたので，ご自身で試したい方はこちらをご利用ください．Google colaboratoryを用いていますので，ネットワーク環境さえあればPCによらずご利用できます．
#1.はじめに
以前，SentinelのAPIを用いて，2020年初頭のオーストラリアの森林火災の人工衛星による撮像画像を取得し，そのGIFアニメーションを作成しました．
オーストラリア 森林火災の衛星画像のタイムラプスを作ってみた．
ここでは，Sentinel-2の標準的な100km四方の撮像画像をAPIを用いて取得し，それより関心部分を切り出しカラー合成処理を行っていました．このとき，撮像タイミングによっては異なるシーンとなるため，一部切り欠きがある画像となっています．異なるシーンであっても，それぞれを取得し結合（Merge)することで関心域全体の画像にすることはできます．しかし，複数の大きなサイズの画像を扱うことは，取得時間を考えるとあまり効率的ではありません．そこで，今回はGoogle Earth Engine(以下，GEEという）より，関心域に切り出された画像を取得しカラー合成画像を作成します．
Sentinelシリーズは，欧州宇宙機関が開発・運用する地球観測シリーズで，すでに１から５まで打上がり運用されています．
今回は，デジカメ画像と同じく光学観測画像を撮像するSentinel-2の衛星画像の取得から，画像処理を紹介します．
Sentinel-2 ©ESA 2000-2018.
Sentinel-2について，その観測計画や観測画像の詳細については，無料で最新の衛星画像を入手する方法．をご参考にしてください．
Sentinel-2は１２の異なる波長の画像を撮像しています．
(Credit: European Space Agency)
一般的な衛星画像（True color)であれば，バンドの4，3，2をR,G,Bとすることで作ることができます．
今回は，森林火災にフォーカスし，赤外波長の観測画像をもちいて，火災の可視化を試してみます．
今回は，以下の記事を参考にしています．
この記事より，火災に感度のある長波長のバンドである12と11，および森林に感度がある8Aの画像を用いて，森林火災現場を見てみます．
#２．衛星画像の取得と評価
##2.1 環境準備（構築）．
GEEでの衛星画像の取得方法については，以下の記事で詳しく説明されていますのでご参考ください．
Google Earth EngineとGoogle Colabによる人工衛星画像解析 〜無料で始める衛星画像解析（入門編）〜
また，GEEより入手した衛星画像の解析方法については，以下の記事をご参考ください．
重複する部分もありますが，衛星画像の取得，およびその処理を紹介します．
import ee
import numpy as np
import matplotlib.pyplot as plt
ee.Authenticate()
ee.Initialize()
最初にこちらを実行し，GEEの接続認証を行います．実行するとリンクが帰ってきますので，そちらをクリックし認証手続きを行い，アクセスコードをコピーして入力ください．
次に，Google Driveへの接続認証を行います．こちらも，GEEの認証とフローは同じですね．
from google.colab import drive
drive.mount('/content/drive')
次に，取得した衛星画像の閲覧や数値化して解析するために必要なモジュールのインストール等の作業を行います．
# パッケージのインストール&インポート
!pip install rasterio
import numpy as np
import matplotlib.pyplot as plt
import rasterio
import json
import os
import glob
import time
from datetime import datetime
from dateutil.parser import parse
よく利用されるモジュールはGoogle Colaboratoryにすでにインストールされているために追加の作業は必要ありませんが，今回は地図情報が付加された画像であるGeotiffを用いるため，画像処理のために必要なRasterioをイントールします．
次に，設定した対象地域を地図上で確認するためにfoliumというモジュールをインストールします．
!pip install folium
import folium
これで環境の準備ができましたので，衛星画像をGEEより取得します．
##2.2 関心域（対象域）の設定
GEEから衛星画像を取得するためには，自分が関心がある対象地域の緯度・経度情報を入力する必要があります．
今回も，簡単に関心域の緯度経度が調べかつ入手するために作成した以下のサイトを用いて取得します．
#関心領域のポリゴン情報の取得．
from IPython.display import HTML
HTML(r'<iframe width="1000" height="580" src="https://gispolygon.herokuapp.com/" frameborder="0"></iframe>')
こちらを実行すると，以下の画面が表示されます．
関心地域を拡大した後，左のアイコンより四角のポリゴンを選んで，関心域のポリゴンを表示させます．その後，Show Featuresをクリックすると，右のウィンドウにポリゴンの地理情報が表示されます．このあと，下部のCopyをクリックして，この地理情報をコピーします．
今回は米国オレゴン州のポートランドの東側エリアの衛星画像を取得します．
そしてコピーした地図情報を下記にペーストして入力します．
A = {"type":"FeatureCollection","features":[{"properties":{"note":"","distance":"210025.76 m","drawtype":"rectangle","area":"520197.34 ha"},"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[-122.60860204696657,44.80966973731369],[-122.60860204696657,45.400920907537866],[-121.60266637802125,45.400920907537866],[-121.60266637802125,44.80966973731369],[-122.60860204696657,44.80966973731369]]]}}]}
この地理情報をGEEへの入力フォーマットに，そしてFoliumでの表示用に加工します．
#今後使用する任意のファイル名をセットする． 例えば，地域の名前など．
object_name = 'Portland'
with open(str(object_name) +'_2.geojson', 'w') as f:
json.dump(A, f)
json_file = open(str(object_name) +'_2.geojson')
json_object = json.load(json_file)
#jsonから関心域の緯度・経度情報のみを抽出する．
AREA = json_object["features"][0]["geometry"]['coordinates'][0]
それでは，設定した関心域を確認します．
m = folium.Map([(AREA[0][1]+AREA[len(AREA)-2][1])/2,(AREA[0][0]+AREA[len(AREA)-3][0])/2], zoom_start=9)
folium.GeoJson(str(object_name) +'_2.geojson').add_to(m)
m
出力
##2.3 GEEから衛星画像の取得
GEEには多くの衛星画像や，すで解析された多くの情報がセットされています．詳しくはデータカタログをご参考んください．Sentinel-1と2は以下となります．
Sentinel-2 MSI: MultiSpectral Instrument, Level-1C
このページから，Sentinel-2の観測画像は2015年6月23日から，のデータが利用できます．
では，GEEよりSentinel-2の画像を取得し，Google colaboratoryに保存します．
最初に，GEEに設定する地理情報のフォーマットを準備します．
region=ee.Geometry.Polygon(AREA)
次に，取得する情報のパラメータを設定します．今回は，取得する画像の期間と，取得した画像の保存先を指定しています． 今回は同地域を9月9日に観測していることを確認しましたので，その画像のみを取得するよう期間を設定しました．
# 期間を指定
from_date='2020-09-08'
to_date='2020-09-10'
# 保存するフォルダ名
dir_name_s2 = 'GEE_Sentinel2_' + object_name
では，Sentinel-2の画像条件の設定です．
def ImageExport(image,description,folder,region,scale):
task = ee.batch.Export.image.toDrive(image=image,description=description,folder=folder,region=region,scale=scale)
task.start()
Sentinel2 = ee.ImageCollection('COPERNICUS/S2').filterBounds(region).filterDate(parse(from_date),parse(to_date)).filterMetadata('CLOUDY_PIXEL_PERCENTAGE','less_than', 80).select(['B12','B11','B8A'])
imageList_s2 = Sentinel2.toList(300)
ここで，Sentinel-2の観測画像のうち，カラー合成でもちいるB12, B11およびB8Aのカラー合成画像を取得するよう設定しました．また，Sentinel-2の画像内の雲の割合のフィルターができるのですが，火災による煙も雲と誤認識することをさけるため，通常よりも高い80%としています．
では，上記の衛星画像の取得条件で，Sentinel-2の関心域の画像を取得します．
for i in range(imageList_s2.size().getInfo()):
image = ee.Image(imageList_s2.get(i))
ImageExport(image.reproject(crs='EPSG:4326',scale=10),image.get('system:index').getInfo(),dir_name_s2,region['coordinates'][0],10)
##2.4 衛星画像の表示
取得した衛星画像を表示し確認します．
衛星画像はGoogle driveのmy driveに設定したディレクトリ（フォルダ）に保存されています．それを呼び出し表示します．
# 時系列で可視化
s2_path = '/content/drive/My Drive/' + dir_name_s2 + '/'
files =os.listdir(s2_path)
files.sort()
plt.figure(figsize=(25, 25))
j=0
v = len(files)//5 +1
for i in range(len(files)):
# 画像を1シーンずつ取得して可視化
with rasterio.open(s2_path + files[i]) as src:
arr = src.read()
j+=1# 画像のプロット位置をシフトさせ配置
plt.subplot(v,5,j)
arrayImg = np.asarray(arr).transpose(1,2,0).astype(np.float16)/10000
plt.imshow(arrayImg)
plt.title(files[i][0:8])# ファイル名から日付を取得
plt.tight_layout()
同地域を観測した4つの画像を取得しました．4つに分割されているのは，Sentinel-2は100km四方に画像が分割されて登録されており，この関心域には分割された4つの画像が含まれるためです．1番目と2番めの画像を合成することで，関心域がカバーされますので，この2つの画像の合成画像を作成します．
はじめに，それぞれの画像を確認します．
# 画像の確認
s2_path = '/content/drive/My Drive/' + dir_name_s2 + '/'
files =os.listdir(s2_path)
files.sort()
# データの読み込み
n = 0
with rasterio.open(s2_path + files[n]) as src:
arr = src.read()
print(files[n][0:8])
# 可視化
arrayImg0 = np.asarray(arr).transpose(1,2,0).astype(np.float16)/10000
plt.imshow(arrayImg0)
# 画像の確認
s2_path = '/content/drive/My Drive/' + dir_name_s2 + '/'
files =os.listdir(s2_path)
files.sort()
# データの読み込み
n = 1
with rasterio.open(s2_path + files[n]) as src:
arr = src.read()
print(files[n][0:8])
# 可視化
arrayImg1 = np.asarray(arr).transpose(1,2,0).astype(np.float16)/10000
plt.imshow(arrayImg1)
次に，これら2つの画像を結合します．
**（更新）**画像の結合方法は複数の手法があると思いますが，以前はnumpyにて結合をもちいましたが，コメントにてご助言いただいたrasterioのmergeコマンドにより合成します．
import rasterio.merge
dest, out_transform = rasterio.merge.merge([s2_path + files[0], s2_path + files[1]])
合成画像は，RGBチャンネル，Vertical,Horizontalの順番のためtranspose関数にて変換します．
image = dest/10000
img_transformed = image.transpose((1, 2, 0))
img_transformed.shape
出力
(6596, 11199, 3)
画像を出力します．
plt.imshow(img_transformed)
うまく合成できていますね．（更新前は結合箇所に黒い部分ができていました．）
では，この合成画像をjpegにて保存します．
import cv2
new_image = ((img_transformed - img_transformed.min()) * (1/(img_transformed.max() - img_transformed.min()) * 255)).astype('uint8')
im_rgb = cv2.cvtColor(new_image, cv2.COLOR_BGR2RGB)
cv2.imwrite(str(object_name) +'.jpg', im_rgb )
img = Image.open(str(object_name) +'.jpg')
plt.figure(figsize=(15, 20))
plt.imshow(img)
plt.show()
次に，この画像に人工衛星の撮像日および画像のクレジットを載せます．
##2.5 衛星画像に撮像日およびクレジットを記載する．
次に，この画像に衛星画像の撮影日およびクレジットを記載します．ここで使用するフォントはネットにてフリーで公開されているものを利用します．これもお好みですね．
#フォントファイルのダウンロードと設定
!wget https://osdn.net/dl/mplus-fonts/mplus-TESTFLIGHT-063a.tar.xz
!xz -dc mplus-TESTFLIGHT-*.tar.xz | tar xf -
fontfile = "./mplus-TESTFLIGHT-063a/mplus-1c-bold.ttf"
では，画像に文字を書き込みます．
date = files[n][0:8]
img = Image.open(str(object_name) +'.jpg')
img = img.convert('RGB')
x = int(img.size[0]/1.3) #日付の記載位置の設定
y = int(img.size[1]/20) #日付の記載位置の設定
fs = int(img.size[0]/70) #日付のフォントサイズの設定
obj_draw = ImageDraw.Draw(img)
obj_font = ImageFont.truetype(fontfile, fs)
obj_draw.text((x, y), 'Observation Date: '+str(date), fill=(255, 255, 255), font=obj_font)
obj_draw.text((img.size[0]/1.6, img.size[1]-y*0.1 - img.size[1]/30 ), 'Contains modified Copernicus Sentinel data (2020)', fill=(255, 255, 255), font=obj_font)
img = img.resize((int(img.size[0] / 2) , int(img.size[1] / 2)))
img.save(str(object_name) +'.jpg')
plt.figure(figsize=(15, 20))
plt.imshow(img)
plt.show()
これで画像に必要な情報を記載することができました．文字の位置やサイズ，色などは上記のコードを変更して，好みに調整してください．
作成した画像ファイルは，左のFile一覧より対象ファイル名を右クリックで選択肢，Downloadにて自分のPCに取得することができます．
#3. 最後に
Google が提供するGoogle Earth Engineを用いて，Sentinel-2の衛星画像の取得，および解析例としてアメリカ西海岸の森林火災の画像処理の方法を紹介しました．
森林火災は今も続いており，多大な被害が報告されています．消防隊の皆様の活動が実り，はやく普段の暮らしに戻ることをお祈りします．
ご意見や質問などありましたら，お気軽にコメントください．嬉しいので．
##参考記事
無料で最新の衛星画像を入手する方法．
人工衛星（Sentinel-2)の観測画像をAPIを使って自動取得してみた．
PyTorchによる人工衛星画像から車の推定分布地図を作成してみる．
衛星画像より駐車場の利用状況を調べてみた．
オーストラリア 森林火災の衛星画像のタイムラプスを作ってみた．
【続編】Google Earth EngineとGoogle Colabによる人工衛星画像解析 〜無料で始める衛星画像解析（実践編）〜
Google Earth EngineとGoogle Colabによる人工衛星画像解析 〜無料で始める衛星画像解析（入門編）〜
6. PythonによるローカルからのGEE実行
衛星データのキホン~分かること、種類、頻度、解像度、活用事例
~
人工衛星から人は見える？~衛星別、地上分解能・地方時まとめ~