人流データでみる東京

こんにちは, 22 新卒のdaishinというものです.

ひょんなことから, 来年の 4 月より東京で働くことになってしまったのですが, 生まれてこのかた都会で暮らした経験が皆無な田舎者なため, 大都会での暮らしに戦々恐々としながら残りの大学生活を過ごしています.
できれば, 少しでも東京に詳しくなり, 華々しいシティボーイデビューを飾りたいのですが, 如何せん, 時間がなかったり, お金がなかったりで, 東京を肌で感じる機会になかなか恵まれませんでした.
そこで, 一般に公開されているオープンデータなどを活用し, 東京について詳しくなってやろうと思い？？, 手始めに人の流れから見てみることにしました.
これで東京のナウでホットなスポットを知ることができますね.

---

さて, ここからが本編です.

はじめに

この記事は DeNA 21 新卒 ×22 新卒内定者 Advent Calendar 2021 の 20 日目の記事となっています.

この記事で, 冒頭でもあったように, オープンデータとして公開されている人流データを簡単に可視化してみようという内容になっています（複雑なデータ処理などはしていません）.
使用言語は python, データの処理・可視化のために geopandas というライブラリを利用しています.
geopandas や GIS データの詳細な紹介を行う記事ではないため, その辺はご了承ください.

使用データ

全国の人流オープンデータ

国土交通省がG 空間情報センターにてオープンデータとして配布している 2019・2020 年の全国の人流データを利用しました.
https://www.geospatial.jp/ckan/dataset/mlit-1km-fromto
携帯電話端末等の位置情報データを元に集計されており, 市町村単位, 1km メッシュ単位のデータとして提供されています.

本記事では, 東京都のデータのみをダウンロードし, 以下の 2 つのファイルを利用しています.

• monthly_mdp_mesh1km：1km メッシュ単位の換算人口値が格納されたデータ
• attribute：1km メッシュデータにおけるメッシュ ID の座標を示すデータ

データ定義書によると, monthly_mdp_mesh1km は以下のような属性を持っています.

• mesh1kmid : メッシュコード
• prefcode : 都道府県コード
• citycode : 市町村コード
• year : 年（2019 or 2020）
• month : 月
• dayflag : 平休日（"0":休日, "1":平日, "2":全日）
• timezone :　時間帯（ "0":昼, "1":深夜, "2":終日）
• population : 平均滞在人口（10 以上）

東京都の境界データ

各府省の様々な統計データをまとめているポータルサイトであるe-statから東京都の境界データを入手しました.

2015 年の国勢調査の shapefile 形式のものをダウンロードしました.
https://www.e-stat.go.jp/gis/statmap-search?page=1&type=2&aggregateUnitForBoundary=A&toukeiCode=00200521&toukeiYear=2015&serveyId=A002005212015&coordsys=1&format=shape&datum=2000

GIS について

GIS とは, Geographic Information System（地理情報システム）の略で, 地理空間情報（地理的な位置情報を持ったデータ）の管理や加工, 可視化などの分析を可能にするための技術のことをいいます.

GIS データには様々なデータ形式が存在するのですが, 描画方法で大別すると, ラスターデータとベクターデータに分けることができます.
ラスターデータとは, 地図を格子状に分割し, その格子状の中にデータが入ったもの, ベクターデータは点・線・ポリゴンの 3 つの要素で表現されるデータになります.
先ほど境界データの説明で出てきた shapefile 形式というのはベクターデータを扱うためのデータ形式となっています.

データ処理

前提条件

本記事でのデータのディレクトリ構造は以下のようになっています.

data
├ monthly_mdp_mesh1km
│ └ tokyo
│   └ {year}_{month}.csv  # 1kmメッシュあたりの滞在人口データ
├ attribute               # 座標データ（csv形式）
└ tokyo_boundary          # 東京の境界データ（shapefile形式）

データの処理には python を採用しています.
また, 利用したライブラリは以下のとおりです

pandas : データ処理用のライブラリ
geopandas : gisデータを扱うためのライブラリ
shapely : ジオメトリを取り扱うためのライブラリ

ざっくりとデータを確認

データがどんな感じなのかを確認するために, 2020 年 1 月のものを可視化してみます.

