こんにちは。つくばに住む研究者です。

今回は茨城県の各市区町村の昼夜間人口比率を使いながら、地図を使ったデータ分析の初歩をやってみたいと思います。

開発環境は今回もColabを使います。利用するデータは国勢調査の結果に加えて、市区町村の地図上での境界を示すポリゴンデータを使います。境界データは国土交通省の整備する国土数値情報ダウンロードサイトの’行政区域データ’から手に入ります。今回は茨城県の令和5年のデータをダウンロードしました。前回までと同じように、解凍したファイルを国勢調査のデータを置いたGoogle Driveのフォルダにおきました。

まずは国勢調査のデータを茨城県の市区町村に絞り込んで、昼夜間人口比率を取り出して見て見ましょう。

import pandas as pd
import numpy as  np
from google.colab import drive

drive.mount('/content/drive')
dir_path = '/content/drive/MyDrive/Colab Notebooks/Data/2023_06/'
df_all = pd.read_excel(dir_path+'major_results_2020.xlsx')

#全国レベル,政令指定都市下の区のデータを除去
df_crop = df_all[df_all['市などの別（地域識別コード）']!='a']

#茨城県の自治体だけに絞りこみ
df_crop = df_all[df_all['都道府県名']=='08_茨城県']

#不完全なデータの除去
df_crop = df_crop[df_crop['昼夜間\n人口比率']!='-']

#計算できるようにfloat型に変換する
df_crop[['昼夜間\n人口比率']] = df_crop[['昼夜間\n人口比率']].astype(float)

#市区町村名の前の自治体コードを除去
df_crop['都道府県・市区町村名'] = df_crop['都道府県・市区町村名'].str.split('_', expand=True)[1]

df_crop[['都道府県名', '都道府県・市区町村名','昼夜間\n人口比率']].head(10)

茨城県全体の昼夜間人口比率は97.6%です。県全体では東京などへの通勤や通学などにより昼間は人口が減るようです。各市区町村の比率を地図上で色分けしながら見ていきましょう。

まずはColabで地図データを扱うためのライブラリをインポートします。Geopandasという地図データに特化したライブラリを使いますが、標準のColab環境にはインストールされていないため、pipを使ってインストールします。

!pip install --upgrade geopandas --q

インストールできたらインポートし、先ほど国土数値情報ダウンロードサイトから入手した'.shp'ファイルを読み込みます。

import shapely
import geopandas as gpd
import folium

gdf = gpd.read_file(dir_path+'/ibaraki_shapes/N03-23_08_230101.shp',encoding='shift-jis')

それでは茨城県の全ての市区町村について境界を書きます。

fig = Figure(width=1000, height=800)

# 地図の中心座標を設定する　ポリゴンの始点と中間点のX,Yの座標の平均値を中心に設定する
x = (gdf.iloc[idx].geometry.exterior.xy[1][0]+gdf.iloc[0].geometry.exterior.xy[1][int(len(gdf.iloc[0].geometry.exterior.xy[1])/2)])/2
y = (gdf.iloc[idx].geometry.exterior.xy[0][0]+gdf.iloc[0].geometry.exterior.xy[0][int(len(gdf.iloc[0].geometry.exterior.xy[0])/2)])/2
map = folium.Map(location=[x,y] ,zoom_start=11)

fig.add_child(map)
folium.GeoJson(gdf["geometry"]).add_to(map)
map

茨城県では南に行くほど市の面積が減る傾向があるようです。つくば市単独の境界を見て見ましょう。

idx = int(gdf[gdf['N03_004']=='つくば市'].index[0])

fig = Figure(width=800, height=600)

# 地図の中心座標を設定する　ポリゴンの始点と中間点のX,Yの座標の平均値を中心に設定する
x = (gdf.iloc[idx].geometry.exterior.xy[1][0]+gdf.iloc[idx].geometry.exterior.xy[1][int(len(gdf.iloc[idx].geometry.exterior.xy[1])/2)])/2
y = (gdf.iloc[idx].geometry.exterior.xy[0][0]+gdf.iloc[idx].geometry.exterior.xy[0][int(len(gdf.iloc[idx].geometry.exterior.xy[0])/2)])/2
map = folium.Map(location=[x,y] ,zoom_start=11)

fig.add_child(map)
folium.GeoJson(gdf.iloc[idx]["geometry"]).add_to(map)
map

つくば市の境界を地図上に描くことができました。格都市について、比率に基づいて色を塗っていきましょう。まずは昼夜間人口比率が高いほど濃い青色、比率が低いほど濃い赤色、100％に近くと白色になるようなカラーマップとしてcoolwarmをを選びました。

注目する数字の範囲に色を対応させます。今回は90~110の範囲の中で色を対応させていきます。90以下や110以上の場合は一番濃い赤色または一番濃い青色になります。

