１．はじめに

高い場所から見る夜景は私たちをロマンティックな気分にさせてくれます。
そして、宇宙からも夜景を見ることができます。

衛星から夜の地球を観測したデータは「夜間光データ」と呼ばれています。
現実的な表現になってしまいますが、明るい場所ほど多くの人間が活動しており、活発な地域ということになります。

以下の図は世界200ヶ国の夜間光面積とGDPの関係を表したグラフです。このように夜間光とGDP(国内総生産)には強い相関があり、実際に経済指標の一つとして扱われています。

本記事では、全国の市区町村を対象に夜間光データと土地価格の相関を調査します。

筆者は衛星データに関しては初心者で、宙畑の記事等を読んで絶賛勉強中です。

衛星データに関する情報を集めていたところ、国のGDP（国内総生産）と夜間光には相関があるという話題を見つけました。それならば国の自治体レベルの総生産、つまりは繁栄具合とも相関があるのか気になり、調査してみることにしました。

自治体が繁栄しているということは、それだけその自治体の土地を購入したいと思う人が多いはず。それならば土地価格に着目すると自治体ごとの繁栄具合がわかるのではないか、と考えて、今回は夜間光と土地価格の関係に興味を持ちました。

直感的には、明るい場所は都心部のため土地価格は高い、暗い場所は郊外のため土地価格が低いといったように相関がありそうです。また、夜間光の変化と土地価格の変化にも相関があるかもしれません。

そこで今回は「年間の夜間光の平均値と土地価格の相関」と「夜間光の平均値と土地価格の変化率の相関」の二つの項目について調査を行います。

２．夜間光データとは

夜間光データは簡単にいうと夜の地球を撮影したデータです。

現在最もよく使われている夜間光データの一つが、アメリカの衛星であるSuomi-NPPに搭載されているVisible/Infrared Imager and Radiometer Suite(VIIRS)によって撮影されたものです。

しかし、撮影したときに雲がかかっていたり、稲妻やオーロラが発生していると正しく地域の光量を観測することができません。そこでNOAAでは、VIIRSの画像を月間もしくは年間で集約することで雲などのノイズを除去した夜間光データを作成してくれています。

VIIRSは地上分解能が500 m程度と通常の光学データに比べて粗いですが、先ほど説明したように月間・年間でのノイズの除去された夜間光データとして配信されているため、よく利用されています。

夜間光データについては、宙畑でもこれまでに記事で紹介してきているのでぜひそちらもご覧ください。

３．解析の流れ

まず、解析環境を準備します。今回の解析にはGoogle Earth Engine(GEE)とGoogle Colabを使います。

GEEは、Googleが提供する人工衛星画像をサーバー上で解析をすることが出来るサービスです。研究・教育・非営利目的であれば無償で利用することができ、複雑なデータフォーマットへの対応や前処理が不要となり、初学者でも簡単に衛星データを扱うことができます。

今回利用するVIIRSのデータもGEEで扱うことができます。

Google ColabはGoogleが提供するブラウザからPython を実行できるサービスです。Google ColabからGEEのデータを扱うこともできます。

解析は以下の流れに沿って進めます。

① 市区町村ごとのポリゴンデータを使って夜間光の平均値を取得し、夜間光テーブルを作成　
② 市区町村ごとの土地価格の平均値のテーブルを作成する
③「夜間光の平均値と土地価格」と「夜間光の平均値と土地価格の変化率」の相関係数を算出

４．市区町村ごとの夜間光の平均値を求める

それではまずは全国の市区町村ごとに夜間光の平均値を求めてみます。

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

import numpy as np
import pandas as pd
import warnings
import json
from scipy import stats
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
   from matplotlib import cm

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
   from matplotlib import cm

次に、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").select('avg_rad')

今回は年間のデータで比較を行いたいため、月間コンポジットの中央値から年間コンポジットを作成します。VIIRSのデータは2014年から公開されているため、2014年から2021年までのデータを作成します。

def viirs_annual_median_reduce(year):
    return viirs.filter(ee.Filter.calendarRange(year,year,"year")).select("avg_rad").median().set('year',year)

start = 2014
end = 2021
years = ee.List.sequence(start, end)

viirs_year = ee.ImageCollection.fromImages(years.map(viirs_annual_median_reduce))

次に市区町村ごとのポリゴンデータを取得します。

ポリゴンデータとは、簡単に言えば緯度経度を持った形状のデータのことで、日本の市区町村ごとのポリゴンデータは国土交通省のサイトからダウンロードすることができます。サイトからダウンロードした全国の市区町村のポリゴンデータをGeopandasで読み込みます。
https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v3_0.html

