はじめに

データ分析や機械学習では、csv ファイルのような表形式のデータで得られると何かと都合がよいです。リレーショナルデータベースのテーブルも表形式だから、そのまま使えそうです。

しかし、分析対象によっては、特殊なデータ形式を扱わなければならない場合があります。今回はその一例として、空間データを扱う PostgreSQL 拡張 PostGIS を紹介します。

空間データと、機械学習の例としてよく使われるアヤメのデータセットと比べてみましょう。

対して、図形データは以下のようになります。

PostGIS のテーブルには、図形を表す geometry カラムがありますが、それは上記のようにそのままではあまりよくわかりません。

実際の図形の情報が格納されていて、そこから長さなどを計算することはできますが、そのままでは学習用データとしても扱いにくいです。そこで、よく知られた数値に値を繰り出してくる必要があるわけです。

今回は、図形を表す geometry からよく知られた数字に繰り出す例を見ていきます。

環境づくり

GIS データをどのように扱うかを見てみる前に、実際に手を動かせる環境を作りましょう。ここでは、PostgreSQL の PostGIS という拡張を使います。

OpenStreetMapのPostGIS/Installationを参考に、インストールしておきます。

交差点の道路のなす角を求める

具体的に特徴量をひねり出す例として、ノードとエッジのなす角を求めてみたいと思います。

ここで、ノードは平面上の点で、x, y の座標で定まる点とします。エッジは、辺と辺を結ぶ線分とします。

具体例として、ノードは交差点、エッジは道路ととらえることで、地図として応用できるでしょう。

テーブル定義とデータ

辺 Edge と、辺が交わる点 Node からなる地図を考えます。Node は、2 点をつなぐ LINESTRING (折れ線は考えない) とします。

create extension postgis;
create table Node (
  id INTEGER PRIMARY KEY
  , geom GEOMETRY
);
create table Edge (
  id INTEGER PRIMARY KEY
  , node_from_id INTEGER
  , node_to_id INTEGER
  , geom GEOMETRY
  , FOREIGN KEY (node_from_id) REFERENCES Node(id)
  , FOREIGN KEY (node_to_id) REFERENCES Node(id)
);

INSERT INTO Node (id, geom) values
    (1, ST_Point(0, 0)),
    (2, ST_Point(1, 1)),
    (3, ST_Point(1, 2)),
    (4, ST_Point(3, 0));

INSERT INTO Edge (id, node_from_id, node_to_id, geom) values
    (1, 1, 2, ST_GeomFromText('LINESTRING(0 0, 1 1)')),
    (2, 2, 3, ST_GeomFromText('LINESTRING(1 1, 1 2)')),
    (3, 2, 4, ST_GeomFromText('LINESTRING(1 1, 3 0)'));

角度を求める

それでは、ひとつの節に3つの辺が接続するようなときの、辺同士のなす角を求めるクエリを実際に組み立ててみましょう。

ほしいもののイメージ

まずは、どのような属性の組が欲しいかイメージを持ちましょう。最終的に欲しいテーブルはこのようなイメージとします。

使用する関数

PostGIS で角度を求めるには、ST_Azimuthという関数があります。

第1引数を始点、第2引数を終点とするような線分と、y軸正方向の線分との時計回りのなす角を返す関数です。この関数を使って角度を求めていくことにします。

エッジの方向を求める

ST_Azimuth は2つの点で角度 (ラジアン) を求めるので、Edge の geom にある LINESTRING の始点と終点を取り出しましょう。

-- 片方向の辺を点に分解
with dumppoints as (
    select id, node_from_id, node_to_id, geom, st_dumppoints(geom) dump from edge
),
edge_points as (
    select
       id, node_from_id, node_to_id, geom, (dump).path[1] path, (dump).geom point
    from
        dumppoints
),

ST_DumpPoints を使うと、辺を構成する点が得られます。path には点の並び順が、geom には点の geometry の値が入っています。

次で、機械的に隣り合う点同士で、辺の方向を求めていきます。

-- 辺ごとの方向を求める
edge_azimuth as (
select
    ep1.id
    , ep1.node_from_id
    , ep1.node_to_id
    , ep1.geom
    , ep1.point point1
    , ep2.point point2
    , st_azimuth(ep1.point, ep2.point) / (2 * PI()) * 360 azimuth
from
    edge_points ep1
    inner join edge_points ep2
        on ep1.id = ep2.id and ep1.path = ep2.path - 1
),

エッジとエッジのなす角を求める

続けて、ノードと結合します。直前の azimuth の計算では、point1 (ep1.point) を始点として考えていたので、point1 を始点ととらえた場合は角度はそのままですが、point2 を始点ととらえた場合には、角度を180度ひっくり返します。

-- 逆向きの辺も考え、角度を補正する
bilateral_edges as(
select
    e.id, e.node_from_id, e.node_to_id, e.geom, e.point1, e.point2, e.azimuth, n.geom node_geom
from edge_azimuth e
inner join node n on e.node_from_id = n.id
where ST_Equals(point1, n.geom)

    union all
select
    e.id, e.node_to_id node_from_id, e.node_from_id node_to_id, e.geom, e.point2 point1, e.point1 point2
    , case when e.azimuth < 180 then e.azimuth + 180 else e.azimuth - 180 end azimuth
    , n.geom node_geomd
from edge_azimuth e
inner join node n on e.node_to_id = n.id
where ST_Equals(e.point2, n.geom)
),

いよいよ隣り合う角の角度を計算します.

