夜間光データで2016年熊本地震の復興状況を調査してみた

はじめに

本記事では夜間光データを用いて、2016年に起きた熊本地震の復興状況を調査してみます。
今回もGoogle Earth Engine(GEE)とGoogle Colabを用いて解析を行っていきます。
「まずそれなに？」という方は、以前初学者向けに書いた登りたい山を探す企画の記事があるので、ぜひご覧ください。
また、夜間光データについても、以前書いた夜間光データでコロナの影響を調査する企画があるため、ぜひご覧ください。

※2016年熊本地震
4月14日、熊本県熊本地方においてマグニチュード6.5の地震が発生し、熊本県益城町で震度7を観測しました。また、16日にはマグニチュード7.3の地震が発生し、益城町及び西原村で震度7を、熊本県を中心にその他九州地方の各県でも強い揺れを観測しました。

熊本の夜間光画像を取得する

熊本地震が発生した2016年4月の1ヶ月前の2016年3月から、2016年7月の夜間光画像を取得します。

まずは必要なライブラリをインポートします。

import numpy as np
import pandas as pd
import warnings
import json
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のボックスへコピペしてください。

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

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

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

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'])

prefectures_shp_path =  'ここに貼り付ける'

prefecture = '熊本県'
kumamoto = get_prefectures_geo(prefectures_shp_path, prefecture)
aoi = kumamoto

今回は熊本地震が発生した2016年4月の1ヶ月前の2016年3月から、2016年7月の夜間光画像を取得します。
取得範囲を定めて、forループで画像を取得していきます。

fdatelist = ['2016-03-01', '2016-04-01','2016-05-01','2016-06-01','2016-07-01']
ldatelist = ['2016-03-31', '2016-04-30','2016-05-31','2016-06-30','2016-07-31']

img_list = []
for i in range(5):
 img = ee.ImageCollection("NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG").filterDate(fdatelist[i], ldatelist[i]).select('avg_rad').median().clip(aoi)
 img_list.append(img)

取得した画像を表示します。
img_list[0]となっている部分を0~4で変えて実行すると、2016年3月〜2016年7月の画像が表示されます。

map1 = geemap.Map(center=[32.78972,130.74167],zoom=9)
map1.addLayer(img_list[0], name='img')
map1.addLayerControl()
map1

それぞれの画像を出力した結果です。
地震の起きた4月は少し明るい箇所が減ったように見えます。しかしこのように画像を見比べても変化がわかりづらいので、次の章では差分画像を取って可視化をしてみます。

図1　2016年3月〜2016年7月の熊本県の夜間光画像

差分画像の取得

1ヶ月ごとの差分画像を作成し、変化の様子を見てみます。

衛星リモートセンシングの変化検出に関して、以下のような手法が宙畑の記事で紹介されていましたので、そちらを参考にします。

「衛星リモートセンシングの変化検出の手法のうち、最もシンプルな方法として、差分画像に対しμ±2σ（μ:平均値; σ:標準偏差）を上限・下限として、それらを超える画素値を変化とみなす方法があります。正規分布においては、μ±2σの区間に全データの95%が含まれ、上限・下限を超えるデータは2.5%ずつになります。2.5%は全体の傾向に比べて異常に大きい・小さいことから、何かしらの大きな変化が地上にあったとみなします。これにより対象領域の全体的なバイアスを除くことができます。」(宙畑記事より抜粋)

2時期の差分画像を引数として与え、以下の条件でピクセル値を置き換えた画像を返す関数を定義します。
・μ+2σより大きい値を持つピクセル値を2
・μ-2σより小さい値を持つピクセル値を1
・μ-2σ以上、μ+2σ以下の値を持つピクセル値を0

ee.image.reduceRegion()に、reducer、geometry、scaleを与えることで、画像の代表値を取得することができます。reducerには平均値μを算出するmeam()、標準偏差σを算出するstdDev()を指定します。

def make_diff_img(img_before, img_after):
 img_diff = img_after.subtract(img_before).clip(aoi)
 mean = ee.Number(img_diff.reduceRegion(reducer=ee.Reducer.mean(), geometry=aoi, scale=30).get('avg_rad')).getInfo()
 stdDev = ee.Number(img_diff.reduceRegion(reducer=ee.Reducer.stdDev(), geometry=aoi, scale=30).get('avg_rad')).getInfo()
 range_min = mean - 2* stdDev
 range_max = mean + 2* stdDev
 img_diff_mask = img_diff.gt(range_max).add(img_diff.gt(range_max)).add(img_diff.lt(range_min))
 return img_diff_mask

1ヶ月毎の差分を計算し、μ±2σを閾値として変化を検出する画像を取得します。

img_diff_mask_list = []
for i in range(0, 4):
  img_before = img_list[i]
  img_after = img_list[i+1]
  img_diff_mask_list.append(make_diff_img(img_before, img_after))

取得した画像を表示します。
こちらも同じくimg_diff_mask_list[0]となっている部分を0~3で変えて実行すると、それぞれに対応した画像が表示されます。

viz_params = {'min':0,
            'max':2,
            'palette':['000000','1d4877','ee3e32']}

map2 = geemap.Map(center=[32.78972,130.74167],zoom=9)
map2.addLayer(img_diff_mask_list[0].mask(img_diff_mask_list[0]), viz_params, "img_diff")
map2.addLayerControl()
map2

図２はそれぞれの画像を出力した結果です。赤色が増加、青色が現象を表しています。
図３は熊本地震の震源地とその推計震度分布です。
地震の起きた2016年3月〜2016年4月では広い地域で減少している様子が分かります。また図３をみてみると、減少している地域は震源地付近であったことがわかります。やはり震源近くの被害が深刻な地域では、夜間光にも大きな変化が現れるということでしょうか。
それ以降は少しずつ増加している地域が増えており、復興・復旧が進んだと考えられます。しかし同時に広い範囲で減少も続いています。震災前後にはどのような要因で夜間光画像が変化しうるかは今後文献を調査する必要があるかなと思います。

図２　月別夜間光の差分から変化を検出した画像（赤色は増加、青色は減少)

図３　熊本地震の推計震度分布(震度７以上)
（出典）気象庁「平成28年（2016年）熊本地震について（第７報）」（平成28年４月16日）

まとめ

今回は夜間光データを用いて、2016年に起きた熊本地震の復興状況を調査してみました。震源近くの被害が深刻な地域では夜間光の減少も顕著であることがわかりました。それ以降は夜間光の増加も見られ、復興・復旧の様子が現れることもわかりました。しかし、地震発生から1ヶ月後以降も広い地域で減少が続いている様子が見られており、その要因は更なる調査が必要です。申し訳ありませんが、本記事では今後の課題としておきたいと思います。

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

参考文献

宙畑「【コード付き】復興状況を人工衛星からモニタリング！ 北海道地震後の夜間光から確認してみた」
https://sorabatake.jp/10526/

内閣府-防災情報のページ「特集1 - 平成28年熊本地震‐内閣府防災情報のページ」
https://www.bousai.go.jp/kohou/kouhoubousai/h28/83/special_01.html