pgRoutingで札幌市内のすべての駅の到達圏解析を行う #PostgreSQL

はじめに

到達圏とは、ある地点から指定した距離や時間で到達可能な領域のことをいいます。
こちらのサービス（How far can I go?）のように、ある1地点の到達圏を可視化するサービスはありますが、今回はpgRoutingを使用して札幌市内にあるすべての駅からの到達圏を作成する方法について説明します。

目標

pgRoutingとは

pgRoutingとは、PostGIS/PostgreSQLにネットワーク解析機能を追加する拡張機能です。

到達圏解析以外にも、経路探索や巡回セールスマン問題など、ネットワーク解析に関するさまざまな機能が用意されています。

pgRouting

データの取得

今回の解析に必要なデータを取得します。作業用フォルダとして[pgrouting]フォルダを作成します。さらに、[pgrouting]直下に[data]フォルダを作成し、そこにダウンロードしたデータを解凍し格納しておきます。

鉄道データ（全国令和3年）

https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N02-v3_0.html
行政区域データ（北海道平成31年）

https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v2_3.html
道路データ

今回は、地理院ベクトルタイルの道路中心データを使用するため、QGISのプラグイン「GSI-VTDownloader」を使用して、札幌市内の駅がすべて含まれるように矩形範囲を設定しダウンロードします。後ほど、QGISの機能を使用してデータベースにインポートするため、ここでshpなどにエクスポートしておく必要はありません。

pgRoutingのインストール

Dockerを使用してpgRoutingをインストールします。まずは、作業用フォルダ[pgrouting]にdocker-compose.ymlを作成します。こちらのQiita記事を参考にさせていただきました。

version: "3"

services:
  pgrouting3:
    image: pgrouting/pgrouting:13-3.0-3.1.1
    ports:
      - "5432:5432"
    volumes:
      - db-data3:/var/lib/postgresql/data
    environment:
      - POSTGRES_USER=postgres
      - POSTGRES_DB=routing
      - POSTGRES_PASSWORD=postgres

volumes:
  db-data3:

コンテナを起動します。

docker-compose up -d

データベースの準備

psqlを使用してPostgreSQLに接続し、データベースにPostGISとpgRoutingの拡張機能を追加します。

psql --host=localhost --user=postgres postgres
\c routing
CREATE EXTENSION postgis;
CREATE EXTENSION pgrouting;

データベースへのインポート

GDAL/OGRを使用して、先ほどダウンロードしたデータをデータベースにインポートします。

cd data

# 駅データのインポート
ogr2ogr -lco GEOMETRY_NAME=geom -f "PostgreSQL" -a_srs "EPSG:4326" PG:"host=localhost user=postgres password=postgres dbname=routing" ./N02-21_GML/N02-21_Station.shp -nln station

# 行政区域データのインポート
ogr2ogr -lco GEOMETRY_NAME=geom -f "PostgreSQL" -a_srs "EPSG:4326" PG:"host=localhost user=postgres password= postgres dbname=routing" ./N03-190101_01_GML/N03-19_01_190101.shp -nln gyouseikuiki

道路データは、QGISのDBマネージャーを使用してデータベースへインポートします。

インポートでエラーが出た場合は、プロセシングツールのマルチパートをシングルパートに変換を実行してから最後実行するとインポートが実行可能かと思います。

なお、QGISとデータベースとの接続の設定については、下記のページを参考にしてください。

QGISでPostGISのデータを見てみよう｜PostGISを入れるところからやってみよう

データの前処理

ここからは、DB管理ツールを使用してクエリを実行していきます。

道路データと駅データには、札幌市外のデータも含まれているため札幌市のみのデータに加工し、このタイミングで平面直角座標系12系（EPSG:6680）に変換しています。

また、駅データについてはラインデータで作成されているためポイントデータに変換しておきます。

--空間インデックスを作成
CREATE INDEX ix_road ON road USING gist(geom);
CREATE INDEX ix_station ON station USING gist(geom);
CREATE INDEX ix_gyouseikuiki ON gyouseikuiki USING gist(geom);

--札幌市内の駅を抽出し、ラインの重心にポイント作成
select 
	a.n02_005 as station_name,
	st_centroid(ST_Transform(a.geom,6680)) as geom
into sapporo_station
from station as a
	left join (
		select *
		from gyouseikuiki
		where n03_003 = '札幌市')as b 
	on st_intersects(a.geom,b.geom)
where n03_001 is not null
;

--札幌市内の道路を抽出
select 
	row_number() over() as id,
	ST_Transform(a.geom,6680) as geom
