真っ直ぐ道を探索するSQL

やりたいこと

「この道を真っ直ぐ」という言葉は誰でも一度は口にした言葉ですが、OpenStreetMapの道路データで、SQLで取得するにはどうしたらよいのか挑戦してみました。

解決すべき課題

探索処理の実現

SQLで次々と道路を結合する計算（探索）を実現してみようと思います。繰り返し処理が基本なので、今年のyoh_chan様のエントリー SpatiaLiteで再帰クエリを使ってポリゴンを半分に分割してみたでも登場した再帰クエリを使います。

交差点での道路選択

道路データは基本的に交差点で分割されています。
真っ直ぐに道路を接続するには、交差点に差し掛かった際に、接続する複数の道路から真っ直ぐ伸びる道路を選択する必要があります。

準備

データの準備

OpenStreetMapのShapeデータをダウンロードし、PostgreSQLにロードします。日本のデータはこちらのサイト:http://download.geofabrik.de/asia/japan.htmlからダウンロード可能です。

データベースの作成

１．PostGISのインストール
ShapeデータをPostgreSQLに投入する際には、空間情報処理を実現するPostGISを事前にインストールする必要があります.
CREATE EXTENSION postgis;

２．Shapeデータのインポート
QGISのSPITやDBマネージャーを利用すると、簡単にデータをPostgreSQLに取り込むことが出来ます。
　 参考 ： http://docs.qgis.org/2.6/ja/docs/training_manual/spatial_databases/import_export.html

３．空間インデックスの作成
　探索時に「道路が接合している」という処理を行います。ある程度の規模（数十万行以上）のデータベースでは、このような処理を行うには、空間インデックスを作成すると、高速化されます。
CREATE INDEX インデックス名 ON テーブル名 USING gist (空間データカラム名);

ストアド・プロシージャ（ユーザ定義関数の作成）

SQLが複雑・長文になるので、繰り返し登場する処理をユーザー定義関数として作成します。今回は、PL/pgSQLのような手続き言語ではなく、ピュアなSQLで記述してみました。（Sqliteでも互換性があればいいのですが。。。）

２線の接合角度を求める関数（ST_LineDiffDegree）

交差点で２つの直線が接合している際に、接合点を中心とした方位差の絶対値を度単位を返却します。一直線の場合は0度、直角に右折または左折している場合は90度という値が得られます。

CREATE_FRUNCTION_ST_LineDiffDegree.sql

CREATE OR REPLACE FUNCTION ST_LineDiffDegree( linestring_geom_4326_a  geometry ,  linestring_geom_4326_b  geometry  )
RETURNS double precision 
AS '
SELECT 
  ABS( CASE WHEN diff_degree_abs > 180 THEN diff_degree_abs -360
  ELSE diff_degree_abs END ) AS diff_degree_abs
FROM 
(
SELECT
    ABS(
        DEGREES( ST_Azimuth(   ST_PointN(base_geom , ST_NpointS(base_geom) -1 )  ,ST_EndPoint(base_geom) ))
        -
        ( CASE WHEN ST_StartPoint(geom) = ST_EndPoint(base_geom) THEN 
            DEGREES( ST_Azimuth(  ST_EndPoint(base_geom) , ST_PointN(geom ,2 ) ))
           ELSE 
            DEGREES( ST_Azimuth(  ST_EndPoint(base_geom) , ST_PointN(ST_Reverse( geom ) ,2 ) ))
           END 
         )
    ) AS diff_degree_abs
FROM 
 ( SELECT $1 AS base_geom , $2 AS geom ) AS base
) AS t1
;'
LANGUAGE SQL
IMMUTABLE
RETURNS NULL ON NULL INPUT;

ラインオブジェクトAの向きを保ったままラインオブジェクトBを結合する関数(ST_KDLineMerge)

　ラインオブジェクトは有向グラフです。また、複数のラインオブジェクトを接合するST_Linemerge関数では、接合するラインオブジェクトAとBのうちBの向きに影響を受ける特性がありました。
　真っ直ぐ探索するには、ラインオブジェクトAの向きを保たないと、途中で進行方向が変わってしまいます。

CREATE_FRUNCTION_ST_KDLineMerge.sql

