夜間光データでコロナによる影響を調査してみた〜①日本全体編〜

## はじめに
今回は夜間光のデータを用いて、日本のコロナによる影響を調査してみます。
今回もGoogle Earth Engine(GEE)とGoogle Colabを用いて解析を行っていきます。
「まずそれなに？」という方は、以前初学者向けに書いた[登りたい山を探す企画の記事](https://qiita.com/oz_oz/items/dd5a69323f499a6dfd20)があるので、ぜひご覧ください。
今回は、新型コロナウイルス(COVID-19)によって、日本の夜間光にどのような変化があったかを調査していきます。
統計学的な知識に疎いので、今回は簡単に夜間光の推移を見る程度ですが、いずれは相関など詳細な調査を行いたいと考えています。
また、夜間光データに関しては[World Bank](https://worldbank.github.io/OpenNightLights/welcome.html)がチュートリアルを公開していますので、興味を持った方はぜひそちらで勉強してみてください。


## 夜間光データについて
夜間光のデータには、主にDMSPとVIIRSと呼ばれるデータがあります。
GEEではDMSPの1992年から2013年の年間コンポジット、VIIRSの2014年以降の月間コンポジットが公開されています。
今回はCOVID-19による影響を見たいため、2019年から2021年のデータがあるVIIRSのデータを使います。
GEEでも以下のURLのようにVIIRSのデータが公開されています。
https://developers.google.com/earth-engine/datasets/catalog/NOAA_VIIRS_DNB_MONTHLY_V1_VCMSLCFG

## 夜間光データと経済
夜間光のデータは経済活動との相関があると言われており、既に経済指標の一つとして用いられています。
以下は夜間光が観測された面積とGDPの関係を表した図[[3]](#参考文献)です。
この図から、夜間光と経済活動には正の相関がありそうだとわかると思います。
![スクリーンショット 2022-05-08 17.48.23.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2152045/c8d846be-f145-19ef-f786-29d1ec0e4f67.png)
Source :Christopher D Elvidge, Marc L Imhoff, Kimberly E Baugh, Vinita Ruth Hobson, Ingrid Nelson, Jeff Safran, John B Dietz, and Benjamin T Tuttle. Night-time lights of the world: 1994–1995. ISPRS Journal of Photogrammetry and Remote Sensing, 56(2):81–99, 2001.

また、GEEで2021年12月の日本周辺の夜間光データを可視化してみました。
実際に東京、愛知、大阪、福岡、札幌等大きな都市があるあたりは光が強い様子が確認できます。
![スクリーンショット 2022-04-20 21.23.52.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2152045/c9a083ef-d210-8af5-d443-8f2b6a9a4aa4.png)

## 日本の夜間光の推移
日本の夜間光の推移を見ていきたいと思います。
今回は2014年から2021年の日本の夜間光の合計値Sum Of Light(SOL)の推移を見てみます。

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

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

# reminder that if you are installing libraries in a Google Colab instance you will be prompted to restart your kernal

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のボックスへコピペしてください。

```python
try:
        ee.Initialize()
except Exception as e:
        ee.Authenticate()
        ee.Initialize()
```
今回扱うデータ、VIIRSの月間コンポジットのImageCollectionを取得します。
```python
viirs = ee.ImageCollection("NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG").filterDate("2014-01-01","2021-12-31").select('avg_rad')
```
日本のgeometryを取得します。
shpファイルやgeojsonファイルがある場合は、それをインポートしてgeometryとすることができます。
しかし、GEEにはFeatureCollectionとして日本のgeometryがあるため、今回はそれを利用します。

```python
japan = ee.Feature(ee.FeatureCollection("FAO/GAUL/2015/level0").filter(ee.Filter.eq('ADM0_NAME', 'Japan')).first()).geometry()
```

日本のSOLを算出する関数を定義します。
引数で受け取ったimage(引数ではimgと記載)に対して、ee.Reducer.sum()を適用し、imgに日付と合計値をsetして返すようになっています。
reduceRegion()に、reducer、geometry、scaleを与えることで、image内の値を操作できます。操作の種類はsum()、max()、min()、meam()、median()などがあります。
今回は光の合計値を取得したいのでsum()を適用します。
![スクリーンショット 2022-04-20 21.24.12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2152045/59487828-905e-0429-d415-ab7ce4b6c3a5.png)

Source : https://developers.google.com/earth-engine/guides/reducers_reduce_region

```python
def get_japan_sum(img):
   sum = img.reduceRegion(reducer=ee.Reducer.sum(), geometry=japan, scale=500, maxPixels=1e9).get('avg_rad')
   return img.set('date', img.date().format()).set('sum',sum)
```

先ほど定義した関数をviirsのImageCollectionに対して適用します。
リスト内の要素に同じ処理をしたい際、map関数を用いることで全てのImageに同じ処理を行うことができます。GEEでは基本的にfor文を使うことはないようです。ImageCollectionはimageの集まりなので、それぞれにfor文で上の処理を行いたくなりますが、GEEではfor文は推奨されていません。高度な並列処理ができなくなるということが理由らしいです。
詳細は[公式のドキュメント](https://developers.google.com/earth-engine/tutorials/tutorial_js_03)を参照してください。

ImageCollectionはImageの集まりで、その一つ一つのImageにmap関数を用いてget_japan_sum関数を適用しています。
この方法を適用したImageCollectionからdate(日付)とsum(合計値)を取り出して、リストにします。ImageCollectionから、指定した列のみを抽出したい場合、reduceColumnsを利用します。reduceColumnsでは出力は辞書型になるので、ee.Reducer.toListを用いてリスト型に変換し、最後にgetInfo()で値のみを抽出しています。

```python
japan_sum = viirs.map(get_japan_sum)
nested_list = japan_sum.reduceColumns(ee.Reducer.toList(2), ['date','sum']).values().get(0)
df = pd.DataFrame(nested_list.getInfo(), columns=['date','sum'])
```
後ほどグラフにするときのために、dateの表記をわかりやすく変えておきます。

```python
df['date'] = pd.to_datetime(df['date'])
df = df.set_index('date')
df
```

出力は以下のようになりました。
2014年から2021年の月ごとの合計値のリストが取得できました。
![スクリーンショット 2022-04-20 21.32.38.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2152045/51c4c74d-8ef9-af7c-8b86-d0cf5bb27c47.png)

折れ線グラフで出力してみます。

```python
fig, ax = plt.subplots(figsize=(15,7))
sns.lineplot(data=df, ax=ax)
ax.set_ylabel('SOL',fontsize=20)
ax.set_xlabel('Date',fontsize=20)
ax.set_title('Monthly Sum Of Lights (SOL) for Japan (2014 to 2021)',fontsize=20)
```
グラフを見てみると、周期的なピークが見られます。ざっくりですが、どの年も12月1月あたりでピークとなっていそうです。月ごとに集計しているので、冬至(一番夜が長い)があるため夜間光が増えているのでしょうか。また、クリスマスや年末年始の影響なのでしょうか。これも詳細に考察してみると面白そうですね。

今回のメインテーマであるCOVID-19の影響を見てみます。
2020年のピークが2019年のピークに比べて下がっていることがわかります。これはCOVID-19による緊急事態宣言等の影響と考えられます。
![スクリーンショット 2022-04-20 21.32.50.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2152045/93885af0-8f29-7224-ec2f-43f1ba613eb0.png)

周期的なピークが見られているので、時系列的に平滑化して全体的な傾向を見てみます。
今回は12ヶ月間の移動平均を見てみます。

```python
fig, ax = plt.subplots(figsize=(15,7))
window=12
sns.lineplot(data=df.rolling(window).mean(), ax=ax)
ax.set_ylabel('SOL',fontsize=20)
ax.set_xlabel('Date',fontsize=20)
ax.set_title(f'Monthly SOL ({window} moving avg.) for Japan (2014 to 2021)',fontsize=20)
```

全体的な傾向としても、2020年は大きく減少している様子が確認できます。
COVID-19の影響で夜間光が減少したと言っても良いのではないでしょうか。
その他にも2018年には大きく増加している様子が確認されます。
何があったのでしょうか…。ここも詳しく考察してみると面白いかもしれません。
![スクリーンショット 2022-04-20 21.33.05.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2152045/4795fc24-3597-34be-caef-c9c9b64faab2.png)

下のグラフは帝国データバンクが公開している全国の景気DIのグラフです。
景気DIとは、景気全体の現状を把握したり、将来の動向を予測したりするときに使われる経済指標です。
これを見てみると、2018年に景気DIは過去最高となり、2020年の緊急事態宣言時には大きく減少している様子がわかります。
12ヵ月間の移動平均のグラフと比較してみると、概ね同じような推移をしているように思えます。今回は統計的な検証は行いませんが、これも今後詳細に考察すると面白いかもしれません。
![スクリーンショット 2022-04-20 21.33.57.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2152045/3ace9598-7e93-5f9f-11ff-bd2c57f02746.png)
Source :https://www.tdb-di.com/economic-trend-survey/ets202202.php

## 都道府県ごとの増減の可視化
COVID-19が蔓延しはじめた2020年と、蔓延直前の2019年の夜間光の増減を都道府県ごとのコロプレスを作成して見てみます。コロプレスとは、下の画像のように領域ごとに色を塗ったような図です。
![chizu_color.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2152045/4937de63-5c15-deee-dae0-ec08f888921b.png)
2019年と2020年の月間コンポジットのImageCollectionを取得し、それらのピクセルごとの中央値をもったimageを作成することでそれぞれの年間コンポジットを作成します。
下の画像[[6]](#参考文献)のような処理をします。
![スクリーンショット 2022-04-20 21.34.50.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2152045/7477da7e-4402-6ada-6d28-8540ac56989e.png)
Source :https://developers.google.com/earth-engine/guides/reducers_image_collection

ImageCollectionの後ろに.median()をつけることで、上で述べたような処理ができます。
```python
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()
```
取得した2019年と2020年の差をとります。
ピクセル毎に2020年と2019年の差を計算した値を持つimageを作ります。
よって、これの値がプラスであれば2019年から2020年にかけて夜間光は増加、値がマイナスであれば2019年から2020年にかけて夜間光は減少していることになります。
```python
viirs_19_20 = viirs_2020.subtract(viirs_2019)
```
次に都道府県毎のgeometryを取得します。
上記同様、GEEにはFeatureCollectionとして日本の都道府県レベルでのgeometryがあるため、今回はそれを利用します。
```python
japan_prefectures = ee.FeatureCollection("FAO/GAUL/2015/level1").filter(ee.Filter.eq('ADM0_NAME', 'Japan'))
```
先ほど取得した2019年と2020年の差のimageで、都道府県ごとにピクセルの平均値を計算します。
上での日本の合計値のように、一つのgeometry内の代表値を取得する場合にはreduceRegionを使いました。
今回のように複数のgeometry内の代表値を取得する場合にはreduceRegionsを使います。
reducerにはee.Reducer.mean()を与えて各都道府県ごとに平均値を計算します。
その後、geometryごとに平均値を持ったee.FeatureCollectionに変換します。
```python
scaleFactor=100

japan_pref = ee.FeatureCollection(viirs_19_20.reduceRegions(reducer=ee.Reducer.mean(),
                                       collection=japan_prefectures,
                                       scale=scaleFactor))
```

コロプレスを作成します。
上で作成したee.FeatureCollectionはgeometryの集まりなので、これをimageに変換します。
都道府県ごとにgeometry内のピクセルに平均値の値を与えていく、塗り絵のような処理だと考えるとわかりやすいかもしれません。
reducerには一番最初の値を取得するee.Reducer.first()を与えています。今回はピクセルに与えたいmeanは単一の値で格納されているため、ee.Reducer.first()が適切です。仮にここにリストのような複数の値が格納されている場合は、ee.Reducer.sum()やee.Reducer.mean()のようなreducerを与えることで、単一の値を取得することができます。
```python
radiance_img = japan_pref.reduceToImage(properties=['mean'],reducer=ee.Reducer.first())
```
増減をわかりやすくするため、値が0よりも大きいピクセルは全て値を1に、値が0よりも小さいピクセルは全て値を-1にします。
```python
radiance_img1 = radiance_img.gt(0).subtract(radiance_img.lt(0))
```

それでは可視化をしてみます。
```python
japanMap = geemap.Map()
japanMap.centerObject(japan_prefectures, zoom=5)
viz_params = {'min':-1,
             'max':1,
             'palette':['1d4877','1b8a5a','ee3e32']}
japanMap.addLayer(japan_pref, {}, "Prefecture boundaries", opacity=1)
japanMap.addLayer(radiance_img1, viz_params, 'VIIRS Japan (2019 to 2020) avg rad by prefecture',opacity=1)
japanMap.addLayerControl()
japanMap
```

青い場所は減少、赤い場所は増加を表しています。
東京、愛知、大阪、福岡等大きな都市は減少している様子が分かります。
しかし、今回は増減と言っても増加か減少かしかみておらず、どの程度増減したのかについては可視化できていません。
この点はそれぞれの都道府県単位でSOLや輝度平均の推移をみて行く必要がありそうです。
![スクリーンショット 2022-04-20 21.33.15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2152045/0c935b98-5c44-3405-c5ec-f51157a01d8c.png)

## まとめ
今回は夜間光データを用いて、日本全体のCOVID-19による影響を調査してみました。
詳細な考察はできていませんが、実際にCOVID-19によって夜間光が減少している様子が分かりました。また、経済指標の一つである景気DIと夜間光SOLが似たような推移になっている様子もわかりました。
次回は東京と大阪の影響をみてみたいと思います。
[夜間光データでコロナによる影響を調査してみた〜②東京大阪編〜](https://qiita.com/oz_oz/items/42ce245ba16b75595a13)

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

[2]World Bank  ”Open Night Lights tutorial”
https://worldbank.github.io/OpenNightLights/welcome.html

[3]Christopher D Elvidge, Marc L Imhoff, Kimberly E Baugh, Vinita Ruth Hobson, Ingrid Nelson, Jeff Safran, John B Dietz, and Benjamin T Tuttle. Night-time lights of the world: 1994–1995. ISPRS Journal of Photogrammetry and Remote Sensing, 56(2):81–99, 2001.

[4]宙畑「衛星が撮影した夜の地球「夜間光」がお金に変わる!? 概要と利用事例」
https://sorabatake.jp/240/

[5]帝国データバンク「2022年2月の景気動向調査」https://www.tdb-di.com/economic-trend-survey/ets202202.php

[6]ImageCollection Reductions
https://developers.google.com/earth-engine/guides/reducers_image_collection