into sapporo_road
from road as a
	left join (
		select *
		from gyouseikuiki
		where n03_003 = '札幌市')as b 
	on st_intersects(a.geom,b.geom)
where n03_001 is not null
;

sapporo_stationとsapporo_roadをQGISで表示した結果

ネットワーク・トポロジーの作成

道路データ（sapporo_road）から、ネットワーク・トポロジーを構築します。

道路データに対して、pgr_createTopology()を実行することにより、道路間の接続関係の情報が構築されます。具体的には、道路データのsourceとtargetカラムにラインの端点のノード番号が付与され、さらに<道路テーブル>_vertices_pgrという名前でノードテーブルが生成されています。

また、このタイミングで道路データのcostカラムにエッジを通過する時間を追加しておきます。今回は、分速80mを前提として計算しました。

--カラムを作成する
ALTER TABLE sapporo_road ADD "source" integer NULL;
ALTER TABLE sapporo_road ADD "target" integer NULL;
ALTER TABLE sapporo_road ADD "cost" double precision NULL;
ALTER TABLE sapporo_road ADD "reverse_cost" double precision NULL;

--コストとして、分速80m前提でエッジの通過時間を計算
UPDATE sapporo_road SET cost = st_length(geom) / 80;
UPDATE sapporo_road SET reverse_cost = cost;

--トポロジーデータ作成
SELECT pgr_createTopology('sapporo_road', 0.001, 'geom', 'id');

sapporo_road（エッジ）

sapporo_road_vertices_pgr（ノード）

sapporo_roadとsapporo_road_vertices_pgrを表示してみると、交差点（とタイルの境）にノードが作成されていることがわかります。

到達圏解析

到達圏解析を行う準備ができたので、まずは到達圏を作成する始点となる各駅の最近傍のノードを取得します。

その後、pgr_drivingDistance()を使用してノードから指定の時間内（今回は徒歩10分）に通過可能なエッジをすべて取得します。

--各駅から最近傍のノードIDを取得
select 
	a.station_name,
	c.id as node_id,
	c.the_geom as geom
into start_node
from 
	sapporo_station as a
cross join lateral(
	select *
	from sapporo_road_vertices_pgr as b
	order by a.geom <-> b.the_geom
	limit 1
	)as c
;

--start_nodeから指定したコストで通過可能なエッジを取得する
with driving_distance as(
select a.*,b.geom
from pgr_drivingDistance
	('select id,source,target,cost,reverse_cost from sapporo_road',
	(SELECT array_agg(node_id) FROM start_node), --スタート地点のノードIDを配列で渡す
	10, --到達圏を作成するコスト
	false)as a
	left join sapporo_road as b
	on a.edge = b.id
) 

select *
into driving_distance_10min
from driving_distance;

driving_distanceを表示してみると、各駅を中心としたエッジが取得されていることがわかります。

最後に、すべてのエッジを囲む凹となるポリゴン(到達圏)を作成します。

pgRoutingにはpgr_drivingDistance()の結果をもとに凹包ポリゴンを生成するpgr_alphaShape()が用意されていますが、今回はPostGISのST_ConcaveHull()を使用して到達圏を作成しました。

-- すべてのエッジを囲む凹となるポリゴン(到達圏)を取得する
with concave_hull as (
select ST_ConcaveHull(ST_Collect(geom), 0.9,false) AS geom
from driving_distance_10min
group by from_v
)

-- ポリゴンの重なりが見づらいためディゾルブ
select st_union(geom)
into concave_hull
from concave_hull;

concave_hullを表示してみると、各駅を中心とした到達圏が作成されていました。

なお、地下鉄や市電の駅については問題なく到達圏が作成されていますが、JR白石駅などいくつかの駅については、線路の片側にしか到達圏が作成されていないことがわかります。

理由としては、それぞれの駅で一つのノードを始点として到達圏を作成しているため、線路によって道路ネットワークが分断されているようなケースだと、線路の片側にしか到達圏が作成されなかったりします。
より実態に即した到達圏を取得したい場合は、道路ネットワークを適宜修正したり、駅のすべて出入り口にポイントデータ作成した上で解析をすることで、より望ましい到達圏が取得可能かと思います。

以上、pgRoutingを使用した到達圏解析の紹介でした。

おまけ

OpenStreetMapのPOIデータを使用して、札幌市内のセイコーマートからの到達圏(徒歩5分)も取得してみました。

使用データ
http://download.geofabrik.de/asia/japan/hokkaido.html