【QGIS+Python】データを使って「駅の都会度論争」に決着をつける

背景

「〇〇駅と××駅とではどっちが都会か」という「駅の都会度論争」を見聞きすることはしばしばありますが、論理的に決着がついている例を見たことはほぼありません (たぶん)。
というわけで、データを使ってこの論争に決着をつけることを試みようと思いました。

前提

本記事では「都会度＝駅周辺の卸売業・小売業事業所の数」(以下事業所数) と定義して分析します。(この定義の妥当性については深く検討していません。何らかの指標に基づいて集計・分析することが目的です。)

使用したデータ

2016年「経済センサスー活動調査」の「産業 (大分類) 別事業所数及び従業者数」の4次メッシュ (500mメッシュ) データ (以下事業所数データ) を使用しました。e-Statに公開されています。境界データ (メッシュがどこで区切れるかを示すデータ) もここから取っています。
鉄道駅データは国土数値情報から取りました。

手順

①データのダウンロード

先述のとおり、e-Statから事業所数データと境界データをダウンロードしました。境界データはダウンロードの際に測地系やデータ形式を選択できますが、今回は「世界測地系平面直角座標系・Shapefile」を選択しました。また、対象地域は首都圏とし、メッシュ番号5239・5339・5340を選択しました (各メッシュの位置はこちら参照)。
ダウンロードしたzipファイルを解凍するとファイル形式が.txtだったので、拡張子を.csvに変更しました。
また、国土数値情報から鉄道駅のデータをダウンロードしました。

②データの下処理・味見

QGISを起動し、新規プロジェクトを立ち上げ、データを全て読み込み、各データに以下の下処理を施しました。

• 事業所データのデータ型を文字列型から整数型に変換
• 鉄道駅データの座標系を投影座標系に変換 (後のバッファ演算の際に都合がよいため)
• 境界データと事業所データとをテーブル結合し、境界データに事業所数のカラムを追加
• 3つのメッシュデータ (5239・5339・5340) を結合

以上の下処理を実施したうえで、データの味見をしました。下の写真は、メッシュごとの事業所数を色分けしたものです。味見が目的なので凡例等は省いていますが、黒い長方形が駅で、メッシュの色が濃いほど事業所数が多いです。新宿・渋谷・池袋の副都心が圧倒的かと思いきや、山手線の東側もかなり色が濃くなっていますね。

③Pythonで駅ごとに集計する

さて、ここから駅ごとの事業所数を集計していきます。やや厳密さには欠けますが、「駅Aの事業所数」を以下のように定義します。

(駅Aの事業所数)\equiv(駅Aからl[m]以内に位置するメッシュ内の事業所数の総計)

なお、駅Aからl[m]以内の少なくとも1つの地点PがメッシュBに含まれるとき、「メッシュBは駅Aからl[m]以内に位置する」と言うこととします。

以上の定義のもとでやるべきこととしては、以下が挙げられます。

1. 各駅からl[m]以内の地点を列挙する
   →ベクタの空間演算ツールである「バッファ」で実現可能
2. 列挙した地点をもとに、各駅からl[m]以内のメッシュを列挙する
   →ベクタの調査ツールである「場所による選択」で実現可能
3. 列挙したメッシュ内の事業所数を総和を計算する
   →「フィールド計算機」で実現可能？

上記の処理は集計対象の駅が少なければ手作業でもできますが、なんとなくスマートではないですよね。今回はPython (PyQGIS) を利用し、集計したい駅を列挙したCSVファイルを入力として、駅ごとの事業所数をCSVファイルとして出力するプログラムを以下のように書いてみました。プログラムはGitHubにも公開しています。

station_power.py

import processing
import pandas as pd
from qgis.core import QgsProject, QgsFeatureRequest, QgsVectorFileWriter, NULL

dir = '~/Desktop/station_power/'
input_fname = 'data/target_stations.csv'
output_fname = 'data/result.csv'
l = 500

# データの読み込み
df = pd.read_csv(dir + input_fname, encoding='cp932')
df['number_of_stores'] = 0

# レイヤーの選択
station_layer = QgsProject.instance().mapLayersByName('N02-20_Station_XY')[0]
mesh_layer = QgsProject.instance().mapLayersByName('MESH')[0]

