衛星画像より街の開発状況を推測してみる．

概要

　これまでSentinel衛星画像のハブであるCopernicus hubより観測画像を取得し，東京ディズニーランドの駐車場の駐車状況や，東京湾の船舶数のトレンドを求めてきました．

衛星画像より駐車場の利用状況を調べてみた．
衛星画像より東京湾の船数のトレンドを求めてみた．

今回は，前回と同様にGoogle Earth Engineより衛星画像を入手し，東京近郊のこの5年間の開発状況を衛星画像より推測してみます．

Google Earth Engine（以下，GEE)の利用方法は，上記のサイトをご参考にしてください．GoogleのアカウントがあればGoogle colaboratoryで利用することがでるため，大変な環境構築を行う必要もないので比較的簡単に体験することができます．
　
こちらが，Sentinel-1の2020年9月12日の観測画像と，およそ5年前の2015年9月21日の観測画像の合成画像になります．過去の画像を赤色に，最近の画像を青色にしています．
　これより，お台場や東京ディズニーランドに比較的大きな青色の部分がいくつか見られます．これらの場所はこの5年で開発が進んだと推測されます．

　この合成画像の説明も含めて作成方法をご紹介します． 最後のアウトプットは上記の衛星画像をGoogle Earthにて閲覧するためのkmzファイルを作成します．このファイルはこちらにありますので，お手元のGoogle Earthでまずは試したい方はこちらをダウンロードしてください．
　また，記事で紹介したコードは Githubに置きましたので，ご自身で試したい方はこちらをご利用ください．Google colaboratoryを用いていますので，ネットワーク環境さえあればPCによらずご利用できます．

#1.はじめに
　今回は，衛星画像の中から，電波センサと呼ばれるSAR衛星が撮影した画像を用います．
SAR衛星の撮像画像は，構造物の難しいのですが，多数の構造物があるかどうかは定性的には評価することができ，かつ雲の影響を受けないため安定して求めることができます．
　ただし，同じSAR画像でも分解能に差があり，できれば高分解能なSAR画像を利用したいのですが，今回も無料でありかつ定期的に観測しているSentinel-1衛星の観測画像を用います．
　SAR画像による構造物の密集具合は画像からある程度推量することができます．絵心がないので，関連する記事を参照させてもらいます．

衛星データのキホン~分かること、種類、頻度、解像度、活用事例
~

　SAR衛星は電波を地球に照射し，その反射信号の強度を画像化しています．このとき，電波を地表面に対して斜め方向から放射するため，図のように地表面が平坦だと多くのものは電波の入射方向，つまり放射した電波を受信する衛星とは逆方向に反射します．そのため，対象の地表面が平坦であれば受信する信号が弱いため暗い画像となります．湖や海などの水面は極めめて平坦なためより暗く写ります．
　一方，建物やビル，車などの構造物がある場合は，対象物が凹凸しているため，入射方向と同じ方向へ反射する電波が多くなり，その為明るく写ります．
　例えば，東京ディズニーランドの電波衛星の観測画像は以下のようになります．

　ディズニーランドでの施設がよくわかりますね．また，ディズニーランドに行かれた方がご存じだと思いますが，周囲に駐車場があり，そこは広く平坦なため暗く写っているのがわかるかと思います．
　ディズニーランドのホームページをみると， 2015年は駐車場であった場所に東京ディズニーシーの大規模拡張プロジェクトの開発予定となっており，現時点でここにアトラクションが建設されていることがわかりました．その他には，お台場や湾岸部は東京オリンピックなどのイベントなどにあわせて開発が進められており，その模様が衛星画像よりわかるかもしれません．
　今回利用する人工衛星は欧州宇宙機関が開発運用しているSentinel-1と2になります．それぞれの詳細については，以下の記事をご参考下さい．

無料で最新の衛星画像を入手する方法．

解析評価の観測画像はSAR衛星であるSentinel−1です．
Sentinel-1は2機の衛星で構成されており，1機では12日間の回帰日数（同一地域を同じ条件で通過するのが12日毎，という意味）ですが，２機構成のため多くの地域を６日間に１回観測しています．Sentinel衛星の大きな特徴であり利用する側の利点は，この定期観測になります．

