はじめに

前回に引き続き、夜間光を用いてCOVID-19による影響を見てみたいと思います。
前回は日本全体のコロナによる影響の調査のため、日本全体のSum Of Light(SOL)の推移、各都道府県の平均輝度の増減を可視化しました。

今回は東京と大阪に絞って、夜間光の推移を見てみます。

東京と大阪の夜間光の推移

東京と大阪の夜間光の推移を見ていきます。
今回は2014年から2021年の東京と大阪の夜間光の輝度平均の推移を見てみます。

前回の日本全体のときには合計値であるSOLの推移を見ました。しかし合計値は面積が広いほど大きくなってしまい比較地域の面積が異なる場合には比較ができません。今回のように比較地域の面積が異なる場合は輝度平均が良さそうです。

まずはモジュールをインポートします。

import numpy as np
import pandas as pd
import warnings
warnings.filterwarnings('ignore')

try:
    import geemap, ee
    import seaborn as sns
    import matplotlib.pyplot as plt
    import geopandas as gpd

except ModuleNotFoundError:
    if 'google.colab' in str(get_ipython()):
        print("package not found, installing w/ pip in Google Colab...")
        !pip install geemap seaborn matplotlib
        !pip install geopandas
    else:
        print("package not found, installing w/ conda...")
        !conda install mamba -c conda-forge -y
        !mamba install geemap -c conda-forge -y
        !conda install seaborn matplotlib -y
    import geemap, ee
    import seaborn as sns
    import matplotlib.pyplot as plt
    import geopandas as gpd

GEEのデータを使用するため、GEEの認証を行います。
GEE認証を実行すると認証に必要なURLが表示されるので、URLへアクセスしてGoogleアカウントを指定し、表示された認証コードをGoogle Colabのボックスへコピペしてください。

try:
        ee.Initialize()
except Exception as e:
        ee.Authenticate()
        ee.Initialize()

今回扱うデータ、VIIRSの月間コンポジットのImageCollectionを取得します。

viirs = ee.ImageCollection("NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG").filterDate("2014-01-01","2021-12-31").select('avg_rad')

次に都道府県のgeometryの取得を行うファイルパスと関数を定義します。
今回は以下のサイトで公開されている都道府県ごとのgeojsonからgeometryを取得します。
https://japonyol.net/editor/article/47-prefectures-geojson.html
国土交通省のサイトでも行政区域のshpファイルやgeojsonが公開されていますが、市区町村レベルでのgeometryとなっており、今回のように都道府県レベルのgeometryにするには少し手間がかかるので、今回は簡単のためこちらのデータを用いることにしました。

まず、サイトからgeojsonをダウンロードし、Google Driveにアップロードします。
アップロードしたファイルのパスを取得し、下記コードに貼り付けてください。
geometryを取得する関数は、引数にgeojsonのファイルパスと都道府県名を与えることで、
指定した都道府県のgeometryを返すようになっています。
geopandasを用いてgeojsonを読み込み、['properties']['name']の中身が指定した都道府県と一致するもののみを取り出して、GEE用のgeometryオブジェクトに変換しています。

prefectures_shp_path =  'prefectures.geojsonのファイルパス'

def get_prefectures_geo(prefectures_shp_path, prefectures):
 geodf = gpd.read_file(prefectures_shp_path, crs='EPSG:4326')
 df = json.loads(geodf.to_json())
 features = [f for f in df['features'] if f['properties']['name']==prefectures]
 return ee.Geometry.MultiPolygon(features[0]['geometry']['coordinates'])

平均輝度を算出する関数を定義します。
引数で受け取ったimage(引数ではimgと記載)の指定された関心領域(aoi)に対して、ee.Reducer.mean()を適用し、imgに日付(date)、都道府県名(prefecture)、平均値(mean)をsetして返すようになっています。
reducerについては前回記事でも少し触れているのでご覧ください。

def aoi_mean(img):
    mean = img.reduceRegion(reducer=ee.Reducer.mean(), geometry=aoi, scale=500).get('avg_rad')
    return img.set('date', img.date().format()).set('prefecture', prefecture).set('mean',mean)

東京都の輝度の平均値を取得します。
prefectureに”東京都”を指定して、上で定義した関数get_prefectures_geo()を使って東京都のgeometryを得て関心領域(aoi)とします。
その後のreducerの適用〜データフレーム作成までの流れについても前回記事で触れているのでご覧ください。