for i in range(len(df)):
    # 対象駅名
    target = df['station_name'].iloc[i]
    
    # 対象駅を別レイヤとして出力
    station_layer.selectByExpression(('"N02_005"=\'' + target + '\''))
    target_station_layer = station_layer.materialize(QgsFeatureRequest().setFilterFids(station_layer.selectedFeatureIds()))
    QgsVectorFileWriter.writeAsVectorFormat(target_station_layer,dir + 'data/rail/target_station.shp', 'cp932', driverName='ESRI Shapefile')

    # 対象駅をl[m]でバッファリング
    buffer_output = processing.run("native:buffer", {'INPUT':target_station_layer,'DISTANCE':l,
    'SEGMENTS':5,'END_CAP_STYLE':0,'JOIN_STYLE':0,'MITER_LIMIT':2,'DISSOLVE':False,'OUTPUT':dir + 'data/rail/buffer_output.shp'})

    # バッファと交差しているメッシュの選択
    processing.run("native:selectbylocation", {'INPUT':dir + 'data/mesh/MESH.shp','PREDICATE':[0],
    'INTERSECT':dir + 'data/rail/buffer_output.shp','METHOD':0})

    # 選択中のメッシュ内の事業所数の総和を計算
    num_stores = 0
    for mesh_feature in mesh_layer.selectedFeatures():
        if mesh_feature['retail'] != NULL:
            num_stores += mesh_feature['retail']
    df['number_of_stores'].iloc[i] = num_stores


df.to_csv(dir + output_fname, encoding='cp932')

コードの詳しい説明はしませんが、入力データのCSVファイルをPandasで読み込み、入力データ記載の各駅について繰り返し、上記の1.~3.の処理をしました。1.と2.はがっつりPyQGISを使い、3.は普通にPythonに足し算させています。

結果

l = 500[m]のときの結果をグラフ化したものを以下にご共有します (字が小さいので適宜拡大をお願いします)。集計対象の駅は完全に私の独断で選定しました。

ここから私の考察というより感想を書きますが、あくまで個人的な感想なので悪しからず。
トップは新宿・渋谷・池袋になると思い込んでいましたが、東京・上野・新橋あたりもかなり強いのですね。東京駅は駅ナカが大きい&駅自体が大きいことが影響したようにも思います。
続いて横浜・大宮と、神奈川県・埼玉県の代表駅が続きます。これは納得です。(千葉…。)
そのあとの吉祥寺・錦糸町はかなり栄えている印象はありますが、ここまで上位に食い込んでくるのは意外でした。そして千葉・川崎と続く。正直川崎>千葉だと思っていましたが、さすがに県庁所在地の意地を見せてくれました。
その次が蒲田というのも意外でしたが、日暮里にはかなり驚きました。駅周辺そんなに何かありましたっけ、、
六本木・品川・北千住を挟みつつ立川・町田・八王子の多摩地域代表駅がランクイン。このあたりは納得です (品川はもっと下だと思っていましたが)。自由が丘は降りたことがないため不明。大森はそんなに大きかったでしたっけ、、。
これ以降は郊外駅が名を連ねますが、あまり意外なところがなかったため省略。特筆するとすれば、船橋・藤沢あたりが思ったより小さいなーということくらいですかね。

「この駅も調べてほしい」という希望や、結果に対する考察、手法等に関するツッコミ等ありましたら、ぜひコメントくださると嬉しいです！！

(2023/11/25追記)
対象駅を追加した結果を以下にご共有します。主に首都圏以外の駅を追加しています。

大阪が第4位・難波が第6位とさすがの結果。博多が第7位、三ノ宮が第10位にランクインするなど、都市規模的にはやや劣るものの求心力の高い駅が予想以上に上位に食い込んできました。名古屋市からは栄が第12位、名古屋が第15位。やはり名古屋市の中心は栄なのですね。仙台は渋谷に続き第14位。仙台駅前は初めて行ったときにとても都会で驚いたことを覚えています。札幌は横浜に続いて第19位。もう少し上かと思っていたのでこれは意外でした。札幌を境に事業所数が大幅に落ち、京都河原町・天王寺・小倉・広島などがランクイン。このあたりはおおむね予想どおりの結果ですね。

ちなみに、各市区町村の中で最も順位が高い駅で比較すると、「東京23区>大阪市>福岡市>神戸市>名古屋市>仙台市>横浜市>札幌市>さいたま市>京都市>武蔵野市>北九州市>広島市>千葉市>川崎市>立川市>…」という順番となります。面白い結果になりましたね。

参考資料

国勢調査結果をQGISで表示する手順を画像で解説：無料でGISを使ってみる
QGISの操作「QGISに取り込んだデータの座標参照系を変換する方法」
pyQGISでできることのまとめ（自分用）
QGIS + Python その1
QGIS + Python その2