[Sentinel-1 Observation Scenario]
(https://sentinel.esa.int/web/sentinel/missions/sentinel-1/observation-scenario)

また，人工衛星の観測時刻はほぼ決まっており，光学衛星であれば午前10時30分頃，SAR衛星であれば午前6時頃と午後6時頃になります．
人工衛星から人は見える？~衛星別、地上分解能・地方時まとめ~

先程紹介しましたSentinel-1は，A衛星が午前6時頃が，B衛星が午後6時頃が観測時刻になります．それぞれの衛星で電波の照射方向が異なるため，両方の衛星画像を用いて比較することが，その差異が捉えられるため適切ではありません．そこで，今回はB衛星の画像のみを用いて，最近の衛星画像と比較する過去の画像を取得し，その合成画像を作成します．合成画像を作成する目的は，その変化を視覚的にわかりやすくするためです．最近の画像を青色，過去の画像を赤色にすることで，青く写る部分は信号が強くなったことを，赤く写る場所は信号が弱くなったことを示す．信号の強弱は建物の有無および密集度に依存するため，青い部分は開発が進んだことを，赤い部分は構造物がなくなったこと（平坦化したこと）を示しています．
　次に，GEEより衛星画像の取得，およびその評価について紹介します．

#２．衛星画像の取得と評価
##2.1　環境準備（構築）．
　これまでも紹介していますが，この記事で初めての方もいらっしゃるため，重複しますが載せておきます．ご存知の方は　2.4までスキップしてください．
　GEEでの衛星画像の取得方法については，以下の記事で詳しく説明されていますのでご参考ください．
　Google Earth EngineとGoogle Colabによる人工衛星画像解析　〜無料で始める衛星画像解析（入門編）〜
Googleのアカウントさえあれば無料で簡単に利用できます．すごい世界になりました．
　取得した衛星画像などのデータはGoogle Driveに保存します．Google Colaboratoryは利用毎にコードやデータが削除されるため，ご自身のGoogle Driveにデータがあると便利です．ただし，無料利用のGoogle dirveの容量は20GBですので，大きな衛星画像をダウンロードされるとすぐに不足しますので，適宜削除するなどご注意ください．
　では，コードを紹介します．

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":"56641.31 m","drawtype":"rectangle","area":"38085.35 ha"},"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[139.7144973278046,35.51642373317183],[139.7144973278046,35.6731886200871],[139.95559573173526,35.6731886200871],[139.95559573173526,35.51642373317183],[139.7144973278046,35.51642373317183]]]}}]}

この地理情報をGEEへの入力フォーマットに，そしてFoliumでの表示用に加工します．

#今後使用する任意のファイル名をセットする． 例えば，地域の名前など．
object_name = 'TokyoBay_RGB'

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]

area = pd.DataFrame(AREA,　columns=['longtitude', 'latitude'])

area_d =[[area['longtitude'].min(), area['latitude'].max()],
 [area['longtitude'].max(), area['latitude'].max()],
 [area['longtitude'].max(), area['latitude'].min()],
 [area['longtitude'].min(), area['latitude'].min()],
 [area['longtitude'].min(), area['latitude'].max()]]

AREA = area_d

それでは，設定した関心域を確認します．

m = folium.Map([(AREA[0][1]+AREA[len(AREA)-2][1])/2,(AREA[0][0]+AREA[len(AREA)-3][0])/2], zoom_start=11)

folium.GeoJson(str(object_name) +'_2.geojson').add_to(m)
m

出力

##2.3　GEEから衛星画像の取得
　GEEには多くの衛星画像や，すで解析された多くの情報がセットされています．詳しくはデータカタログをご参考んください．Sentinel-1と2は以下となります．

Sentinel-1 SAR GRD: C-band Synthetic Aperture Radar Ground Range Detected, log scaling
Sentinel-2 MSI: MultiSpectral Instrument, Level-1C

　このページから，Sentinel-1の観測画像は2014年10月3日から，Sentinel-2は2015年6月23日からのデータが利用できます．Sentinel-2の観測画像は，大気のモヤを削除した大気補正のLevel2Aの画像も用意されていますが，標準的に解析がされるようになった2017年3月28日以降の画像のみが対象となります．そのため，より古い観測画像を用いたい場合は，今回の画像をご利用ください．

　では，GEEよりSentinel-1,2の画像を取得し，Google Colaboratoryに保存します．

　最初に，GEEに設定する地理情報のフォーマットを準備します．