prefecture = '東京都'
tokyo = get_prefectures_geo(prefectures_shp_path, prefecture)
aoi = tokyo

aoi_reduced_imgs = viirs.map(aoi_mean)
nested_list = aoi_reduced_imgs.reduceColumns(ee.Reducer.toList(3), ['date','prefecture','mean']).values().get(0)
df_tokyo = pd.DataFrame(nested_list.getInfo(), columns=['date','prefecture','mean'])

df_tokyo

2014年から2021年までの月別での東京都の平均輝度が得られました。

次に大阪の輝度の平均値を取得します。
prefectureに”大阪府”を指定して、東京都のときと同じ処理を行います。

prefecture = '大阪府'
osaka = get_prefectures_geo(prefectures_shp_path, prefecture)
aoi = osaka

aoi_reduced_imgs = viirs.map(aoi_mean)
nested_list = aoi_reduced_imgs.reduceColumns(ee.Reducer.toList(3), ['date','prefecture','mean']).values().get(0)
df_osaka = pd.DataFrame(nested_list.getInfo(), columns=['date','prefecture','mean'])

df_osaka

同じく2014年から2021年までの月別での大阪府の平均輝度が得られました。

東京都と大阪府のデータフレームを結合します。
また、データフレームの中身に対して以下の2点の処理を行います。
・dateの形式をyyyy-mm-ddに変更
・都道府県名を漢字からローマ字に変更
この処理は行わなくてもデータに影響はありませんが、後ほど行うグラフの表示を美しくするために行っておきます。

df = pd.merge(df_tokyo, df_osaka, how="outer")

df['date'] = pd.to_datetime(df['date'])
df.set_index('date', inplace=True)

df.replace('東京都', 'Tokyo', inplace=True)
df.replace('大阪府', 'Osaka', inplace=True)

df

無事結合することができました。

東京都と大阪府の平均輝度の推移を折れ線グラフで表示します。

fig, ax = plt.subplots(figsize=(15,7))

for city in ['Tokyo','Osaka']:
   data = df.loc[df['prefecture']==city]
   sns.lineplot(data=data, x=data.index, y='mean', label=city, ax=ax)

ax.set_ylabel('mean radiance',fontsize=20)
ax.set_xlabel('date',fontsize=20)
ax.legend(fontsize=20)
ax.set_title(f'Monthly mean radiance',fontsize=20)

グラフを見てみると、周期的なピークが見られます。前回記事で見た日本のグラフと同じように、どの年も12月1月あたりでピークとなっていそうです。
東京都は2020年に大きく減少しています。これはコロナによる影響と思われます。大阪府は2020年で大きく減少したという様子は見られません。

時系列的に平滑化して全体的な傾向を見てみます。
今回は12ヶ月間の移動平均を見てみます。

fig, ax = plt.subplots(figsize=(15,7))
window=12

for city in ['Tokyo','Osaka']:
   data = df.loc[df['prefecture']==city]
   sns.lineplot(data=data.rolling(window).mean(), x=data.index, y='mean', label=city, ax=ax)

ax.set_ylabel('mean radiance',fontsize=20)
ax.set_xlabel('date',fontsize=20)
ax.legend(fontsize=20)
ax.set_title(f'Monthly mean radiance',fontsize=20)

全体的な傾向としても、2020年は東京都では大きく減少しており、大阪府ではあまり減少していない様子が分かります。
東京都と大阪府でここまで変化に差がある理由はなんでしょうか…。
大阪工業大学で行われている夜間光データを用いた人間活動の把握の研究でも、コロナ前後で「東京都では輝度が減少した地域が多い一方で，大阪府や愛知県では増加しているなど，西高東低となる地域的な偏りが顕著となった。」と記載がありました。

この理由も今後詳細に調査してみると面白いかもしれません。

東京と大阪の増減の可視化

COVID-19が蔓延しはじめた2020年と、蔓延直前の2019年の東京都の夜間光の増減を表した画像を作成してみます。
まずは2019年と2020年の月間コンポジットのImageCollectionを取得し、それらのピクセルごとの中央値をもったimageを作成することでそれぞれの年間コンポジットを作成します。

viirs_2019 = ee.ImageCollection("NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG").filterDate("2019-01-01", "2019-12-31").select('avg_rad').median()
viirs_2020 = ee.ImageCollection("NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG").filterDate("2020-01-01", "2020-12-31").select('avg_rad').median()