filepath_geo = 'geojsonのファイルパス'

gdf = gpd.read_file(filepath_geo)

# GeoJSONファイルから、属性データの読み込み
with open(filepath_geo) as f:
 df = json.load(f)

夜間光データをまとめるテーブルを作成します。

夜間光テーブル(df_nl)と夜間光変化率テーブル(df_pnl)の２種類を用意し、それぞれ行政区域コード(code)、市区町村名(muni)、2014年から2021年までの夜間光の平均値、2014年から2021年までの夜間光の平均値の年間変化率を格納します。

cols_nl = ['code', 'muni', '2014-NL', '2015-NL', '2016-NL', '2017-NL', '2018-NL', '2019-NL', '2020-NL', '2021-NL']
cols_pnl = ['code', 'muni', '2014-2015_PNL', '2015-2016_PNL', '2016-2017_PNL', '2017-2018_PNL', '2018-2019_PNL', '2019-2020_PNL', '2020-2021_PNL']
df_nl = pd.DataFrame(index=[], columns=cols_nl)
df_pnl = pd.DataFrame(index=[], columns=cols_pnl)

市区町村ごとの夜間光データを取得していきます。

まず、対象地域の平均値を取得する関数(aoi_mean)を定義します。

次に読み込んだGeopandasから、市区町村名(munis)と市区町村コード(codes)のリストを取得します。

取得したリストにある市区町村に対して、ループ処理で夜間光の平均値を取得していき、夜間光テーブル(df_nl)と夜間光変化率テーブル(df_pnl)に追加していきます。

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

munis = gdf['N03_004'].unique()
codes = gdf['N03_007'].unique()

for i in range(len(munis)):
 code = codes[i]
 muni = munis[i]
 features = [[f['geometry']['coordinates'][0]] for f in df['features'] if f['properties']['N03_007']==code]
 if len(features) > 1000:
   continue
 aoi = ee.Geometry.MultiPolygon(features)
 poi_reduced_imgs = viirs_year.map(aoi_mean)
 nested_list = poi_reduced_imgs.reduceColumns(ee.Reducer.toList(2), ['year','mean']).values().get(0)
 df_temp = pd.DataFrame(nested_list.getInfo(), columns=['year','mean'])
 PNL_14_15 = df_temp.at[1, 'mean'] / df_temp.at[0, 'mean']
 PNL_15_16 = df_temp.at[2, 'mean'] / df_temp.at[1, 'mean']
 PNL_16_17 = df_temp.at[3, 'mean'] / df_temp.at[2, 'mean']
 PNL_17_18 = df_temp.at[4, 'mean'] / df_temp.at[3, 'mean']
 PNL_18_19 = df_temp.at[5, 'mean'] / df_temp.at[4, 'mean']
 PNL_19_20 = df_temp.at[6, 'mean'] / df_temp.at[5, 'mean']
 PNL_20_21 = df_temp.at[7, 'mean'] / df_temp.at[6, 'mean']
 df_pnl = df_pnl.append(pd.Series([code, muni, PNL_14_15, PNL_15_16, PNL_16_17, PNL_17_18, PNL_18_19, PNL_19_20, PNL_20_21], index=df_pnl.columns), ignore_index=True)
 df_nl = df_nl.append(pd.Series([code, muni, df_temp.at[0, 'mean'], df_temp.at[1, 'mean'], df_temp.at[2, 'mean'], df_temp.at[3, 'mean'], df_temp.at[4, 'mean'], df_temp.at[5, 'mean'], df_temp.at[6, 'mean'], df_temp.at[7, 'mean']], index=df_nl.columns), ignore_index=True)

夜間光テーブル(df_nl)と夜間光変化率テーブル(df_pnl)が作成できました。

５．市区町村ごとの土地価格の平均値を求める

続いて市区町村ごとの土地価格のテーブルを作成します。

土地価格のデータも国土交通省のサイトからダウンロードすることができます。全国の市区町村の土地価格のデータをGeopandasで読み込みます。
https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-L02-v3_0.html

fpath = '土地価格のgeojsonのファイルパス'

# ファイル読込（encoding='SHIFT-JIS'とすることで日本語データにも対応）
df_gjs = gpd.read_file(fpath,encoding='SHIFT-JIS')

土地価格データをまとめるテーブルを作成します。

土地価格テーブル(df_lp)と土地価格変化率テーブル(df_plp)の２種類を用意し、それぞれ行政区域コード(code)、2014年から2021年までの土地価格の平均値、2014年から2021年までの土地価格の平均値の年間変化率を格納します。