CREATE OR REPLACE FUNCTION ST_KDLineMerge( linestring_geom_4326_a  geometry ,  linestring_geom_4326_b  geometry  )
RETURNS geometry
AS '
SELECT 
 ( CASE 
  WHEN 
   ST_EndPoint( $1 ) = ST_StartPoint( $2) THEN ST_linemerge( ST_Union($1 , $2 ) )
  WHEN 
   ST_EndPoint( $1 ) = ST_EndPoint( $2 )  THEN ST_linemerge( ST_Union($1 , ST_Reverse($2) ))
  ELSE 
   null
  END ) AS geo
;'
LANGUAGE SQL
IMMUTABLE
RETURNS NULL ON NULL INPUT;

真っ直ぐ探索SQL

準備編で用意したデータベースと、自作関数を駆使してSQLで探索処理を記述します。

クエリ設計

・再帰クエリの応用
　WITH RECURSIVE が再帰クエリの宣言です。何回繰り返すかは、WITH句の外のSELECT句のLIMITで指定します。

・GEOMETRY型の対応
　　再帰クエリのデータ型は全ての型が対応しているわけではなく、GEOMETRY型も再帰するデータセットincrement_n(n ,rank, a_id , base_geom , a_diff_degree) に記述する事は出来ません。今回はTEXT型にキャストすることで乗り越えました。

・配列型の活用
データを繰り返し接続していく過程を記録するために、配列を多用しています。

・出発地点の指定
道路リンクのID（gid）で指定します。
今回は荏田綱島線（神奈川県道102号)の終点リンクから出発してみます。

© OpenStreetMap contributors

SQL

GoStraightDown.sql

WITH RECURSIVE
increment_n(n ,rank, a_id , base_geom , a_diff_degree) AS (
    SELECT  
        1 AS n
        , 1::bigint AS rank
        , array[gid] AS a_id
        , the_geom::text AS base_geom
        , array[1:: double precision  ] AS a_diff_degree
        FROM osm_roads AS osm
        WHERE gid = 27075 -- 出発地点のリンクID
    UNION ALL
        SELECT 
         n+1
         , row_number()OVER( PARTITION BY n+1 ORDER BY ST_LineDiffDegree( base_geom , osm.the_geom )  ) AS rank
         , array_append( a_id , osm.gid ) AS a_id
         , ST_KDLineMerge( base_geom , ST_Reverse(osm.the_geom)  )::text  AS base_geom 
         , array_append( a_diff_degree , ST_LineDiffDegree( base_geom , osm.the_geom ) ) AS a_diff_degree
        FROM 
         (SELECT n , a_id , base_geom::geometry AS base_geom ,a_diff_degree FROM  increment_n WHERE rank= 1 ) AS t1 
         LEFT JOIN 
         osm_roads AS osm
        ON  
         ST_InterSects( ST_EndPoint( base_geom ) ,osm.the_geom ) 
         AND a_id[array_length(a_id , 1 )] != osm.gid
         AND ST_LineDiffDegree( base_geom , the_geom ) < 15 -- 直進の判断基準：２線の接合分前後の方位差が１５度以内
)
SELECT
 n , rank , a_id , ST_AsText( base_geom )
FROM
    (
        SELECT  * 
        FROM increment_n
        LIMIT 100 -- 道路を連結する上限数
    ) AS t1
WHERE base_geom IS NOT NULL
ORDER BY n DESC LIMIT 1 ;

実行結果

３カラムのデータが出力されます。

ｎ	a_id	ST_AsText
１４	27075,6464038...	LINESTRING(139.6415109 35.5432942...

n : 結合したリンクの数
a_id : 結合したリンクのgidの配列。
ST_AsText : 結合済みの道路形状のWKT

QGISでST_AsTextの文字列を表示してみます

© OpenStreetMap contributors

終点の分析

接合する角度を１５度に限定した結果、上下線分離リンクから上下線非分離のリンクに変更する地点で、閾値を超えてしまったようです。（真っ直ぐの定義を超えてしまったようです）

角度差の閾値設けないと、丁字路になるまで伸び続けます

© OpenStreetMap contributors

首都高３号線で日本縦断にチャレンジ

© OpenStreetMap contributors

御殿場で降りて,中央道にのって戻ってきてしいました。。。

皆様良いお年を