import matplotlib.cm as cm
import numpy as np

#色情報の基礎を作成
c = cm.coolwarm(range(0,200,10))
v= np.arange(90, 110, 1)

colormapはRGBAの4次元配列で与えられます。ポリゴンに色を塗る際は色をhex表記で表す必要があるため、RGBをHEXに変換する関数を作ります。

def rgb_to_hex(rgb):
    return '%02x%02x%02x' % rgb

ポリゴンデータのデータフレームに昼夜間人口の数字をコピーしておきます。

for i in range(0,len(gdf)):
  gdf.loc[i, '昼夜間人口比率'] = float(df_crop[df_crop['都道府県・市区町村名'] == gdf.iloc[i]['N03_004']]['昼夜間\n人口比率'])

昼夜間人口の数字に応じた色を新しい列に保存しておきます。通常、matplotlibを使った描写の塗る際に、予め色を作ってデータフレームに保存しておくことはありませんが、地図の描写につかっているfoliumではポリゴン毎に色を変える処理を実現するには少し難しい書き方をしなければならないため、今回はこの方法を選びました。

試しに素直なForループを使い、人口比率ごとに設定された色を使って全ての市区町村を描いて見ます。

fig = Figure(width=1000, height=800)

x = (gdf.iloc[idx].geometry.exterior.xy[1][0]+gdf.iloc[idx].geometry.exterior.xy[1][int(len(gdf.iloc[idx].geometry.exterior.xy[1])/2)])/2
y = (gdf.iloc[idx].geometry.exterior.xy[0][0]+gdf.iloc[idx].geometry.exterior.xy[0][int(len(gdf.iloc[idx].geometry.exterior.xy[0])/2)])/2
map = folium.Map(location=[x,y] ,zoom_start=11)

for i in range(0,len(gdf)):
  folium.GeoJson(gdf.iloc[i]["geometry"],
                style_function = lambda feature:{"fillColor": gdf.iloc[i]['色'],
                                                  "color": '#000000',
                                                  "weight": 1,
                                                  "fillOpacity": 0.5}).add_to(map)

map

全ての市区町村の色が同じになってしまいました。foliumでForループを使ってポリゴンごとのスタイルを変更するには、少し難解な書き方が必要のようです。今回は力技を使います。Forを使ってポリゴン毎の描写ができないのであれば、いっそのことForを使った場合に実行される全てのコードをテキストで用意してしまい、1つずつ実行させてしまいましょう。そのために、Forループを使って実行したいコードを生成します。

sentence = ''
for i in range(0,len(gdf)):
  sentence = sentence+'folium.GeoJson(gdf.iloc['+str(i)+'][\"geometry\"],style_function = lambda feature:{\"fillColor\": gdf.iloc['+str(i)+'][\'色\'],\"color\": \'#000000\',\"weight\": 1,\"fillOpacity\": 0.8}).add_to(map) \n'

こうして作られたsentence変数の中には、通常ならForループで実行されるであろう、約200行分のコードが入っています。pythonでは文字列で作られたコードをexec()で実行することができます。つまりexec(sentence)を実行すれば、200個分のポリゴンを描写するためのコードをForを使わず逐次処理で実行できます。

idx = 0
fig = Figure(width=800, height=600)
x = (gdf.iloc[aid].geometry.exterior.xy[1][0]+gdf.iloc[idx].geometry.exterior.xy[1][int(len(gdf.iloc[idx].geometry.exterior.xy[1])/2)])/2
y = (gdf.iloc[aid].geometry.exterior.xy[0][0]+gdf.iloc[idx].geometry.exterior.xy[0][int(len(gdf.iloc[idx].geometry.exterior.xy[0])/2)])/2
map = folium.Map(location=[x,y] ,zoom_start=9)
fig.add_child(map)

exec(sentence)

map

実行結果は次のようになります。

東京にアクセスのよい南西側の地域は総じて低い昼夜間人口比率です。通学や通勤で県外に移動する人口が多いことが窺えます。県内では、赤の地域が青の地域を取り囲むように形成されています。水戸市、日立市、土浦市、つくば市などの昼間人口が多い大規模な都市に、近隣都市から流入する構造が窺えます。もっとも西側に見える濃い青は五霞市で、昼夜間人口比率は驚異の158％でした。

茨城県の市区町村の昼夜間人口比率については、茨城県のホームページから確認できる統計データでもまとめられています（本記事とほとんど同一の結果となっています）。https://www.pref.ibaraki.jp/kikaku/tokei/fukyu/tokei/betsu/jinko/kokucho/kokucho2020-3/index.html

次回も地図を使った分析などを行う予定です。

それでは。