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

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

[衛星画像より駐車場の利用状況を調べてみた．](https://qiita.com/nigo1973/items/c3be4aca9d882cb5a0e1)
[衛星画像より東京湾の船数のトレンドを求めてみた．](https://qiita.com/nigo1973/items/eeface437e7be513605a)

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

 Google Earth Engine（以下，GEE)の利用方法は，上記のサイトをご参考にしてください．**GoogleのアカウントがあればGoogle colaboratoryで利用することがでる**ため，大変な環境構築を行う必要もないので**比較的簡単に体験**することができます．
　![Screenshot from 2020-09-22 10-44-47.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/90741c0e-8437-81f1-f26e-cbabdae77452.png)
こちらが，Sentinel-1の**2020年9月12日の観測画像**と，およそ5年前の**2015年9月21日の観測画像**の合成画像になります．過去の画像を赤色に，最近の画像を青色にしています．
　これより，お台場や東京ディズニーランドに比較的大きな青色の部分がいくつか見られます．これらの場所は**この5年で開発が進んだ**と推測されます．

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


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

![Screenshot from 2020-08-16 10-34-57.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/a2f77ea6-e39a-555c-9169-59ea59a238a4.png)
[衛星データのキホン~分かること、種類、頻度、解像度、活用事例
~](https://sorabatake.jp/279/)

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

![Screenshot from 2020-08-16 10-41-01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/da5a82de-f99d-69a1-4df2-939aabe88a8c.png)

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

[無料で最新の衛星画像を入手する方法．](https://qiita.com/nigo1973/items/9bb6a11caac8e3e1e850)

 解析評価の観測画像はSAR衛星であるSentinel−1です．
 Sentinel-1は2機の衛星で構成されており，**1機では12日間の回帰日数**（同一地域を同じ条件で通過するのが12日毎，という意味）ですが，２機構成のため多くの地域を**６日間に１回観測**しています．Sentinel衛星の大きな特徴であり利用する側の利点は，この**定期観測**になります．
![Screenshot from 2020-08-16 10-50-29.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/31b45153-f406-18b8-129c-a95972ffd79e.png)
 [Sentinel-1 Observation Scenario]
(https://sentinel.esa.int/web/sentinel/missions/sentinel-1/observation-scenario)

また，人工衛星の観測時刻はほぼ決まっており，**光学衛星であれば午前10時30分頃，SAR衛星であれば午前6時頃と午後6時頃**になります．
[人工衛星から人は見える？~衛星別、地上分解能・地方時まとめ~](https://sorabatake.jp/441/)
![Screenshot from 2020-09-22 11-03-12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/4f044606-6985-e30d-efe9-803b7a929335.png)


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

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

```python
import ee
import numpy as np
import matplotlib.pyplot as plt

ee.Authenticate()
ee.Initialize()
```
　最初にこちらを実行し，GEEの接続認証を行います．実行するとリンクが帰ってきますので，そちらをクリックし認証手続きを行い，アクセスコードをコピーして入力ください．

　次に，Google Driveへの接続認証を行います．こちらも，GEEの認証とフローは同じですね．

```python
from google.colab import drive
drive.mount('/content/drive')
```
次に，取得した衛星画像の閲覧や数値化して解析するために必要なモジュールのインストール等の作業を行います．

```python
# パッケージのインストール&インポート
!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**というモジュールをインストールします．

```python
!pip install folium

import folium
```
これで環境の準備ができましたので，衛星画像をGEEより取得します．

##2.2　関心域（対象域）の設定
　GEEから衛星画像を取得するためには，**自分が関心がある対象地域の緯度・経度情報**を入力する必要があります．
　学校の地理の授業や地球儀から日本の緯度・経度はなんとなく覚えていますが，自宅の緯度・経度は？と言われても，検索して調べてようやくわかると思います．（もちろん，私も覚えていません．）そのため，簡単に関心域の緯度経度が調べかつ入手できるように以下を作ってみました．

```python
#関心領域のポリゴン情報の取得．
from IPython.display import HTML
HTML(r'<iframe width="1000" height="580" src="https://gispolygon.herokuapp.com/" frameborder="0"></iframe>')
```
こちらを実行すると，以下の画面が表示されます．

![Screenshot from 2020-08-16 11-17-18.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/0a0d758f-bbbb-bded-ecee-6ba70a8318da.png)

　関心地域を拡大した後，左のアイコンより四角のポリゴンを選んで，関心域のポリゴンを表示させます．その後，**Show Features**をクリックすると，右のウィンドウにポリゴンの地理情報が表示されます．このあと，下部の**Copy**をクリックして，この地理情報をコピーします．
例えば，今回の対象の東京ディズニーランドの駐車場の場合は下記のようになります．
![Screenshot from 2020-09-22 11-07-07.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/09e789a0-a44e-5fcc-9500-b78acb0c13bc.png)

　そしてコピーした地図情報を下記にペーストして入力します．

```python
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での表示用に加工します．

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

```python
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
```
それでは，設定した関心域を確認します．

```python
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
```
出力
![Screenshot from 2020-09-22 11-05-54.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/451824b6-05f2-bf58-ee17-a8d8ab397324.png)


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

[Sentinel-1 SAR GRD: C-band Synthetic Aperture Radar Ground Range Detected, log scaling](https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S1_GRD#description)
[Sentinel-2 MSI: MultiSpectral Instrument, Level-1C](https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2)

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

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

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

```python
region=ee.Geometry.Rectangle(area['longtitude'].min(),area['latitude'].min(), area['longtitude'].max(), area['latitude'].max())
```
　次に，取得する情報のパラメータを設定します．今回は，取得する画像の期間と，取得した画像の保存先を指定しています．
　今回は，異なる日の衛星画像を取得するため，まずは最近の衛星画像を取得します．ここでは，2020年9月20日としています．

```python
#期日を指定
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の画像を用いませんが，参考のために取得します．

```python
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の関心域の画像を取得します．

```python
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)
```
```python
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の観測画像からです．

```python
# 時系列で可視化
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()
```
![Screenshot from 2020-09-22 11-13-15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/a64fa5e9-dc74-a73b-b6c5-bd4de98df27f.png)


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

```python
# 時系列で可視化
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()
```
![Screenshot from 2020-09-22 11-39-26.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/42327912-09c1-6724-a036-bf415f4c0f14.png)
　これより，同一の観測日の2つの観測画像が取得されています．同じエリアで2つの観測画像に別れるのは，Sentinel-1の衛星画像の標準サイズが100km☓100kmであり，今回の関心域がそれを横断していたため，と思われます．
　そこで，2つに分割した衛星画像をまずは一つにまとめます．

```python
# データの読み込み
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])
```

```python
# データの読み込み
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つの画像を結合します．

```python
image_BG = image_1 + image_2
plt.imshow(image_BG, cmap='gray')
```
![Screenshot from 2020-09-22 11-43-36.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/791e41b0-ecfc-4877-b3ed-b0fb8d18c1e2.png)

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

```python
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ピクセル内の平滑化処理を行いました．

```python
image_BG_lee = lee_filter(image_BG, 3)

plt.imshow(image_BG_lee,cmap='gray')
```
![Screenshot from 2020-09-22 11-48-40.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/03a22889-0acd-8702-0344-f82f5b7aa2af.png)

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

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

```python
# 　観測期間

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
```
```python
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)
```
　では，取得した画像を表示してみます．

```python
# 時系列で可視化
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()
```
![Screenshot from 2020-09-22 11-51-34.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/38aba076-b112-60c5-0b9a-71726cfc9610.png)

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

```python
image_R_lee = lee_filter(image_BG, 3)
plt.imshow(image_R_lee,cmap='gray')
```
　次に，最近の衛星画像と過去の衛星画像を合成した**RGB画像**を作成します．今回は，**過去画像を赤色**に，**最近の画像を青色および緑色**としています．

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

```python
#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()
```
![Screenshot from 2020-09-22 12-00-05.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/7e57f3d9-5e40-aa0b-369a-1793f5e710ef.png)

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

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

```python
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()
```
![Screenshot from 2020-09-22 13-03-09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/b11d52a9-3c0b-ef51-f245-ca81da9e3a7a.png)
　このままでは全体に白い画像となっており，変化している青や赤が見づらくかんじます．ここで，下記を実行しコントラストを変更します．ここはお好みですね．

```python
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()
```
![Screenshot from 2020-09-22 13-04-41.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/4658f102-4ce0-4cf4-36b5-1b426c5e9435.png)
　先程よりも青と赤の変化が見やすくなっていすね．
　次に，この画像に衛星画像の撮影日およびクレジットを記載します．ここで使用するフォントはネットにてフリーで公開されているものを利用します．これもお好みですね．

```python
#フォントファイルのダウンロードと設定
!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"
```
　では，画像に文字を書き込みます．

```python
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()
```
![Screenshot from 2020-09-22 13-07-23.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/5d3ee9dd-949f-f678-73e1-1131af9da233.png)
　これで画像に必要な情報を記載することができました．文字の位置やサイズ，色などは上記のコードを変更して，好みに調整してください．
　次に，ここで取得した画像をローカルのPCで**GISアプリ**(例えば，**Google Earth**)で閲覧できるよう**kmzファイル**を作成します．

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

[人工衛星画像の無料環境なプラットフォーム「Tellus」にてデータ解析． -Google Earth編-](https://qiita.com/nigo1973/items/4bee0280b2569709d836)

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

```python
!pip install simplekml

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

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

north = area["latitude"].max()
south = area["latitude"].min()
east = area["longtitude"].max()
west = area["longtitude"].min()
```
　次に，上記の関心域のkmlファイル(緯度・経度情報）を作成します．

```python
#地理情報を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ファイル**を作成します．

```python
#地理情報と画像とをあわせた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**が表示されますので，ファイルをダウンロードします．

![Screenshot from 2020-09-22 13-17-24.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/03d72256-8ded-824f-e1e5-37463cbb42e5.png)

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

![Screenshot from 2020-09-22 13-19-43.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187602/b63b2cf3-ef2b-a1ce-d0f2-5152ab46a431.png)
　画像から青や赤くなっている場所を確認し，kmzファイルの選択を外すと該当する場所がどこか確認できます．**お台場の青い部分はコンテナの置き場所？**であることがわかりました．それ以外には，**東京ディズニーランドの駐車場で一部が青く**なっており，ここは計画されていたディズニーシーの**新しいアトラクションの建設場所**と推測されます．（Google Earthの航空写真は時間スライダーがありますので，該当箇所の撮影日や過去のものなど確かめてみてください．）


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

　[災害時の人工衛星活用ガイドブック水害版・衛星基礎編](https://www.mlit.go.jp/common/001227724.pdf)

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

##参考記事
[衛星画像より駐車場の利用状況を調べてみた．](https://qiita.com/nigo1973/items/c3be4aca9d882cb5a0e1)
[衛星画像より東京湾の船数のトレンドを求めてみた．](https://qiita.com/nigo1973/items/eeface437e7be513605a)
[無料で最新の衛星画像を入手する方法．](https://qiita.com/nigo1973/items/9bb6a11caac8e3e1e850)


[衛星データのキホン~分かること、種類、頻度、解像度、活用事例
~](https://sorabatake.jp/279/)
[人工衛星から人は見える？~衛星別、地上分解能・地方時まとめ~](https://sorabatake.jp/441/)

[Sentinel-1 SAR GRD: C-band Synthetic Aperture Radar Ground Range Detected, log scaling](https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S1_GRD#description)
[Sentinel-2 MSI: MultiSpectral Instrument, Level-1C](https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2)

　[災害時の人工衛星活用ガイドブック水害版・衛星基礎編](https://www.mlit.go.jp/common/001227724.pdf)