全国の人流オープンデータを触る(その２)

import geopandas as gpd
import pandas as pd
import zipfile

# CRSのセット
crs = 'EPSG:6677'

# 世界遺産のデータを読み込み
gdf1 = gpd.read_file('shp/A34a-180316.geojson')

# 日光を抽出
heritage = gdf1[gdf1['A34a_003']=='日光の社寺']
heritage = heritage.to_crs(crs) # 座標系の変換
heritage.to_file('output/heritage.shp', encoding='cp932') # 日光のシェープ保存

# 栃木県の1kmメッシュデータを読み込み
gdf2 = gpd.read_file('shp/09tochigi1km.geojson')
gdf2 = gdf2.to_crs(crs) # 座標系の変換

# 日光の社寺付近にかかるメッシュ番号を抽出
mesh_list = []
for index1, h in heritage.iterrows():
    for index2, m in gdf2.iterrows():      
        if m['geometry'].intersects(h['geometry']): 
            mesh_list.append(m[0])

# メッシュデータの抽出
mesh = gdf2[gdf2['code'].isin(mesh_list)]
mesh.to_file('output/mesh.shp', encoding='cp932') # メッシュ4カ所のシェープ保存

# 人流データの読み込みと結合
pieces = []
file = 'pop/monthly_mdp_mesh1km_09.zip'

with zipfile.ZipFile(file) as f:
    pref = f.namelist()[0]
    for y in range(2019, 2022):
        for m in range(1, 13):
            with f.open('{}{}/{}/monthly_mdp_mesh1km.csv.zip'.format(pref, str(y), str(m).zfill(2))) as ff:
                with zipfile.ZipFile(ff).open('monthly_mdp_mesh1km.csv') as fff:
                    frame = pd.read_csv(fff)
                    frame['code'] = frame['mesh1kmid'].astype(str) # 人流データはコードがint型なので
                    frame = frame[frame['code'].isin(mesh_list)]
                    pieces.append(frame)

pop = pd.concat(pieces, ignore_index=None)
pop = pop.drop(['mesh1kmid', 'prefcode', 'citycode'], axis=1)
pop.to_csv('output/pop.csv', encoding='cp932', index=None) # 人流データcsv保存

import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# 人口データ読み込み
pop = pd.read_csv('output/pop.csv', dtype={'code': 'str'})

# ＜単純に4カ所のメッシュのみで集計＞
# グラフの準備
pop1 = pop.groupby(['year', 'month', 'dayflag', 'timezone']).sum().reset_index()
print(pop1['population'].max()) # グラフの目盛りの目安に
order = {1:[1, 0, '平日　昼間'], 2: [1, 1, '平日　深夜'], 3: [0, 0, '休日　昼間'], 4: [0, 1, '休日　深夜']} # グラフの配置と集計フラグを整理
h = 12000 # 目盛り

# グラフ作成
sns.set(font='Hiragino Maru Gothic ProN')
fig = plt.figure(figsize=(20,10))
plt.subplots_adjust(wspace=0.2, hspace=0.3)
for k, v in order.items():
    df = pop1[(pop1['dayflag']==v[0]) & (pop1['timezone']==v[1])]
    ax = fig.add_subplot(2, 2, k)
    sns.barplot(data=df, x='month', y='population', hue='year', ci=None)
    plt.ylim(0, h)
    plt.title(v[2])

fig.savefig('output/graph1.png')


# ＜面積按分による集計＞
# 世界遺産（日光の社寺）のシェープファイルを読み込み
heritage = gpd.read_file('output/heritage.shp', encoding='cp932')

# メッシュのシェープファイルを読み込み
mesh = gpd.read_file('output/mesh.shp', encoding='cp932')
mesh['area_all'] = mesh.area # 面積を算出（分母）

# オーバーレイ
over = gpd.overlay(mesh, heritage, how='intersection')
over['area_cut'] = over.area #  面積を算出（分子）
over['rate'] = over['area_cut'] / over['area_all'] # 比率

# 人口に比率をかけ、按分人口を算出
pop2 = pd.merge(pop, over[['code', 'rate']], on='code')
pop2['population'] = pop2['population'] * pop2['rate']

# グラフの準備
pop2 = pop2.groupby(['year', 'month', 'dayflag', 'timezone']).sum().reset_index()
print(pop2['population'].max()) # グラフの目盛りの目安に
h = 1500 # 目盛り

# グラフ作成
fig = plt.figure(figsize=(20,10))
plt.subplots_adjust(wspace=0.2, hspace=0.3)
for k, v in order.items():
    df = pop2[(pop2['dayflag']==v[0]) & (pop2['timezone']==v[1])]
    ax = fig.add_subplot(2, 2, k)
    sns.barplot(data=df, x='month', y='population', hue='year', ci=None)
    plt.ylim(0, h)
    plt.title(v[2])

fig.savefig('output/graph2.png')