衛星画像より駐車場の利用状況を調べてみた．

【更新】作成した東京TDLの駐車場利用状況のトレンドを2020年8月31日時点に更新しました。

概要

　これまでSentinel衛星画像のハブであるCopernicus hubより観測画像を取得し，その画像の解析を行っていました．

無料で最新の衛星画像を入手する方法．
人工衛星（Sentinel-2)の観測画像をAPIを使って自動取得してみた．

観測画像の過去画像との比較や，その時系列変化のGifアニメなどは作りましたが，APIで取得できる過去画像が過去１年間程度であり，より長期なトレンドをみたいときは，別途Copernicusの管理者に連絡して用意してもらう必要がありました．最新の観測画像は取得できるのでこれでよいかと考えていましたが，知り合いから「長期のトレンドが知りたいんだけど」と相談があり，マニュアルでの作業を考えていました．
　そんなとき，以下の記事を知りみてみると，Google Earth Engineを使えばSentinel衛星の過去の観測画像がオンラインで取得し解析できることを知りました．ありがとうございます．

　【続編】Google Earth EngineとGoogle Colabによる人工衛星画像解析　〜無料で始める衛星画像解析（実践編）〜

Google Earth Engine（以下，GEE)は知っていましたが，GEEの環境でしか使えないものであり，言語がPythonではないので敬遠していました．しかし，上記の記事からGoogle colaboratoryで利用でき，これまでのようにSentinelの観測画像を取得してから処理しなくてもよいため，これは便利じゃない！と調べたところ比較的簡単だったので紹介します．

解析例は，以下の東京ディズニーランドの駐車場の利用状況のトレンドです．

このグラフの説明も含めてこの記事で紹介します．
また，記事で紹介したコードは Githubに置きましたので，ご自身で試したい方はこちらをご利用ください．Google colaboratoryを用いていますので，ネットワーク環境さえあればPCによらずご利用できます．

1.はじめに

　衛星画像から駐車場の利用状況を見積もり方法は，例えば以下の記事が参考になるかと思います．

いつ空いてるの!? 無料衛星データでディズニーランドの混雑予想チャレンジ（前編）

これより，衛星画像でも電波センサの衛星であるSAR衛星の観測画像が使われてるのがわかります．それは，光学観測画像は車を確認するためには50cm以下の高分解能の画像が必要であり，かつ対象地域が雲のない晴天であることが必要になります．アジア地域の雲がかかっている被雲率は約30%と言われており，場所によっては何日も雲がかかって衛星画像ではみえないことがあります．
　それと比べてSAR画像は，車かトラックかなどの車種の特定は難しいのですが，駐車場に多数の車が駐車されているかどうかであれば定性的には評価することができ，かつ雲の影響を受けないため安定して評価することができます．
　今回は航空写真から車の台数分布を推定したときとはことなり，ある日時での車の利用状況を基準として，利用状況がどのように推移しているのか相対的なトレンドを求めます．そのためSAR画像であっても分析が可能です．
　とはいえ，同じSAR画像でも分解能に差があり，できれば高分解能なSAR画像を利用したいのですが，無料でありかつ定期的に観測しているSentinel-1衛星の観測画像を用います．
　SAR画像による駐車場の車などの構造物の密集具合は画像からある程度推量することができます．絵心がないので，関連する記事を参照させてもらいます．

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

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

ディズニーランドでの施設がよくわかりますね．また，ディズニーランドに行かれた方がご存じだと思いますが，周囲に駐車場があり，そこは広く平坦なため暗く写っているのがわかるかと思います．この駐車場に車が駐車されると明るく映るため，駐車場エリアの明るさを評価することで，相対的に駐車場の利用状況を推量することが今回の評価方法になります．
　今回利用する人工衛星は欧州宇宙機関が開発運用しているSentinel-1と2になります．それぞれの詳細については，以下の記事をご参考下さい．

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

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

Sentinel-1 Observation Scenario

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

先程紹介しましたSentinel-1は，A衛星が午前6時に，B衛星が午後6時が観測時刻になります．そのため，A衛星を用いると早朝の駐車状況を評価することになります．この時間から駐車場で待機している方は稀だと思いますので，B衛星の観測画像のみを今回用いました．また，今回の駐車場評価以外でSentinel−1衛星を利用し時系列的な評価をされる方は，A衛星かもしくはB衛星のどちらかの観測画像のみを用いることをお薦めします．それは，衛星によって観測方向が異なるため，建物などの高さのある構造物を評価する場合，電波の照射方向が異なるため，陰影が変わるため，単純な前後比較では時間成分以外の差が生じます．
　次に，GEEより衛星画像の取得，およびその評価について紹介します．

２．衛星画像の取得と評価