region=ee.Geometry.Rectangle(area['longtitude'].min(),area['latitude'].min(), area['longtitude'].max(), area['latitude'].max())

　次に，取得する情報のパラメータを設定します．今回は，取得する画像の期間と，取得した画像の保存先を指定しています．
　今回は，異なる日の衛星画像を取得するため，まずは最近の衛星画像を取得します．ここでは，2020年9月20日としています．

#期日を指定
import datetime

date = datetime.date(2020, 9, 20)

dt_now = datetime.date.today()
dt = dt_now - date

if dt.days < 6:
  dt_nd = datetime.timedelta(days=12)
  dt_pd = datetime.timedelta(days=0)
else:
  dt_nd = datetime.timedelta(days=6)
  dt_pd = datetime.timedelta(days=6)

# 　観測期間
from_date= date - dt_nd
to_date= date + dt_pd

from_date = from_date.strftime('%Y-%m-%d')
to_date = to_date.strftime('%Y-%m-%d')


# 保存するフォルダ名
dir_name_s1 = 'GEE_Sentinel1_' + object_name
dir_name_s2 = 'GEE_Sentinel2_' + object_name

　ここで，Sentinel-1の観測頻度は12日に1回のため，設定した日が解析時に近い場合は過去の画像を，ある程度前の画像であればその前後6日間の画像を取得することとしました．

　では，Sentinel-1と2の画像条件の設定です．今回はSentinel-2の画像を用いませんが，参考のために取得します．

def cloudMasking(image):
    qa = image.select('QA60')
    cloudBitMask = 1 << 10  
    cirrusBitMask = 1 << 11
    mask = qa.bitwiseAnd(cloudBitMask).eq(0).And(qa.bitwiseAnd(cirrusBitMask).eq(0))
    return image.updateMask(mask).divide(10000)

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()

#同一衛星での観測画像を比較するために，AscendingかDecsendingのどちらかのみとする．
Sentinel1 = ee.ImageCollection('COPERNICUS/S1_GRD').filterBounds(region).filterDate(parse(from_date),parse(to_date)).filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING')).select(['VV'])
Sentinel2 = ee.ImageCollection('COPERNICUS/S2').filterBounds(region).filterDate(parse(from_date),parse(to_date)).filterMetadata('CLOUDY_PIXEL_PERCENTAGE','less_than', 10).map(cloudMasking).select(['B4','B3','B2'])

imageList_s1 = Sentinel1.toList(300)
imageList_s2 = Sentinel2.toList(300) 

　ここで，Sentinel-1の観測画像は，B衛星の午後6時の観測画像のみを用いるため，'orbitProperties_pass'で **'ASCENDING'**を選択しています．”Descending”にすると午前6時の観測画像になります．

　では，上記の衛星画像の取得条件で，Sentinel-1 およびSentinel-2の関心域の画像を取得します．

for i in range(imageList_s1.size().getInfo()):
    image = ee.Image(imageList_s1.get(i))
    ImageExport(image.reproject(crs='EPSG:4326',scale=10),image.get('system:index').getInfo(),dir_name_s1,region['coordinates'][0],10)

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　衛星画像の表示とフィルターによる平滑化処理
　それでは，2.3で取得した衛星画像を表示し確認します．
　衛星画像はGoogle driveのmy driveに設定したディレクトリ（フォルダ）に保存されています．それを呼び出し表示します．最初にSentinel-2の観測画像からです．

# 時系列で可視化
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.float32)*2 #輝度を2倍に明るさを補正．
  plt.imshow(arrayImg)
  plt.title(files[i][0:8])# ファイル名から日付を取得
  #plt.tight_layout()

　Sentinel-2の観測画像の分解能は10mですが，大きな構造物であれば建設されたかどうかわかるかもしれません．光学衛星の観測画像はGoogle Mapにて普段から見慣れているため，このサービスを利用するときは目的（の場所）は特定されているので簡単に見つけることができます．ただし，これはそこに構造物が建てられたであろう，という情報があれば注目してみることができますが，この範囲のどこかに構造物が建設されているか見つけて，では指標がないため難しくなります．
　そこで，今回はSentinel-2の観測画像を評価には利用せず，Sentinel-1の観測画像を比較した変化を示すことで，差がある箇所を容易に見つけることができます．
　では，Sentinel-1の観測画像を表示します．