-- ノードに流入する辺ごとに、角度で順序付けする
azimuth_order as (
select
    row_number() over (partition by node_from_id order by azimuth) azimuth_ord
    , count(*) over (partition by node_from_id) edge_count
    , be.id, be.node_from_id, be.node_to_id, azimuth
from bilateral_edges be
),
-- 隣り合う辺のなす角を求める
edge_angle as (
select
    ao1.node_from_id
    , ao1.azimuth_ord
    , ao1.id edge1_id, ao2.id edge2_id
    , case
        when ao2.azimuth_ord <> 1 then ao2.azimuth - ao1.azimuth
        else 360 + (ao2.azimuth - ao1.azimuth)
    end angle
from
    azimuth_order ao1
    inner join azimuth_order ao2
        on ao1.node_from_id = ao2.node_from_id
            and ((ao1.azimuth_ord + 1 = ao2.azimuth_ord) or (ao1.azimuth_ord = ao1.edge_count and ao2.azimuth_ord = 1))
where ao1.edge_count > 1
)

azimuth では、時計回りに角度が増加します。そこで、角度で順序を入れて、隣り合う角と JOIN して角度を計算します。

横持ち化

最後に、ほしい形に横持ち化します。

select
   e1.node_from_id node_id
   , e1.edge1_id edge1_id
   , e2.edge1_id edge2_id
   , e3.edge1_id edge3_id
   , e1.angle angle1
   , e2.angle angle2
   , e3.angle angle3
from
    edge_angle e1
    inner join edge_angle e2
        on e1.node_from_id = e2.node_from_id and e1.azimuth_ord=1 and e2.azimuth_ord=2
    inner join edge_angle e3
        on e1.node_from_id = e3.node_from_id and e1.azimuth_ord=1 and e3.azimuth_ord=3

得られる結果は次のようになります:

以上、今回は図形から特徴量を取り出す例を見てみました。

図形から特徴量を取り出す、ということを示す目的のため、単純化した想定でクエリが書かれており、場面によってはもう少し工夫が必要になるでしょう。

たとえば Edge が 2点からなる折れ曲がりのない LINESTRING である仮定や、ノードには高々3つのエッジのみが接続するという仮定のもとにクエリを書いていました。実際の地図などでは、4叉路、5叉路など、どういうものがあるのかあらかじめ調べてからクエリを組み立てていくことになるでしょう。

完成したクエリ

それでは、最後に完成したクエリを通してみてみましょう.

-- 片方向のエッジを点に分解
with dumppoints as (
    select id, node_from_id, node_to_id, geom, st_dumppoints(geom) dump from edge
),
edge_points as (
    select
       id, node_from_id, node_to_id, geom, (dump).path[1] path, (dump).geom point
    from
        dumppoints
),
-- 辺ごとの方向を求める
edge_azimuth as (
select
    ep1.id
    , ep1.node_from_id
    , ep1.node_to_id
    , ep1.geom
    , ep1.point point1
    , ep2.point point2
    , st_azimuth(ep1.point, ep2.point) / (2 * PI()) * 360 azimuth
from
    edge_points ep1
    inner join edge_points ep2
        on ep1.id = ep2.id and ep1.path = ep2.path - 1
),
-- 逆向きの辺も考え、角度を補正する
bilateral_edges as(
select
    e.id, e.node_from_id, e.node_to_id, e.geom, e.point1, e.point2, e.azimuth, n.geom node_geom
from edge_azimuth e
inner join node n on e.node_from_id = n.id
where ST_Equals(point1, n.geom)

    union all
select
    e.id, e.node_to_id node_from_id, e.node_from_id node_to_id, e.geom, e.point2 point1, e.point1 point2
    , case when e.azimuth < 180 then e.azimuth + 180 else e.azimuth - 180 end azimuth
    , n.geom node_geomd
from edge_azimuth e
inner join node n on e.node_to_id = n.id
where ST_Equals(e.point2, n.geom)
),
-- ノードに流入する辺ごとに、角度で順序付けする
azimuth_order as (
select
    row_number() over (partition by node_from_id order by azimuth) azimuth_ord
    , count(*) over (partition by node_from_id) edge_count
    , be.id, be.node_from_id, be.node_to_id, azimuth
from bilateral_edges be
),
-- 隣り合う辺のなす角を求める
edge_angle as (
select
    ao1.node_from_id
    , ao1.azimuth_ord
    , ao1.id edge1_id, ao2.id edge2_id
    , case
        when ao2.azimuth_ord <> 1 then ao2.azimuth - ao1.azimuth
        else 360 + (ao2.azimuth - ao1.azimuth)
    end angle
from
    azimuth_order ao1
    inner join azimuth_order ao2
        on ao1.node_from_id = ao2.node_from_id
            and ((ao1.azimuth_ord + 1 = ao2.azimuth_ord) or (ao1.azimuth_ord = ao1.edge_count and ao2.azimuth_ord = 1))
where ao1.edge_count > 1
)
-- 3つの流入を前提に
select
   e1.node_from_id
   , e1.edge1_id edge1_id
   , e2.edge1_id edge2_id
   , e3.edge1_id edge3_id
   , e1.angle angle1
   , e2.angle angle2
   , e3.angle angle3
from
    edge_angle e1
    inner join edge_angle e2
        on e1.node_from_id = e2.node_from_id and e1.azimuth_ord=1 and e2.azimuth_ord=2
    inner join edge_angle e3
        on e1.node_from_id = e3.node_from_id and e1.azimuth_ord=1 and e3.azimuth_ord=3

DATUM STUDIOでは様々なAI/機械学習のプロジェクトを行っております。
詳細につきましてはこちら

詳細/サービスについてのお問い合わせはこちら