# 必要なデータのみ格納するための変数を宣言
df_plp = pd.DataFrame(columns = ['code', '2014-2015_PLP', '2015-2016_PLP', '2016-2017_PLP', '2017-2018_PLP', '2018-2019_PLP', '2019-2020_PLP', '2020-2021_PLP'])
df_lp = pd.DataFrame(columns = ['code', '2014-LP', '2015-LP', '2016-LP', '2017-LP', '2018-LP', '2019-LP', '2020-LP', '2021-LP'])

取得したデータは番地単位での土地価格(例：千代田区３丁目等)になっているため、行政区域コードごとの平均値を取得し、土地価格テーブル(df_lp)と土地価格変化率テーブル(df_plp)に値を追加していきます。

# 市区町村ごとの平均値を取得
group_df_gjs = df_gjs.groupby('L02_020', as_index=False).mean().loc[:,['L02_020', 'L02_086', 'L02_087', 'L02_088', 'L02_089', 'L02_090', 'L02_091', 'L02_092', 'L02_093']]

df_plp['code'] = group_df_gjs['L02_020']
df_plp['2014-2015_PLP'] = group_df_gjs['L02_087'] / group_df_gjs['L02_086']
df_plp['2015-2016_PLP'] = group_df_gjs['L02_088'] / group_df_gjs['L02_087']
df_plp['2016-2017_PLP'] = group_df_gjs['L02_089'] / group_df_gjs['L02_088']
df_plp['2017-2018_PLP'] = group_df_gjs['L02_090'] / group_df_gjs['L02_089']
df_plp['2018-2019_PLP'] = group_df_gjs['L02_091'] / group_df_gjs['L02_090']
df_plp['2019-2020_PLP'] = group_df_gjs['L02_092'] / group_df_gjs['L02_091']
df_plp['2020-2021_PLP'] = group_df_gjs['L02_093'] / group_df_gjs['L02_092']

df_lp['code'] = group_df_gjs['L02_020']
df_lp['2014-LP'] = group_df_gjs['L02_086']
df_lp['2015-LP'] = group_df_gjs['L02_087']
df_lp['2016-LP'] = group_df_gjs['L02_088']
df_lp['2017-LP'] = group_df_gjs['L02_089']
df_lp['2018-LP'] = group_df_gjs['L02_090']
df_lp['2019-LP'] = group_df_gjs['L02_091']
df_lp['2020-LP'] = group_df_gjs['L02_092']
df_lp['2021-LP'] = group_df_gjs['L02_093']

土地価格テーブル(df_lp)と土地価格変化率テーブル(df_plp)が作成できました。

６．夜間光と土地価格の相関を調べてみる

作成した夜間光のテーブルと土地価格のテーブルを使って、二つのデータの相関を調査していきます。

まずは夜間光テーブル(df_nl)と土地価格テーブル(df_lp)、夜間光変化率テーブル(df_pnl)と土地価格変化率テーブル(df_plp)を結合します。

df_nl_lp = pd.merge(df_lp, df_nl, how='inner')
df_pnl_plp = pd.merge(df_plp, df_pnl, how='inner')

夜間光と土地価格の相関係数を算出してみます。

今回は2014年と2020年の散布図と相関係数、p値の表示をします。

df_nl_lp.plot.scatter(x='2014-NL', y='2014-LP')

x1, y1 = df_nl_lp["2014-NL"], df_nl_lp["2014-LP"]
print(stats.pearsonr(x1, y1))

df_nl_lp.plot.scatter(x='2020-NL', y='2020-LP')

x2, y2 = df_nl_lp["2020-NL"], df_nl_lp["2020-LP"]
print(stats.pearsonr(x2, y2))

2014年と2020年ともに相関係数が0.7以上、p値が0に限りなく近いことから、相関があると言えそうです。

次に夜間光と土地価格の変化率の相関係数を算出してみます。変化率は同じ年同士の関係だけでなく異なる年との関係も見たいため、相関係数の一覧で表示します。

df_pnl_plp.corr()

同じ年の変化率で見てみると、相関係数は0前後となっており相関がなさそうです。また、2017から2018にかけて夜間光が増加していたら、2018から2019にかけて土地価格が上がるなど、夜間光の変化が少し遅れて土地価格に反映される可能性があるかと考ていましたが、異なる年の変化率で見てみても相関はなさそうです。

いくつかの市区町村の夜間光と土地価格の推移を比較してみましょう。

日本経済新聞のサイトで、土地価格について取り上げられていた「東京都千代田区」と「東京都稲城市」の推移を可視化してみます。
https://vdata.nikkei.com/newsgraphics/land-value-map_analysis2022-k6r8kj/