# 時系列で可視化
s1_path = '/content/drive/My Drive/' + dir_name_s1 + '/'
files =os.listdir(s1_path)
files.sort()

plt.figure(figsize=(20, 40))
j=0


v = len(files)//5 +1 
for i in range(len(files)):
  # 画像を1シーンずつ取得して可視化
  with rasterio.open(s1_path + files[i]) as src:
      arr = src.read()
  #print(arr[0].shape)
  j+=1# 画像のプロット位置をシフトさせ配置
  plt.subplot(v,5,j)
  plt.imshow(arr[0], cmap='gray')
  plt.title(files[i][33:41])# ファイル名から日付を取得
  date0 = files[i][33:41]
  plt.tight_layout()

　これより，同一の観測日の2つの観測画像が取得されています．同じエリアで2つの観測画像に別れるのは，Sentinel-1の衛星画像の標準サイズが100km☓100kmであり，今回の関心域がそれを横断していたため，と思われます．
　そこで，2つに分割した衛星画像をまずは一つにまとめます．

# データの読み込み
n = 0

with rasterio.open(s1_path + files[n]) as src:
    arr = src.read()

print(files[n][33:41])
# 可視化
plt.imshow(arr[0], cmap='gray')

type(arr[0])
image1 =np.zeros(arr[0].shape)

image_1 = np.nan_to_num(arr[0])

# データの読み込み
n = 1

with rasterio.open(s1_path + files[n]) as src:
    arr = src.read()

print(files[n][33:41])
# 可視化
plt.imshow(arr[0], cmap='gray')

type(arr[0])
image2 =np.zeros(arr[0].shape)

image_2 = np.nan_to_num(arr[0])

　ここで，関心域内で取得されていな領域（白い部分）は，データなし**（Nan)となっているため，ここを0**に置き換えています．次に，この2つの画像を結合します．

image_BG = image_1 + image_2
plt.imshow(image_BG, cmap='gray')

　これで，2020年9月12日のSAR衛星の衛星画像を準備することができました．
　SAR衛星の画像は地上の構造物からの反射によってスペックルノイズとよばれる多数の明るい点のノイズが生じます．これはSAR衛星で用いる電波が波であるため，反射時の干渉によって生じるものです．このノイズがこのあとの処理に影響することがあるため，平滑化処理が必要となります．今回は，一般的なSAR画像の平滑化処理でありLeeフィルタを用いました．準備したLeeフィルタの関数は以下となります．

from scipy.ndimage.filters import uniform_filter
from scipy.ndimage.measurements import variance

def lee_filter(img, size):
    img_mean = uniform_filter(img, (size, size))
    img_sqr_mean = uniform_filter(img**2, (size, size))
    img_variance = img_sqr_mean - img_mean**2

    overall_variance = variance(img)

    img_weights = img_variance / (img_variance + overall_variance)
    img_output = img_mean + img_weights * (img - img_mean)
    return img_output

　では，平滑化処理を行います．ここでは，3ピクセル内の平滑化処理を行いました．

image_BG_lee = lee_filter(image_BG, 3)

plt.imshow(image_BG_lee,cmap='gray')

　この衛星像のサイズではフィルタの効果はわかりませんね．作成した画像をダウンロードするなどして，その差異を確認してみてください．
　次に，比較する過去の衛星画像を取得します．ここでは5年前の撮像画像としました．

#ベース画像の観測日を指定
date = datetime.date(2015, 9, 20)
print(date)

# 　観測期間

from_date= date - dt_nd
to_date= date + dt_pd

from_date = from_date.strftime('%Y-%m-%d')
to_date = to_date.strftime('%Y-%m-%d')


# 保存するフォルダ名
dir_name_s1 = 'GEE_Sentinel1_R_' + object_name
dir_name_s2 = 'GEE_Sentinel2_R_' + object_name

for i in range(imageList_s1.size().getInfo()):
    image = ee.Image(imageList_s1.get(i))
    ImageExport(image.reproject(crs='EPSG:4326',scale=10),image.get('system:index').getInfo(),dir_name_s1,region['coordinates'][0],10)

　では，取得した画像を表示してみます．

# 時系列で可視化
s1_path = '/content/drive/My Drive/' + dir_name_s1 + '/'
files =os.listdir(s1_path)
files.sort()