from pathlib import Path
import pandas as pd

mesh1km_dir = Path('monthly_mdp_mesh1km/tokyo')
attribute_dir = Path('attribute')

attr_2020_df = pd.read_csv(attribute_dir/'attribute_mesh1km_2020.csv')
mesh1km_2020_01_df = pd.read_csv(mesh1km_dir/'2020_01.csv')

座標データが欲しいため, mesh1km_2020_01_dfとattr_2020_dfを結合し, geopandas の GeoDataFrame に変換します.
ここでは, geopandas のpoints_from_xy関数を用いて, 座標データから Point 形式の geometry に変換しています.
また, 元のデータには小笠原諸島などの離島のデータも含まれていましたが, 今回は離島のデータは排除しておきました.

attr_cols = [col for col in attr_2020_df.columns if col not in ['prefcode', 'citycode']]
sample_df = pd.merge(mesh1km_2020_01_df, attr_2020_df[attr_cols], how='left', on='mesh1kmid')
sample_gdf = gpd.GeoDataFrame(sample_df, geometry=gpd.points_from_xy(sample_df.lon_center, sample_df.lat_center))
sample_gdf = sample_gdf.query('35 < lat_center < 36 & 139 < lon_center < 140')

plot メソッドを使って可視化してみます.

sample_gdf.plot(column='population', cmap='jet', figsize=(25, 15))

なんとなくですが東京の形になってますね.

前処理

滞在人口のデータがどんなものかざっくりと可視化できたので, 全てのデータを読み込んで東京の人の流れを把握していきたいと思います.

まずは, 全てのデータを読み込み, 位置情報のデータを結合します.

df_list = []
for p in mesh1km_dir.glob('2020_*.csv'):
    df = pd.read_csv(p)
    df_list.append(df)

mesh1km_df = pd.concat(df_list)

attr_cols = [col for col in attr_2020_df.columns if col not in ['prefcode', 'citycode']]
mesh1km_df = mesh1km_df.merge(attr_2020_df[attr_cols], how='left', on='mesh1kmid')

続いて不要な情報を除いておきます.

mesh1km_df = mesh1km_df.query('35 < lat_center < 36 & 139 < lon_center < 140')

drop_cols = ['mesh1kmid', 'prefcode', 'year', 'lon_center', 'lat_center']
mesh1km_df.drop(drop_cols, axis=1, inplace=True)

また, 平休日・時間帯でデータを分けておきます.

mesh1km_df_22 = mesh1km_df.query('dayflag == 2 & timezone == 2') # 終日・全日
mesh1km_df_20 = mesh1km_df.query('dayflag == 2 & timezone == 0') # 終日・昼間
mesh1km_df_21 = mesh1km_df.query('dayflag == 2 & timezone == 1') # 終日・深夜

mesh1km_df_02 = mesh1km_df.query('dayflag == 0 & timezone == 2') # 平日・全日
mesh1km_df_00 = mesh1km_df.query('dayflag == 0 & timezone == 0') # 平日・昼間
mesh1km_df_01 = mesh1km_df.query('dayflag == 0 & timezone == 1') # 平日・深夜

mesh1km_df_12 = mesh1km_df.query('dayflag == 1 & timezone == 2') # 休日・全日
mesh1km_df_10 = mesh1km_df.query('dayflag == 1 & timezone == 0') # 休日・昼間
mesh1km_df_11 = mesh1km_df.query('dayflag == 1 & timezone == 1') # 休日・深夜

geometry カラムの追加

先ほどはデータを Point で可視化しましたが, せっかく 1km メッシュ（1km の格子状）としてデータが与えられているので, 格子状でデータを可視化して行きたいと思います.
そのため, shapely ライブラリを用いて, csv で与えられた 4 角の座標から Polygon を作成し, geometry カラムとして付与します.

def get_mesh_point(df: pd.DataFrame):
    lat_min = df['lat_min']
    lat_max = df['lat_max']
    lon_min = df['lon_min']
    lon_max = df['lon_max']

    point1 = (lon_min, lat_min)
    point2 = (lon_min, lat_max)
    point3 = (lon_max, lat_max)
    point4 = (lon_max, lat_min)

    return Polygon([point1, point2, point3, point4])

