City Street Orientations 街の道路の向きの調べ方

9142 ワード
GIS OpenStreetMap PostGIS OpenStreetMap テキストリンク

道路の向きを集計する

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

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

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

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

データ入手

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

.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でもそれっぽいことができました。