2020年と2019年の画像の差を取り、東京のgeometryで切り取ります。
その差の画像に対して、増減をわかりやすくするため、値が0よりも大きいピクセルは全て値を1に、値が0よりも小さいピクセルは全て値を-1にします。

img_slope_tokyo = viirs_2020.subtract(viirs_2019).clip(tokyo)
img_slope_tokyo1 = img_slope_tokyo.gt(0).subtract(img_slope_tokyo.lt(0))

それでは可視化をしてみます。

viz_params = {'min':-1,
             'max':1,
             'palette':['1d4877','1b8a5a','ee3e32']}

map1 = geemap.Map(center=[35.68944,139.69167],zoom=10)
map1.addLayer(img_slope_tokyo1, viz_params, "2019-2020 VIIRS-DNB year rate of change")
map1.addLayerControl()
map1

東京都の地域ごとの増減がわかる画像が得られました。
青い場所は減少、赤い場所は増加を表しています。
都心寄りの多くの場所で輝度が減少している様子がわかります。
しかし、今回は増減と言っても増加か減少かしかみておらず、どの程度増減したのかについては可視化できていません。
この点はそれぞれの市区町村単位でSOLや輝度平均の推移をみて行く必要がありそうです。

個人的に気になった二つの地域について少し調査してみました。
一つ目は「稲城市周辺」(図中左)です。正直なんでこんなところが増えてるんだ？と思ったのですが、日本経済新聞の記事によると、稲城市は土地区画整理が進んでおり地価が9年連続上昇しているようです。
地域ごとの地価と夜間光の関係について調査してみるのも面白いなと思いました。

二つ目は「江東区周辺」(図中右)です。
この辺りはオリンピックの会場が多く、選手村のあった場所です。2020年といえばオリンピックの準備が行われていたはずなのでその影響でしょうか。とはいえ工事のスケジュールを見てみると、2019年には竣工予定なので予定通り竣工していれば工事の影響ではないようです。これも詳細に調査をしてみると面白そうです。

大阪府についても、東京都と同じく増加した箇所は1、減少した箇所は-1となる画像を作成し、可視化してみます。

img_slope_osaka = viirs_2020.subtract(viirs_2019).clip(osaka)
img_slope_osaka1 = img_slope_osaka.gt(0).subtract(img_slope_osaka.lt(0))

viz_params = {'min':-1,
            'max':1,
            'palette':['1d4877','1b8a5a','ee3e32']}

map2 = geemap.Map(center=[34.68639,135.52],zoom=9)
map2.addLayer(img_slope_osaka1, viz_params, "2019-2020 VIIRS-DNB year rate of change")
map2.addLayerControl()
map2

大阪府の地域ごとの増減がわかる画像が得られました。
青い場所は減少、赤い場所は増加を表しています。
東京都と比較すると、減少している地域が少ないことがわかります。
しかし、上でも述べたように増減と言っても増加か減少かしかみておらず、どの程度増減したのかについては可視化できていません。
東京と大阪でこのような差がある理由について、もし意見がある方がいらっしゃいましたらコメントいただけると大変喜びます。

まとめ

今回は東京都と大阪府に絞って、夜間光を用いてCOVID-19による影響を見てみました。
今回は夜間光の推移や増減を可視化しただけでしたが、これらのデータとその他のデータとの相関をとってみると、よりおもしろい考察ができそうです。

GEEは簡単に衛星データに触れることができ、プログラミングやデータ解析の練習もできるのでオススメです。
これからも勉強しながら成果を記事にしていこうと思うので、よろしくお願いします。

参考文献

[1]Google Earth Engineデータカタログ
https://developers.google.com/earth-engine/datasets/catalog/NOAA_VIIRS_DNB_MONTHLY_V1_VCMSLCFG

[2]大阪工業大学「夜間光データを用いた人間活動の把握」
https://www.research.oit.ac.jp/oitid/seeds/seeds/seeds-11005/

[3]日本経済新聞「地価回復、コロナの爪痕なお」
https://vdata.nikkei.com/newsgraphics/land-value-map_analysis2022-k6r8kj/

[4]江東おでかけ情報局「東京2020オリンピック・パラリンピック江東区内会場配置MAP」
http://www.koto-kanko.jp/special/detail.php?spid=46