2022年8月31日以降、Tellus OSでのデータの閲覧方法など使い方が一部変更になっております。新しいTellus OSの基本操作は以下のリンクをご参照ください。
https://www.tellusxdp.com/ja/howtouse/tellus_os/start_tellus_os.html

1.衛星データで遊園地の混雑予測（概要説明）

本記事は「いつ空いてるの!? 無料衛星データでディズニーランドの混雑予想チャレンジ」(https://sorabatake.jp/349/)の後編となります！

前半ではディズニーの混雑予想をするための方法を列挙しました。今回は、そのうちの一つ、衛星画像から駐車場の混雑具合を観測するという手法で、実際にうまく観測できるかどうかを調べて見ました！

※本記事は宙畑メンバーが気になったヒト・モノ・コトを衛星画像から探す不定期連載「宇宙データ使ってみた-Space Data Utilization-」です。まだまだ修行中の身のため至らない点があるかと思いますがご容赦・アドバイスいただけますと幸いです！

2.解析に使う衛星データ

- 衛星の種類

今回は、雨や曇りの時でも地表が見えるSARデータを使っています。

中でもヨーロッパの衛星Sentinel-1を使いました。何といってもオープンソースデータが豊富なところ、素晴らしいです。ヨーロッパの政府衛星データ Copernicusのサイトへは下記のリンクから飛べます。

ヨーロッパの政府衛星データ Copernicus
http://www.copernicus.eu/

今回は、すでにユーザー登録をしている状態での解析になります。

ユーザー登録の仕方はこちらをご覧ください。

- 偏波

宙畑の他の記事でもご紹介している通り、SAR画像には「偏波」という項目があります。衛星から電波を発する時と受信する時の電波が進行する面を表しています。

EO browserで検索してみると、Sentinel-1ではVVとVHという偏波が多いことが分かります。

それぞれの画像をEO browser上で確認してみると以下のような感じになりました。

実際に画像を見比べてみて、より見やすいように思えたVVを採用しました。

- 軌道

人工衛星には、南から北へ向かうアセンディングと北から南へ向かうディセンディングの2つ軌道があり、衛星画像にはどちらの軌道で撮影した情報かが付加されています。

アセンディングとディセンディングは時系列で画像を比較する場合、撮影条件が異なってくるため、なるべく揃えた方が良い項目と考えられます。

Sentinel-1の場合、太陽同期軌道という撮影地点の時刻がほぼ同じになる軌道を飛んでおり、アセンディングはおよそ18時、ディセンディングはおよそ6時に舞浜周辺を撮影しています。

今回は、ディズニーが営業している時間帯に合わせてアセンディングで観測したデータを取得しました。

因みに回帰日数（観測頻度）は12日で、2機稼働しているので実質6日置きとなります。週別の解析は難しいですが、月別での駐車場状況はみれそうですね…

なお、本章で紹介した衛星データのスペックに関する説明は以下の記事で解説していますので、気になる方はぜひご覧ください。

3.解析の仕方

今回は以下の記事のコードを参照しています。

- 衛星データを見てみる

まず、必要なものをインストール＆インポートします。

!pip install sentinelsat
!pip install rasterio
!apt install gdal-bin python-gdal python3-gdal 
!apt install python3-rtree 
!pip install git+git://github.com/geopandas/geopandas.git
!pip install descartes 
!pip install shapely
!pip install six
!pip install pyproj
!pip install descartes
!pip install geopandas
!pip install folium
!pip install tifffile
!pip install numexpr

import os
import numpy as np

from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt 
import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

from shapely.geometry import MultiPolygon, Polygon
import rasterio as rio
from rasterio.plot import show
import rasterio.mask
import fiona
import folium 
import zipfile
import glob
import gdal
import tifffile
from PIL import Image,ImageDraw,ImageFont
import sys
import shutil

次に、キーン州立大学が公開しているツールを利用し、取得したい領域の座標を入手します。

from IPython.display import HTML
HTML(r'＜iframe width="960" height="480" src="https://www.keene.edu/campus/maps/tool/" frameborder="0"＞＜/iframe＞')
# 実際に実行するときは＜を全角から半角に変更して使用してください。

右クリックで4点選択した後に“Close Shape”を押すと、四角形のポリゴンを入手できるようです。例えばこんな風に。

そして右下の”coordinates”から以下を下記にペーストしてください。これで、取得する範囲を定義できました。(バージョンアップに伴い、上記の経度緯度と表示のされ方が違うのですが、取得する範囲は同じです。)

# ディズニーのエリア
AREA=[
      [
        -220.1284218,
        35.6412547
      ],
      [
        -220.1285934,
        35.621861
      ],
      [
        -220.1038742,
        35.6217214
      ],
      [
        -220.1042175,
        35.6412547
      ],
      [
        -220.1284218,
        35.6412547
      ]
    ]

定義した領域のポリゴン情報を生成し、ファイル名をつけます。

# 経度の値を戻すために360度足す
for i in range(len(AREA)):
    AREA[i][0] = AREA[i][0] +360
from geojson import Polygon
m=Polygon([AREA]) 

#ファイル名を定義. 好きな名称に設定してください. 
object_name ='XXXX'

import json
with open(str(object_name) +'.geojson', 'w') as f:
    json.dump(m, f)
footprint_geojson = geojson_to_wkt(read_geojson(str(object_name) +'.geojson'))

# ポリゴンがきちんと作れているかの確認

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

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

画像に日付を入れたい場合は以下のコードを使います。

# フォントファイルの作成（任意）
# 同じディレクトリ内に該当のファイルをダウンロードしておく
# https://osdn.net/projects/mplus-fonts/downloads/62344/mplus-TESTFLIGHT-063a.tar.xz/
!xz -dc mplus-TESTFLIGHT-*.tar.xz | tar xf -
fontfile = "./mplus-TESTFLIGHT-063a/mplus-1c-bold.ttf"

Sentinel Hubからデータを取得するために、事前にWebサイトで登録したSentinel Hubのユーザ情報を入力します。

Sentinel Hub
https://scihub.copernicus.eu/dhus/#/self-registration

Sentinel-1のAPIリファレンスも載せときます。(https://scihub.copernicus.eu/userguide/OpenSearchAPI)

#Sentinel Hubhttps://scihub.copernicus.eu/dhus/#/self-registrationでアカウントを作成し、アカウント情報を以下へ入れる。
user = '自身のユーザ名を貼り付けてください' 
password = ' 自身のパスワードを貼り付けてください' 
api = SentinelAPI(user, password, 'https://scihub.copernicus.eu/dhus')

画像を取ってくる関数は以下の様に定義します。

def Sentinel1_get():
    #Sentinel-1のデータ取得条件
    products = api.query(footprint_geojson,
                 date =  (Begin_date, End_date), #取得希望期間の入力
                 platformname = 'Sentinel-1',
                 producttype = 'GRD', #SLC, OCN, RAW,
                 Polarisationmode = 'VV'
#                  productlevel = 'L1',
#                  Mode = 'IW',
#                  Polarisation = 'VH', #'VV \n",
#                  Product =  'GRD'
                 )  

#     print(products) 
    products_gdf = api.to_geodataframe(products) 
#     print(len(products))
#     print(products_gdf["orbitdirection"])
  
    #Ascending観測のみを抽出
    product_ged_descending = products_gdf[products_gdf['orbitdirection'] == 'ASCENDING']
    print(product_ged_descending)


    uuid = product_ged_descending.iloc[0]["uuid"]
    print(uuid)
    product_title = product_ged_descending.iloc[0]["title"]

    #データの観測日を確認．
    date = product_ged_descending.iloc[0]["ingestiondate"].strftime('%Y-%m-%d')
    print(date)
    #Sentinel-1 data download
    api.download(uuid)
    file_name = str(product_title) +'.zip'

    # def desneylandy(filename):
    with zipfile.ZipFile(file_name) as zf:
        zf.extractall()

    path = str(product_title) + '.SAFE/measurement'
    files = os.listdir(path)

    path_c = glob.glob(str(product_title) + '.SAFE/measurement/*[v][v]*.tiff')
        #path_b = str(product_title) + '.SAFE/measurement/' + str(files[1])

    path_b = path_c[0]
    b = rio.open(str(path_b))
    print(b)

    gdal.Warp(path_b, path_b, dstSRS="EPSG:4326")

    nReserve_geo = gpd.read_file(str(object_name) +'.geojson')

            #src = rio.open(r\"RGB.tiff\")
    epsg = 'epsg:4326'
    nReserve_proj = nReserve_geo.to_crs({'init': epsg})

     #マスクTiffファイルの一時置き場
    os.makedirs('./Image_tiff', exist_ok=True)  
    
    with rio.open(path_b) as src:
        out_image, out_transform = rio.mask.mask(src, nReserve_proj.geometry,crop=True)
        out_meta = src.meta.copy()
        out_meta.update({"driver": "GTiff",
                         "height": out_image.shape[1],
                         "width": out_image.shape[2],
                         "transform": out_transform})
    with rio.open('./Image_tiff/Masked_' +str(object_name) +'.tif', "w", **out_meta) as dest:
    # with rasterio.open('home/jovyan/work/Image_tiff/Masked_' +str(object_name) +'.tif', "w", **out_meta) as dest:
        dest.write(out_image)

    print('jpeg generation')

    scale = '-scale 0 250 0 60' # '-scale 0 250 0 30'
    options_list = [
            '-ot Byte',
            '-of JPEG',
            scale
          ] 
    options_string = " ".join(options_list)

        #ディレクトリの作成 
    os.makedirs('./Image_jpeg_'+str(object_name), exist_ok=True)

       
    
        #jpeg画像の保存．
    gdal.Translate('./Image_jpeg_'+str(object_name) +'/' + str(Begin_date) + 'Masked_' +str(uuid) +'.jpg',
                   './Image_tiff/Masked_' +str(object_name) +'.tif',
                   options=options_string)

            #画像への撮像日を記載
    img = Image.open('./Image_jpeg_'+str(object_name) +'/' + str(Begin_date) + 'Masked_' +str(uuid) +'.jpg')

    x = img.size[0]/100 #日付の記載位置の設定
    y = img.size[1]/100 #日付の記載位置の設定
    fs = img.size[0]/50 #日付のフォントサイズの設定
    fs1 = int(fs)

    obj_draw = ImageDraw.Draw(img)
    obj_font = ImageFont.truetype(fontfile, fs1)
    obj_draw.text((x, y), str(date), fill=(255),ont=obj_font)
    obj_draw.text((img.size[0]/2, img.size[1]-y - img.size[1]/20 ), 'produced from ESA remote sensing data', fill=(255),font=obj_font)
    img.save('./Image_jpeg_'+str(object_name) +'/' + str(Begin_date) + 'Masked_' +str(uuid) +'.jpg')
#     print(img)
            #不要なファイルの削除
    shutil.rmtree( str(product_title) + '.SAFE')
    os.remove(str(product_title) +'.zip')

    return

ところが、この領域の画像を実際に表示させてみると、

このように、地図で見た範囲とはずれている領域が取得されました。
衛星画像を何枚か見ていくと、画像によって、微妙に位置がずれていることが分かりました。
つまり、取得する衛星画像の位置情報にずれがあり、同じ緯度経度情報で取得すると撮影範囲に差が出てしまうようです。

今回は画像を見ながら手作業で緯度経度を調整しています。

# ディズニーのエリア（調整後）
AREA=[
      [
        -220.11445760000004, # 左上
        35.643927399999995
      ],
      [
        -220.1148008,#  左下
        35.623513
      ],
      [
        -220.09107740000002,#  右下
        35.622887
      ],
      [
        -220.09153500000002,#  右上
        35.64358839999999        
      ],
      [
        -220.11445760000004,#  左上
        35.643927399999995
      ]
    ]

調整した結果がこちら。
このように、ディズニーランドとシー全体の画像を取得できました。

ある期間の画像を取ってくる

sentinel-1の画像データは、現在から1年前の月までの画像であればオンラインで取得できます。例えば、今日が7月24日だとしたら、一年前の7月1日のデータまでは得られるようです。

長期アーカイブ機能も使えば1年前より古いデータが取得できるかもしれません。詳しい話は、こちらから確認してみてください。

# 取得する期間を設定(過去1年分の指定が可能)
Begin_date = '20190901'
End_date = '20191225'

以下のコードで複数月にまたがる画像を持ってきます。

#以下のコードにて月数をカウントする
d_year = int(End_date[2:4]) - int(Begin_date[2:4])
d_year = d_year*12 + int(End_date[4:6]) - int(Begin_date[4:6])


for i in range(d_year + 1):
    if i < 1:
        m = int(Begin_date[4:6])
    else:
        m = int(Begin_date[4:6]) +1
    #print(m)
    if m <10:
        Begin_date = Begin_date[:4] +'0'+ str(m) + Begin_date[6:]
    elif m <13:
        Begin_date = Begin_date[:4] + str(m) + Begin_date[6:]
    else:
        y = int(Begin_date[2:4]) +1
        Begin_date = Begin_date[:2] + str(y) + '01' + Begin_date[6:]
    End_date = Begin_date[:6] +'28' 
    
    Sentinel1_get()

駐車場の形に切り抜く

続いて、取得してきた画像から駐車場の部分だけを切り抜きます。

まずは駐車場の位置を確認します。

必要なものをインポートします。

私は新規ファイルで作成したので、インポートしましたが、すでにインポートしている場合は必要ありません。

import glob
from PIL import Image
import numpy as np
import os, json, requests, math, glob
from skimage import io
from io import BytesIO
import matplotlib.pyplot as plt
from skimage.draw import polygon
%matplotlib inline
import cv2

切り抜く前に、取得した画像をリサイズしておきましょう。

後々の処理がしやすくなります。

新しくディレクトリを作成しリサイズした画像を保存しています。

path3='./resize_image_sisaku/'
os.makedirs(path3, exist_ok=True)

# リサイズ 
for f in glob.glob("./Image_jpeg_disneysisaku/*.jpg"):
    img = Image.open(f)
    img_resize = img.resize((256, 256))
#     title, ext = os.path.splitext(f)
#     img_resize.save(title + '_resize' + ext)
    imgdir3 = os.path.dirname(path3)
    imgname3 = os.path.basename(f)
    newfname3 = imgdir3 + "/resize_" + imgname3
    print(newfname3)
   
    img_resize.save( newfname3)

改めて画像を眺めて見ます。すると、画像によって駐車場のエリアの明るさ（白さ）が異なることが分かります。

コンクリートなどの平らな地面は、電波が全反射するため黒く、複雑な人工物である車などは電波が強く反射して見えます。

つまり、切り抜いた画像から黒のピクセル数をカウントすると駐車場の空き具合が分かり、白のピクセル数をカウントすると駐車場の混み具合が分かるということになります。

切り抜き方はこちらから。

# Tellusの開発環境のトークンを貼り付ける
# トークンの取得は以下の記事を参考に
# https://www.tellusxdp.com/ja/howtouse/dev/install_development_environment.html

TOKEN = "ここには自分のアカウントのトークンを貼り付ける"

# Tellus OSからgeojsonファイルを取ってくる
def  get_shape_geojson_by_name(name):
    URL = "https://api.tellusxdp.com/v1/shapes"

    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    
    r = requests.get(URL , headers=headers)
    print(r)
    shapes = r.json()

# Tellus OSからシェープの一覧を取ってくるときの一覧を表示
#     print(shapes)

    shape = next((shape for shape in shapes["items"] if shape["name"]==name))
    return json.loads(shape["data"])

def get_polygon_image(points, bbox, imgsize):
    def world2pixel(bbox, size, lon, lat):
        """
        経度緯度をピクセル座標に変換
        """
        dist = ((bbox[2] - bbox[0])/size[0], (bbox[3] - bbox[1])/size[1])
        pixel = (int((lon - bbox[0]) / dist[0]), int((bbox[3] - lat) / dist[1]))
     
        return pixel
    
    pixels = []
    for p in points:
        pixels.append(world2pixel(bbox, (imgsize["width"], imgsize["height"]), p[0], p[1]))
    pixels = np.array(pixels)

    poly = np.ones((imgsize["height"], imgsize["width"]), dtype=np.uint8)
    rr, cc = polygon(pixels[:, 1], pixels[:, 0], poly.shape)
    poly[rr, cc] = 0
    return poly

#   (左下経度, 左下緯度, 右上経度, 右上緯度)
bbox=(139.871369, 35.621890, 139.895782, 35.641084 )

shapename = "disneyland17" #自分で作った図形の名前に変更する※OS上に登録された10個以内のものまでしかページネーションの都合で呼べないので注意
shape = get_shape_geojson_by_name(shapename)
points = shape["features"][0]["geometry"]["coordinates"][0]

取得した画像の位置情報に誤差があると、位置情報に基づいて切り取った画像にも当然ズレが生じます…

抜き出す領域である図形の位置やバウンディングボックスを人の手で調整し、それでもダメならその画像は解析に加えないようにしましょう。私の場合もずれが大きく手作業による調整が難しいため、何枚か解析に加えられない画像がありました。

ピクセル数のカウント

そして実際に切り抜き、ピクセル数をカウントします。その後、新しいファイルを作って画像を保存します。

手作業で調整を行い、黒のピクセルを数える閾値を1以上60未満として計算しました。上記は黒ピクセルをカウントするコードです。
白の場合は、60以上256未満で計算しました。

# 閾値の設定
px_low = 1 
px_up  = 60
n = np.random.randint(1, 60, (256, 256))

# jpgのみを取得
imagedata = []
for i in glob.glob("./resize_image_sisaku/*.jpg"):
    jpg22 = cv2.imread( i, cv2.IMREAD_GRAYSCALE)
      
    pxsize2 = {"width": jpg22.shape[1], "height": jpg22.shape[0]}

    mask = get_polygon_image(points, bbox, pxsize2)
    gray_img_masked = jpg22 * (1-mask)
    
    io.imshow(gray_img_masked)
    
    binary = cv2.inRange(gray_img_masked, px_low, px_up)
    cnt = cv2.countNonZero(binary)
    print(cnt)
    
   
#     new file
    pathpark = './parking/'
    os.makedirs(pathpark, exist_ok=True)

    imgdir = os.path.dirname(pathpark)
    imgname = os.path.basename(i)
    newfname = imgdir + "/out_" + imgname
    print(newfname)
    pil_img = Image.fromarray(gray_img_masked)
    pil_img.save(newfname)
    

そうして抜き出した領域が左側の画像です。

右の画像から、抜き出した範囲が分かると思います。

因みに、右の画像は左側の画像を反転させてます。

コードは以下の通りです。

path = './hantai/'
os.makedirs(path, exist_ok=True)

imagedata = []
for t in glob.glob("./resize_image_sisaku/*.jpg"):
    jpg22 = cv2.imread( t, cv2.IMREAD_GRAYSCALE)
    gray_img_masked2_1 = jpg22 * (mask)

    io.imshow(gray_img_masked2_1)
    
    imgdir4 = os.path.dirname(path)
    imgname4 = os.path.basename(t)
    newfname4 = imgdir4 + "/gyaku_" + imgname4
    print(newfname4)
    pil_img4 = Image.fromarray(gray_img_masked)
    pil_img4.save(newfname4 )

これで、どの範囲を抜き取っているか画像ごとに確認できます。

解析に加えられない画像の判断はこれでしました。

例えば、以下の画像、左側は正しく抜き取られるように見えても、抜き取った範囲は違うと右の画像で確認できますね。

同様にして、白のピクセル数もカウントします。

4.解析結果

- オープン前の立体駐車場を見てみる

今回ターゲットにするディズニーランドの新立体駐車場は2019年、7月20日にオープンしました。

ということは、先ほどのこの画像は、オープン前のものとなります。

なので、7月8日の駐車場にて、黒と白のピクセル数を基準とし他の画像のピクセル数を比較することで空き具合、混み具合を判断したいと思います。

つまり、凄く空いていたら黒のピクセル数は7月8日に近い値になり、

逆に凄く混んでいたら白のピクセル数は7月8日のピクセル数を大幅に上回るだろうという想定です。

- 時系列で駐車状況を見てみる

実際に、12枚の画像を解析し時系列でまとめてみたグラフが以下になります。

これだけでは山も谷もランダムにあるように見え、何か傾向を見出すことは難しそうです。

- 月ごとにまとめた棒グラフでどの月が一番空いているのか見てみる

次に、これらの結果を月ごとにまとめてみます。

どちらのグラフでも一番左の7月の値が、オープン前の基準値（車が何も止まっていない）になります。ですが、グラフを見ると、黒いピクセルを数えたグラフでも基準値を上回る場合があったり、白いピクセルを数えた場合でも基準値を下回るケースがあることが分かります。

先に説明した通り、取得範囲のズレなどがあり、上手く意図した範囲が切り取れていないのではないかと考えられます。

実際に混んでいるであろうクリスマス直前の祝日である12/23の画像を見てみると以下のような感じで、駐車場は真っ黒、つまり空いているように見えました。。。このあたりはもう少し検討が必要なようです。

もろもろまだ課題はある＆解析したデータが少ないという点はありますが、無理やり考察をしてみると、以下のようなことはいえそうです。

・平日の方が比較的空いている。休日は混んでいる。
・比較的空いているのは9月、混んでいるのは8月

5.まとめ

今回のチャレンジで、できたことと出来なかったことをまとめると以下の様になります。

＜できたこと＞
-衛星画像の切り抜き
-切り抜いた画像のあるピクセル値の範囲におさまっているピクセル数を数える

＜できなかったこと＞
-Tellus OSでの図形描画から所定の位置の画像を切り抜く（位置ずれがあったため）
-衛星画像の位置ずれの自動補正（手で修正をするのは限界がある&時間が膨大にかかる）
-大量の衛星画像の解析（手作業での修正が発生してしまったため）
-SAR画像からの駐車場の混み具合の判定（上記、位置ずれの補正が難しかった&解析枚数が増やせなかったことによる）

なかなか思うようには行かなかった解析ですが、位置情報がうまく合わせられれば、様々な街の変化を捉えられる解析になるのではないかと考えられます。

とっつきにくいSAR画像データですが、皆さんもぜひ触ってみてください。

最後にpythonのファイルを貼り付けておきます。