上で作成した夜間光テーブルと土地価格テーブルから、対象地域の2014年から2021年までの夜間光と土地価格の値を持ったテーブルを作成する関数を定義します。

def get_chart(code):
 years = [i for i in range(2014, 2022)]
 data = pd.DataFrame(columns = ['NL', 'LP'], index=years)
 data_nl = df_nl.loc[df_pnl['code']==code].drop('code', axis=1).drop('muni', axis=1)
 data_lp = df_lp.loc[df_lp['code']==code].drop('code', axis=1)
 for i in range(len(years)):
   data.at[years[i], 'NL'] = data_nl.iat[0, i]
   data.at[years[i], 'LP'] = data_lp.iat[0, i]
 return data

まずは千代田区のグラフを表示してみます。

data_chiyoda = get_chart(code='13101')

fig, ax1 = plt.subplots()
ax2 = ax1.twinx()
ax1.plot(data_chiyoda["NL"], color=cm.Set1.colors[1], label="Night Light")
ax2.plot(data_chiyoda["LP"], color=cm.Set1.colors[0], label="Land Price")
plt.title("Night Light and Land Price in Chiyoda", fontsize=15)
handler1, label1 = ax1.get_legend_handles_labels()
handler2, label2 = ax2.get_legend_handles_labels()
ax1.legend(handler1 + handler2, label1 + label2, loc=2, borderaxespad=0.)

変化率には相関がありませんでしたが、確かに推移で見ても異なる推移となっているようです。

次に稲城市のグラフを表示してみます。

data_inagi = get_chart(code='13225')
fig, ax1 = plt.subplots()
ax2 = ax1.twinx()
ax1.plot(data_inagi["NL"], color=cm.Set1.colors[1], label="Night Light")
ax2.plot(data_inagi["LP"], color=cm.Set1.colors[0], label="Land Price")
plt.title("Night Light and Land Price in Inagi", fontsize=15)
handler1, label1 = ax1.get_legend_handles_labels()
handler2, label2 = ax2.get_legend_handles_labels()
ax1.legend(handler1 + handler2, label1 + label2, loc=2, borderaxespad=0.)

こちらも推移を見てみると、やはり異なる推移となっています。

しかし、全体的な傾向としては夜間光も土地価格も右肩上がりになっていることから、変化率に相関はなくても、長期的に見れば夜間光が大きく増加した場所は土地価格も増加するというような予測は立てられそうです。

７．夜間光と土地価格の変化率の相関をグループ分けして調べてみる

上では夜間光と土地価格の変化率には相関がなさそうといった結果でした。次に変化率についてもう少し詳細に調査してみます。

ここで二つの仮説を立てて検証してみます。

一つ目は、郊外のような夜間光の値が小さいところでは変化率にも相関があるという説です。東京などの都心部では少々夜間光が変化したところで土地価格への影響は少なそうですが、郊外は開発が進めば夜間光も土地価格も上昇しそうですよね。そこで夜間光の値が高いグループと低いグループに分けて相関係数を調査してみます。

二つ目は、人口密度が低いところでは変化率にも相関があるという説です。人口密度が高い地域はマンションなどの集合住宅が多いと考えられます。人口が密集してすでに十分に明るい中で夜間光が変化したとしても、その変化を捉えるのは難しい可能性があります。人口密度が低い場所はマンションが少ないと考えられるため、住宅が増える＝人口が増える＝土地価格が高くなる（＝税収が増えるため公共投資できるようになる）＝夜間光が明るくなる、という変化を捉えやすいのではないかと考えました。

・夜間光の値でグループ分け
まずは夜間光の値でグループ分けをします。今回は2021年の夜間光の平均値を用いることにします。

まず、夜間光テーブル(df_nl)から市区町村コードと2021年の夜間光の平均値を抜き出したテーブルを作成し、2021年の夜間光の平均値でソートします。

loc_df2 = df_nl.loc[:, ['code','2021-NL']]

loc_df2['code'] = loc_df2['code'].astype(str)
loc_df2 = loc_df2[loc_df2['code'].isin(code_list)]
sorted_loc_df2 = loc_df2.sort_values('2021-NL').reset_index(drop=True)
sorted_loc_df2

次にグループ分けをします。夜間光の平均値が小さい順にソートされているので、前半の半分を低いグループ(low_group)、後半の半分を高いグループ(high_group)のリストに加えます。

high_group2 = []
low_group2 = []
for i in range(len(sorted_loc_df2)):
 if i < len(df_pnl_plp)/2:
   low_group２.append(sorted_loc_df2.at[i, 'code'])
 else:
   high_group２.append(sorted_loc_df2.at[i, 'code'])

