概要
前回のGoogle Earth Engineより入手した衛星画像より東京ディズニーランドの駐車場の車数の推移を求めました．
今回は，前回と同様にGoogle Earth Engineより衛星画像を入手し，東京湾の船舶数の推移を求めてみます．
こちらは，東京湾入り口付近の人工衛星の撮影画像にある船舶をカウントし，そのトレンドを求めたものです．
このグラフの説明も含めてこの記事で紹介します．
また，記事で紹介したコードは Githubに置きましたので，ご自身で試したい方はこちらをご利用ください．Google colaboratoryを用いていますので，ネットワーク環境さえあればPCによらずご利用できます．
#1.はじめに
前回は東京ディズニーランドの駐車場の車の数を，衛星画像のなかの電波センサ（SAR)画像を用いてトレンドを求めました．
SAR画像は車かトラックかなどの車種の特定は難しいのですが，駐車場に多数の車が駐車されているかどうかであれば定性的には評価することができ，かつ雲の影響を受けないため安定して評価することができます．
そのため，今回もSAR画像を用いて，東京湾の入口付近の衛星画像から船舶数を求めてみます．これは，新型コロナウィルスの影響で船舶数がどのように変わったのか，変わっていないのか，その程度が評価できるのではないか，と思ったからです．
SAR画像にて船舶がどのように見られるのか，こちらのサイトが参考になるかと思います．
陸域観測技術衛星2号「だいち2号」（ALOS-2）初期機能確認期間中の観測画像について
（その5・船舶編）
だいち2号のこの撮像画像は，分解能3mの観測モードでの画像であり，大型の船舶であればその形状もわかるかと思います．ここで重要なポイントは，船舶は明るく写り，周囲の海は暗く写る，ということです．
SAR衛星は電波を地球に照射し，その反射信号の強度を画像化しています．このとき，電波を地表面に対して斜め方向から放射するため，地表面が平坦だと多くのものは電波の入射方向，つまり放射した電波を受信する衛星とは逆方向に反射します．そのため，対象の地表面が平坦であれば受信する信号が弱いため暗い画像となります．湖や海などの水面は極めめて平坦なためより暗く写ります．
一方，建物やビル，車や船などの構造物がある場合は，対象物が凹凸しているため，入射方向と同じ方向へ反射する電波が多くなり，その為明るく写ります．
この明るく写る箇所を船とし,衛星画像内の船舶数を数えてみます．
#２．衛星画像の取得と評価
##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":"35937.06 m","drawtype":"rectangle","area":"12713.15 ha"},"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[139.71197873353958,35.33475055335104],[139.71197873353958,35.452996930987666],[139.81840878725052,35.452996930987666],[139.81840878725052,35.33475055335104],[139.71197873353958,35.33475055335104]]]}}]}
この地理情報をGEEへの入力フォーマットに，そしてFoliumでの表示用に加工します．
#今後使用する任意のファイル名をセットする． 例えば，地域の名前など．
object_name = 'tokyobay3'
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=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-1の観測画像は2014年10月3日から，のデータが利用できます．
では，GEEよりSentinel-1の画像を取得し，Google colaboratoryに保存します．
最初に，GEEに設定する地理情報のフォーマットを準備します．
region=ee.Geometry.Polygon(AREA)
次に，取得する情報のパラメータを設定します．今回は，取得する画像の期間と，取得した画像の保存先を指定しています．
# 期間を指定
from_date='2019-01-01'
to_date='2020-08-31'
# 保存するフォルダ名
dir_name_s1 = 'GEE_Sentinel1_' + object_name
では，Sentinel-1の画像条件の設定です．
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'])
imageList_s1 = Sentinel1.toList(300)
ここで，Sentinel-1の観測画像は，B衛星の午後6時の観測画像のみを用いるため，'orbitProperties_pass'で 'ASCENDING'を選択しています．”Descending”にすると午前6時の観測画像になります．
では，上記の衛星画像の取得条件で，Sentinel-1の関心域の画像を取得します．
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)
##2.4 衛星画像の表示と信号強度の評価
取得した衛星画像を表示し確認します．
衛星画像はGoogle driveのmy driveに設定したディレクトリ（フォルダ）に保存されています．それを呼び出し表示します．
# 時系列で可視化
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衛星の観測画像は白黒画像であり，この画像内の明るい部分が船舶と思われます．次に，画像内の明るい部分を船舶としてその数をカウントします．
# データの読み込み
n = 2
with rasterio.open(s1_path + files[n]) as src:
arr = src.read()
print(files[n][33:41])
# 可視化
plt.imshow(arr[0], cmap='gray')
衛星画像を拡大してみると，明るい点が多数あるのが確認されました．では，この画像の強度分布を求めます．
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.hist(arr[0], bins=50)
fig.show()
信号強度は0以下が多いのですが，いくつか0以上の信号も観測されています．この条件より，0以上とそれ以下とで2値化処理をまずは行います．
##2.5 衛星画像より船舶数を算出する．
画像内の物体のカウント方法について，以下のサイトを参考にさせていただきました．
この方法は，任意の信号強度で2値化を行い，広がりがある場合は中心位置を求め，その中心位置が同じ物体かそうでないのかを任意のウィンドウにあるかないかを調べ，物体の数を数えます．
では，具体的にその方法をご紹介します．
まず，サンプル画像は取得した画像内の最新の撮像画像としました．
# データの読み込み
n = len(files) -1
with rasterio.open(s1_path + files[n]) as src:
arr = src.read()
print(files[n][33:41])
# 可視化
plt.imshow(arr[0], cmap='gray')
多数の明るい点があるのがわかるかと思います．次に，二値化処理を行います．ここでは，先の信号強度分布を参考に0以上とそれ以下にて二値化処理を行います．また，求めた2値化画像から中心位置のマップを作成します．
#二値化処理
import cv2
binimg = (arr[0]>0) #0をしきい値とする．
plt.imshow(binimg)
plt.colorbar()
plt.show()
#距離マップ(distance map)を計算する
binimg = binimg.astype(np.uint8)
distmap = cv2.distanceTransform(binimg,1,3)
plt.imshow(distmap)
plt.colorbar()
plt.show()
多数の船が撮像された部分を拡大します．
plt.imshow(distmap[0:200,0:200])
plt.colorbar()
plt.show()
船舶と思われる信号がある方向に伸びているのがわかります．
次に，この画像から信号強度が最大値をとる分布を取得します．このとき，船のサイズを30ピクセルとしました．つまり，30ピクセル内に船は一艘のみである，と定義しています．
#距離マップの極大値をもつ場所を計算する
out = distmap*0
ksize=30 #30ピクセル以内の最大値
for x in range(ksize,distmap.shape[0]-ksize*2):
for y in range(ksize,distmap.shape[1]-ksize*2):
if distmap[x,y]>0 and distmap[x,y]==np.max(distmap[x-ksize:x+ksize,y-ksize:y+ksize]):
out[x,y]=1
plt.imshow(distmap[0:200,0:200])
plt.colorbar()
plt.show()
次に，船と思われる位置の近くの画素で高い信号強度となった場合誤検知する可能性があるため，ある範囲内にある場合は無視するように処理します．
#膨張させ、輪郭検出し、数える
out = cv2.dilate(out,(10,10)) #10*10に膨張させ，同じ範囲に入っていたらカウントしない
contours, _ = cv2.findContours(out.astype(np.uint8), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
arr=[]
for i in contours:
x_=0
y_=0
for j in i:
x_ += j[0][0]
y_ += j[0][1]
arr.append([x_/len(i), y_/len(i)])
arr = np.array(arr)
plt.imshow(distmap[0:200,0:200])
plt.colorbar()
plt.show()
print("船の数: ", len(arr))
print(arr)
ここでは，船舶の数，およびその位置が出力されました．
この方法を用いて，取得した衛星画像より船舶数を求め，その推移をグラフにします．
# 当該エリアの推定された船の数を時系列グラフ化
# 時系列で可視化
s1_path = '/content/drive/My Drive/' + dir_name_s1 + '/'
files =os.listdir(s1_path)
files.sort()
sum_ship = []
label_signal = []
for i in range(len(files)):
label_signal.append(files[i][33:41])
# 画像を1シーンずつ取得して可視化
with rasterio.open(s1_path + files[i]) as src:
arr = src.read()
#二値化処理
binimg = (arr[0]>0)
#距離マップ(distance map)を計算する
binimg = binimg.astype(np.uint8)
distmap = cv2.distanceTransform(binimg,1,3)
out = distmap*0
ksize=30 #30ピクセル以内の最大値
for x in range(ksize,distmap.shape[0]-ksize*2):
for y in range(ksize,distmap.shape[1]-ksize*2):
if distmap[x,y]>0 and distmap[x,y]==np.max(distmap[x-ksize:x+ksize,y-ksize:y+ksize]):
out[x,y]=1
out = cv2.dilate(out,(10,10)) #10*10に膨張させ，同じ範囲に入っていたらカウントしない
contours, _ = cv2.findContours(out.astype(np.uint8), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
num=[]
for i in contours:
x_=0
y_=0
for j in i:
x_ += j[0][0]
y_ += j[0][1]
num.append([x_/len(i), y_/len(i)])
num = np.array(num)
sum_ship.append(len(num))
# 可視化
fig,ax = plt.subplots(figsize=(15,6))
plt.plot(sum_ship, marker='o')
ax.set_xticks(np.arange(0,len(files)))
ax.set_xticklabels(label_signal, rotation=90)
plt.show()
横軸が観測日，縦軸が観測画像より算出した船舶数になります．
2019年3月付近に極大値をとりますが，それ以降は平均30隻程度になります．2020年3月以降のコロナ禍では船による移動が制限もしくは減衰しているのかと思いましたが，これを見る限り東京ティズニーランドの駐車数のように，極端に減少しているようには見えません．
実際，東京湾の船舶の航行数がどうなっているのかインターネットにて調べてみましたが，クルーズ船は運行を控えている記事を見つけましたが，その他の多くの船が減少しているものはありませんでした．コロナ禍でも物流の安定のために，みなさん安全に留意しつつ運行されていたのかと思います．ありがとうございます．
東京湾以外の船舶の航行が多い箇所も同様の方法で調べてみましたが，今回と同様にコロナ禍であっても極端に減少するなどの傾向はなく，船舶を運行されている方々のご努力によって，安定した物流が維持されているのだと実感しました．
#3. 最後に
Google が提供するGoogle Earth Engineを用いて，衛星画像の取得および解析例として東京湾の船舶数の求め方を紹介しました．
これを機会に衛星画像に関心を持つ方が増えると嬉しいです．
ご意見や質問などありましたら，お気軽にコメントください．嬉しいので．
##参考記事
無料で最新の衛星画像を入手する方法．
人工衛星（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