2.1　環境準備（構築）．

　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":"1127.16 m","drawtype":"rectangle","area":"11.71 ha"},"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[139.87487004138532,35.63225823758155],[139.87487004138532,35.6360949926148],[139.87790093757215,35.6360949926148],[139.87790093757215,35.63225823758155],[139.87487004138532,35.63225823758155]]]}}]}

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

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

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=15)

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

次に，取得する情報のパラメータを設定します．今回は，取得する画像の期間と，取得した画像の保存先を指定しています．

# 期間を指定
from_date='2019-01-01'
to_date='2020-08-31'

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

では，Sentinel-1と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()

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', 20).map(cloudMasking).select(['B4','B3','B2'])

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

上記の設定はこちらの記事を参考にしました．
6. PythonによるローカルからのGEE実行
上記をベースに，Sentinel-1の観測画像の設定や，Sentinel-2では雲画像の削除を追加しています．

ここで，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　衛星画像の表示と時系列データの評価

取得した衛星画像を表示し確認します．
衛星画像は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であり，これでは車の車種などの識別は難しいですね．今回はSentinel-2の観測画像を評価には利用しませんが，Sentinel-1の観測画像の取得時期の駐車場がどのような状況だったのか，概要はわかります．
東京ディズニーランドは数年前以上行っていないため今の状況はわかっていませんが，路上駐車場に立体駐車場か何かの建物（左上）が建てられており，観測画像の取得時の2019年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()
  j+=1# 画像のプロット位置をシフトさせ配置
  plt.subplot(v,5,j)
  plt.imshow(arr[0], cmap='gray')
  plt.title(files[i][33:41])# ファイル名から日付を取得
  plt.tight_layout()

　SAR衛星の観測画像は白黒画像であり，これだけでは何かよくわかりませんね．先程確認したSentinel-2の観測画像が，この画像の理解に役立ちます．

　次に，Sentinel-1の各観測画像の信号（明るさ）の積算値を求め，そのトレンドを取得します．

# 当該エリアの反射強度の合計値を時系列グラフ化
sum_signal = []
label_signal = []
for i in range(len(files)):
  # 画像を1シーンずつ取得して可視化
  with rasterio.open(s1_path + files[i]) as src:
      arr = src.read()
  sum_signal.append(arr.sum())
  label_signal.append(files[i][33:41])

# 可視化
fig,ax = plt.subplots(figsize=(15,6))
plt.plot(sum_signal, marker='o')
ax.set_xticks(np.arange(0,len(files)))
ax.set_xticklabels(label_signal, rotation=90)
plt.title('Trend in parking lot usage at TDL.')
plt.xlabel('date')
plt.show()

　横軸が観測日，縦軸が観測画像の信号の積算値（画像全体の信号の積算値）になります．縦軸の値の物理的な意味はありませんが，2019年1月18日の値を基準に評価します．
　2019年1月から3月にかけて値が増加しています．Sentinel-2の観測画像をみると，この期間は駐車場エリアに建物（立体駐車場？）が建設されており，この構造物によってSAR画像の信号が増加（対象エリアの建物の密集度が増加）したためと思われます．その後，2020年1月末までは多少の上下はするもののおおきな変化は見られませんが，2020年2月になると急激に値が低下しそれを維持しています．これは，新型コロナウィルスの感染拡大の影響と考えます．
　実際，東京ディズニーランドは2020年2月29日から運営を停止しており，7月1日より運営が再開されました．運営の再開に伴い，値の増加がみられましたが，8月にはいってまた下がっています．これは感染が再拡大している影響があると思います．夏休みになって多くの方が来場されることが期待されていたと思いますが，難しい状況がつづいていると思います．感染拡大の予防に努めつつ，テーマパークが再び楽しい場所に戻ることを祈っています．

3. 最後に

　Google が提供するGoogle Earth Engineを用いて，衛星画像の取得および解析例を紹介しました．
　今回は，ディズニーランドの駐車場の利用状況を対象としましたが，対象地域の構造物が増えるとSAR信号が増加しますので，駐車場以外の都市域の開発状況のトレンドなどにも応用できると思います．
　これを機会に衛星画像に関心を持つ方が増えると嬉しいです．
　ご意見や質問などありましたら，お気軽にコメントください．嬉しいので．

参考記事

無料で最新の衛星画像を入手する方法．
人工衛星（Sentinel-2)の観測画像をAPIを使って自動取得してみた．
PyTorchによる人工衛星画像から車の推定分布地図を作成してみる．

　【続編】Google Earth EngineとGoogle Colabによる人工衛星画像解析　〜無料で始める衛星画像解析（実践編）〜
　Google Earth EngineとGoogle Colabによる人工衛星画像解析　〜無料で始める衛星画像解析（入門編）〜
6. PythonによるローカルからのGEE実行

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

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