mesh1km_df_22['geometry'] = mesh1km_df_22.apply(get_mesh_point, axis=1)

geopandas の GeoDataFramen 変換します.
この際, 地図の座標系を緯度経度で表現される地理座標系にセットしておきます.
（ここでは地図の座標系についての詳しい説明は割愛します. ）

mesh1km_gdf = gpd.GeoDataFrame(mesh1km_df_22)
mesh1km_gdf.set_crs(epsg=4612, inplace=True)

データの可視化

plot 関数を使って Polygon を可視化してみます.
今の状態だと, 各 Polygon ごとに 12 ヶ月分の行が存在しているため, 今回は全て合計して表示したいと思います.

all_month = mesh1km_gdf[['population', 'geometry']].groupby(by='geometry', sort=False, as_index=False).sum()
all_month.plot(column='population', cmap='jet', edgecolor='w', alpha=0.5, figsize=(15,9))

これだけだと少し味気ないので, 東京の境界データと重ねてみたいと思います.

境界データの shapefile を geopandas で読み込みます.

shape_file = Path('data/tokyo_boundary/h27ka13.shp')
tokyo_gdf = gpd.read_file(shape_file)
tokyo_gdf = tokyo_gdf.query('35 < Y_CODE < 36 & 139 < X_CODE < 140')['geometry']

読み込めたら, 人流データと重ね合わせて表示してみます.

ax = tokyo_gdf.plot(figsize=(15, 9), color='none')
all_month.plot(ax=ax, column='population', cmap='jet', edgecolor='w', alpha=0.5)

少し見にくいですが, なんとなく人が集まっている地域が確認できますね.

他の時間帯のものも確認してみます.
月ごとの違いがわかるように gif アニメーションでまとめてみました.

平日・全日

休日・全日

geocoding

住所やランドマークの名称などを座標に変換することを geocoding, その逆の作業を逆 geocoding と言います.

今回, 地図上に人口のデータを可視化しましたが, 如何せん私は東京の土地勘がないため, どこがどういう場所なのかわかりません.
そのため, geocoding で人が集まっている場所の住所を検索してみたいと思います.

geopandas ではtools.geocodeメソッド・tools.reverse_geocodeメソッドによって geocoding, 逆 geocoding を行うことが可能となっております.
このメソッドでは内部的に geopy というライブラリを用いて外部 API を利用しています.

geocoding に利用するデータを用意し, geocoding を実行します.
今回は休日の夜間の滞在人口が上位 5 位までの場所を検索しようと思います.
geocoding メソッドでは引数で Point 形式のリストを受け取るため, 最初に利用したsample_gdfを使っています.

sample_gdf_02 = sample_gdf.query('dayflag == 0 & timezone == 2')
search_gdf = sample_gdf_02.sort_values(by='population', ascending=False).head(5)

gdf_dam_rgeo = gpd.tools.reverse_geocode(search_gdf['geometry'], provider='nominatim', user_agent='test')

     geometry                    address
443  POINT (139.70617 35.69590)  ホテル・グランシャリオ, 明治通り, 新宿五丁目, 新宿区, 東京都, 160-8430, 日本
164  POINT (139.76852 35.67928)  東京駅, 外堀通り, 銀座一丁目, 八重洲二丁目, 中央区, 東京都, 100-6633, 日本
3546 POINT (139.70630 35.72920)  西池袋公園, 441, 南池袋一丁目, 西池袋三丁目, 豊島区, 東京都, 171-0021...
416  POINT (139.70955 35.68507)  新宿御苑, 明治通り, 内藤町, 渋谷区, 東京都, 160-0014, 日本
3555  POINT (139.71955 35.72883) サンシャインシティ, 3, 文化会館3階高速荷捌所出口, 南池袋二丁目, 東池袋三丁目, 豊...

しっかりと結果が返ってきました.

まとめ

本記事では python を利用して, 2020 年の人流データを可視化してみました.
軽く触れた程度にはなってしまいましたが, 様々な統計情報などと組み合わせることでかなり有用なことができるのではないかと思いました.

他にも様々な地理データや統計データなどが公開されているので, それらを利用して東京についてより詳しくなり, シティボーイになれるようにように研鑽して行きたいと思います.