plt.figure(figsize=(20, 40))
j=0


v = len(files)//5 +1 
for i in range(len(files)):
  # 画像を1シーンずつ取得して可視化
  with rasterio.open(s1_path + files[i]) as src:
      arr = src.read()
  print(arr[0].shape)
  #print(image)
  j+=1# 画像のプロット位置をシフトさせ配
  plt.subplot(v,5,j)
  #image =  arr[0].fillna(0)
  plt.imshow(arr[0], cmap='gray')
  plt.title(files[i][33:41])# ファイル名から日付を取得
  date1 = files[i][33:41]
  plt.tight_layout()

　今回は，先程と異なる画像は分割されておらず，関心域のすべてを含めた一枚の画像となっています．
　こちらの画像も，先程と同じくフィルタにて平滑処理を行います．

image_R_lee = lee_filter(image_BG, 3)
plt.imshow(image_R_lee,cmap='gray')

　次に，最近の衛星画像と過去の衛星画像を合成したRGB画像を作成します．今回は，過去画像を赤色に，最近の画像を青色および緑色としています．

##2.5　最近と過去の衛星画像のRGB合成画像を作成
　では，それぞれの衛星画像よりRGB合成画像を作成します．

#RGB合成
RGB = np.dstack((image_R, np.dstack((image_BG_lee, image_BG_lee))))

print(RGB.shape)
plt.figure(figsize=(15, 20))
plt.imshow(RGB)
plt.show()

　この画像よりお台場や舞浜で青くなっている場所が確認できます．東京湾の赤と青の点は船舶を示しており，撮影日の船が写っています．他には赤く写っているところがあり，これは開発によって建物がなくなるか集合して一つの大きなものになったのかと思われます．赤色はその意味でも開発を示しているとも言えますね．

　では，ここで取得した画像を観測日やクレジットを記載しjpegで保存します．

import cv2
from PIL import Image, ImageDraw, ImageFont

img_bgr = np.dstack((image_R, np.dstack((image_BG_lee, image_BG_lee))))

new_arr = ((img_bgr - img_bgr.min()) * (1/(img_bgr.max() - img_bgr.min()) * 255)).astype('uint8')