上でグループ分けしたリストを基に、変化率のテーブルをグループ分けします。

low_group_df2 = df_pnl_plp[df_pnl_plp['code'].isin(low_group2)]
high_group_df2 = df_pnl_plp[df_pnl_plp['code'].isin(high_group2)]

夜間光の平均値が低いグループの相関係数を算出します。

low_group_df2.corr()

夜間光の平均値が高いグループの相関係数を算出します。

high_group_df2.corr()

図７　2014年-2015年の夜間光と土地価格の変化率

結果を見てみると、同じ年での変化率も異なる年の変化率でも相関はなさそうです。

・人口密度でグループ分け
次に人口密度でグループ分けをします。市区町村ごとの人口密度のデータは、e-statと呼ばれる日本の統計データがまとめられたサイトからダウンロードすることができます。
https://www.e-stat.go.jp/stat-search/database?page=1&query=%E4%BA%BA%E5%8F%A3%E5%AF%86%E5%BA%A6
&layout=dataset&statdisp_id=0003445099

ダウンロードした人口密度をpandasで読み込み、市区町村コードと人口密度を持ったテーブルを作成します。

ダウンロードしたCSVには都道府県ごとの合計や合併した市区町村の合併前の旧名など不要な行があるため、今回解析に使用する市区町村コードの行のみを抽出します。

その後、人口密度でソートします。

fp_pop_den = '人口密度のcsvのファイルパス'

df = pd.read_csv(fp_pop_den, encoding='SHIFT-JIS', skiprows=1)
loc_df = df.loc[:,['全国，都道府県，市区町村（2000年市区町村含む） コード','人口密度【1km2当たり】']]
loc_df.rename(columns={'全国，都道府県，市区町村（2000年市区町村含む） コード': 'code','人口密度【1km2当たり】': 'pop_den'}, inplace=True)
loc_df

code_list = []
for i in range(len(codes)):
 code_list.append(codes[i])

loc_df['code'] = loc_df['code'].astype(str)
loc_df['pop_den'] = loc_df['pop_den'].str.replace(',', '').astype(float)
loc_df = loc_df[loc_df['code'].isin(code_list)]
sorted_loc_df = loc_df.sort_values('pop_den').reset_index(drop=True)
sorted_loc_df

人口密度の小さい順に並んだテーブルを作成できました。

夜間光でグループ分けを行った時と同じようにグループ分けをして、変化率のテーブルを作成します。

high_group = []
low_group = []
for i in range(len(sorted_loc_df)):
 if i < len(df_pnl_plp)/2:
   low_group.append(sorted_loc_df.at[i, 'code'])
 else:
   high_group.append(sorted_loc_df.at[i, 'code'])

low_group_df2 = df_pnl_plp[df_pnl_plp['code'].isin(low_group2)]
high_group_df2 = df_pnl_plp[df_pnl_plp['code'].isin(high_group2)]

人口密度が小さいグループの相関係数を算出します。

low_group_df.corr()

人口密度が大きいグループの相関係数を算出します。

high_group_df.corr()

図８　2014年-2015年の夜間光と土地価格の変化率

結果を見てみると、人口密度の低いグループでも高いグループでも変化率に相関はなさそうです。また、同じ年同士だけでなく異なる年同士で見ても相関は見られません。

８．まとめ

本記事では全国の市区町村を対象に夜間光データと土地価格の相関を調査しました。「年間の夜間光の平均値と土地価格の相関」と「夜間光の平均値と土地価格の変化率の相関」の二つの項目について調査を行いました。

「年間の夜間光の平均値と土地価格」には相関がありそうでした。明るい場所は土地価格が高く、暗い場所は土地価格が低いという結果が得られました。

「夜間光の平均値と土地価格の変化率」には相関がなさそうでした。試しに夜間光の明るいグループと暗いグループ、人口密度の高いグループと低いグループに分けて相関を調査してみましたが、相関は見られませんでした。変化率に相関があれば、夜間光の変化率から土地価格の推移を予測でき、大儲け？と思っていたので、少し残念です(笑)

しかしながら、夜間光データを用いて土地価格を予測するアルゴリズム開発のコンペが開催されていたりと、もっと高度な処理をすれば面白い関係性がわかるかもしれません。
https://prtimes.jp/main/html/rd/p/000000007.000065306.html

土地価格以外にもたくさんのデータが公開されています。今回は土地価格との相関を調査してみましたが、今後その他の様々なデータとも相関を調査して記事にできたらと思っています。みなさんも夜間光と相関のある面白いデータを探ってみてはいかがでしょうか。

関連記事