道路の向きを集計する

元ネタは、下記のBlogです。

Comparing City Street Orientations - Geoff Boeing
http://geoffboeing.com/2018/07/comparing-city-street-orientations/

　GIGAZINEでも掲載されましたが、街の道路の向きを調べると、都市の複雑さを観察することができるというものです。

　記事のなかでは、OpenStreetMapのデータと、Pythonを使っているということですが、PostGISで挑戦してみます。

データ入手

OpenStreetMapのダウンロードサイトがありますので、調査したいエリアを選択して、入手します。
日本は、国全体ではなく地方毎に分割しています。

日本のデータのダウンロード先
http://download.geofabrik.de/asia/japan.html

.shp形式をダウンロードすれば、フリーのGISソフトのQGISで表示可能です。

QGISでの表示例

 © OpenStreetMap contributors

データのロード

PostGISで処理するにために、データベースに読み込む必要があります。

範囲選択

QGISの範囲選択ツールで調査対象エリアを囲みました。

DBへの投入

DBツールのインポート機能を用いて、選択範囲のデータをPostgreSQLに投入します。

集計方針

道路の方向の計算方法

　道路は交差点から交差点までを線で結びますが直線とは限りません。また、OpenStreetMapの道路形状は、必ず交差点で分割されているわけではありません。

　１）道路形状を構成する点で分割し
　２）２点間で線分で結び、方位と線分長の算出
　３）方位毎に線分長を集計する
　　　　　（一方通行は１倍、両方向通行道路は１８０度反転した方位でもカウント）

集計クエリ

count_street_orientations.sql


WITH step2 AS 
(
	-- ２点を結び、線分長と方位を計算
	SELECT
	   id 
	 , path
	 , degrees(ST_Azimuth( geom , LEAD( geom)OVER( PARTITION BY id ORDER BY path ) ) ) AS dir
	 , ST_Distance( geom::geography , LEAD( geom::geography )OVER( PARTITION BY id ORDER BY path ) ) AS dist_meter
	 , oneway
	FROM
	-- 集合体を個々の点にバラす
	(
		SELECT 
		 id , ((pd).path)[1] , (pd).geom AS geom ,oneway
		FROM 
		-- ST_DumpPointsで道路中心線を点に分割した集合体に変換
		(
			SELECT 
				id 
				, ST_DumpPoints(geom) AS pd 
				, oneway
			FROM public.tuduki_road
		) AS base
	) AS step1
) 
-- 順行・逆行の値を集計
SELECT
 degree_10 , SUM( sum )
FROM 
(
	-- 一方通行・両方通行の順方向の集計
	SELECT 
	 TRUNC( dir :: numeric , -1 ) AS degree_10 , SUM( dist_meter )
	FROM step2
	WHERE dist_meter IS NOT NULL 
	GROUP BY 1

	UNION ALL 

	-- 両方向通行の道路方向の反転し集計
	SELECT 
	 TRUNC( CASE WHEN dir :: numeric  < 180 THEN dir :: numeric  + 180  ELSE dir :: numeric  - 180 END , -1 ) AS degree_10 
	 , SUM( dist_meter ) 
	FROM step2
	WHERE dist_meter IS NOT NULL 
	 AND oneway = 'B'
	GROUP BY 1
) AS step3
GROUP BY 1
ORDER BY 1
  ;

結果

鶏頭図が手元のアプリですぐに描画できなかったので、EXCELのレーダーチャートとして出力しました

まとめ

簡単ではありますが、方位と道路長の集計を実施しました。
SQLでもそれっぽいことができました。