im_rgb = cv2.cvtColor(new_arr, 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()

　このままでは全体に白い画像となっており，変化している青や赤が見づらくかんじます．ここで，下記を実行しコントラストを変更します．ここはお好みですね．

import PIL.Image
from PIL import ImageEnhance

IMAGE_PATH = str(object_name) +'.jpg'
CONTRAST = 2.0

img = PIL.Image.open(IMAGE_PATH)

# コントラストを変える
contrast_converter = ImageEnhance.Contrast(img)
contrast_img = contrast_converter.enhance(CONTRAST)

#　画像を保存
contrast_img.save(str(object_name) +'.jpg')

plt.figure(figsize=(15, 20))

plt.imshow(contrast_img)
plt.show()

　先程よりも青と赤の変化が見やすくなっていすね．
　次に，この画像に衛星画像の撮影日およびクレジットを記載します．ここで使用するフォントはネットにてフリーで公開されているものを利用します．これもお好みですね．

#フォントファイルのダウンロードと設定
!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"

　では，画像に文字を書き込みます．

img = Image.open(str(object_name) +'.jpg')

img = img.convert('RGB')


x = int(img.size[0]/1.21) #日付の記載位置の設定
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), 'R: '+str(date1), fill=(255, 255, 255), font=obj_font)
obj_draw.text((x, y+1.1*fs), 'G: '+str(date0), fill=(255, 255, 255), font=obj_font)
obj_draw.text((x, y+2.2*fs), 'B: '+str(date0), fill=(255, 255, 255), font=obj_font)
obj_draw.text((img.size[0]/2.0, 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 = img.convert('L')


img.save(str(object_name) +'.jpg')

plt.figure(figsize=(15, 20))

plt.imshow(img)
plt.show()

　これで画像に必要な情報を記載することができました．文字の位置やサイズ，色などは上記のコードを変更して，好みに調整してください．
　次に，ここで取得した画像をローカルのPCでGISアプリ(例えば，Google Earth)で閲覧できるようkmzファイルを作成します．

##2.6　GISソフトで閲覧するためのkmzファイルを作成
　kmzファイルの作成方法は，以前詳しく紹介した以下の記事をご参考にしてください．

人工衛星画像の無料環境なプラットフォーム「Tellus」にてデータ解析． -Google Earth編-

　kmzファイルの作成方法はいくつかあると思いますが，今回はpythonに提供されているsimplekmlを使います．モジュールをインストールし使用します

!pip install simplekml

import simplekml
import zipfile

　次に，画像の緯度・経度情報を読み込みます．

area = pd.DataFrame(AREA,
                  columns=['longtitude', 'latitude'])

north = area["latitude"].max()
south = area["latitude"].min()
east = area["longtitude"].max()
west = area["longtitude"].min()

　次に，上記の関心域のkmlファイル(緯度・経度情報）を作成します．

#地理情報をkmlで出力．
kml = simplekml.Kml()
ground = kml.newgroundoverlay(name=str(object_name))
ground.icon.href = str(object_name) +'.jpg'
ground.latlonbox.north = north
ground.latlonbox.south = south
ground.latlonbox.east = east
ground.latlonbox.west = west
ground.latlonbox.rotation = 0
kml.save(str(object_name)+".kml")

　そして，作成したkmlファイルと閲覧する画像とのリンクを定義し，kmzファイルを作成します．

#地理情報と画像とをあわせたkmzを出力．
with zipfile.ZipFile(str(object_name)+'.kmz', "w", zipfile.ZIP_DEFLATED) as zf:
    zf.write(str(object_name)+".kml", arcname=str(object_name)+".kml")
    zf.write(str(object_name) +'.jpg', arcname=str(object_name) +'.jpg')

　これで，kmzファイルが作成されました．Google Colaboratoryで作成したkmzファイルが確認できます．このファイルの右側をマウスで選択するとdownloadが表示されますので，ファイルをダウンロードします．

　端末にGoogle Earthがインストールされていましたら，ダウンロードしたkmzファイルをダブルクリックすると自動でGoogle Earthが起動し，地図上に作成した画像が貼り付けられます．

　画像から青や赤くなっている場所を確認し，kmzファイルの選択を外すと該当する場所がどこか確認できます．お台場の青い部分はコンテナの置き場所？であることがわかりました．それ以外には，東京ディズニーランドの駐車場で一部が青くなっており，ここは計画されていたディズニーシーの新しいアトラクションの建設場所と推測されます．（Google Earthの航空写真は時間スライダーがありますので，該当箇所の撮影日や過去のものなど確かめてみてください．）

#3. 最後に
　Google が提供するGoogle Earth Engineを用いて，希望する日の衛星画像を取得し，街がどのように変わっているのか，大きく変化している場所を確認する方法を紹介しました．
　ここでは，東京湾を中心に評価し，お台場や東京ディズニーランドで積極的に開発が進んでいるのがわかりました．同じ方法にて，別の場所を見てみると発展している都市が抽出できるかと思います．また，街の状況の時間変化を見えるかしていますので，開発以外にも災害要因による被害が同様の方法でも識別することができます．

　災害時の人工衛星活用ガイドブック水害版・衛星基礎編

　個人的には，衛星画像の利用の一番の方法は，このような変化を抽出すること，と思っています．衛星画像は面的な情報であり，過去と現在の画像としての差や，その間でのトレンド（例えば，以前の記事にて紹介した駐車場の駐車数など）を長期にわたり得ることができます．目的は色々ありますが，方法は変化抽出になります．衛星画像の取得は閲覧だけでなく，過去と現在などの比較を行い，何が起こっているのか推測してみるのも楽しいです．
　これを機会に衛星画像に関心を持つ方が増えると嬉しいです．
　ご意見や質問などありましたら，お気軽にコメントください．

##参考記事
衛星画像より駐車場の利用状況を調べてみた．
衛星画像より東京湾の船数のトレンドを求めてみた．
無料で最新の衛星画像を入手する方法．

衛星データのキホン~分かること、種類、頻度、解像度、活用事例
~
人工衛星から人は見える？~衛星別、地上分解能・地方時まとめ~

Sentinel-1 SAR GRD: C-band Synthetic Aperture Radar Ground Range Detected, log scaling
Sentinel-2 MSI: MultiSpectral Instrument, Level-1C

　災害時の人工衛星活用ガイドブック水害版・衛星基礎編