Mapillaryのmap_matchingライブラリを使ったPostGISによるマップマッチング

#はじめに
「マップマッチング」とは位置データから道路ネットワーク上の移動した経路を特定することを指します。 プローブデータと呼ばれる車両の位置データ記録を道路データとマップマッチングすることにより、道路の所要時間や走行速度といったパフォーマンス指標を算定します。

一般的な地図サイトでは、「経路探索」はサービスされています。経路探索は、pgRoutingやQGISのprocessingにある経路探索プラグインのように、あらかじめ指定した2点間の最短経路情報や経路候補を出力します。

一方マップマッチングによるプローブデータ解析では始終点の2点ではなく、移動体の走行軌跡とマッチする経路情報を得る必要があります。

目的にあったQGISのprocessingツールには、Offline-MapMatchingというのがありますが、2.18でなんとか動いていたものの3.xではフリーズしてしまいます。
また、範囲が大きくなるとかなりスローになってしまう特徴があります。
https://github.com/jagodki/Offline-MapMatching

ということで色々と調べた結果、今回mapillaryのmap_matchingライブラリを使用して、マップマッチングを行い、通過した道路リンク番号を得るプローブ解析を実装してみます。

#使用するもの

• Mapillary map_matching
  https://github.com/mapillary/map_matching

• python3
• psycopg2
• osm2pgrouting
• postgreSQL
• postGIS
• wget

#プログラム概要
STEP1: プローブの走行軌跡データCSVの取得
STEP2: マップマッチングするマップの取得と加工
STEP3: プローブ走行軌跡 CSVデータのpostgreSQLインポート
STEP4: map_matchingライブラリを使用してプローブデータをマップマッチング
STEP5: マップにマッチしたプローブデータとOSMの道路リンク番号を出力

という感じで進めてゆきます。

#事前作業
###postgreSQLのインストール
説明は省略します。
この辺を参考にしてください。
https://qiita.com/sanak/items/f50e8c5bb97bf9619958

windowsでもLinuxでもいいんですが・・・、とりあえずwindowsユーザーが多そうなので、
windowsベースと仮定して、次の通りにインストールするものとします。

インストール先：c:\pgsql\11.0
データベース名：mapmatching
ユーザー名：postgres
パスワード：merryxmas

###データベース作成
サンプルコードの実行用に下記のデータベースを作成します。

CREATE DATABASE mapmatching
    WITH 
    OWNER = postgres
    ENCODING = 'UTF8'
    LC_COLLATE = 'C'
    LC_CTYPE = 'C'
    TABLESPACE = pg_default
    CONNECTION LIMIT = -1;

#プローブの走行軌跡データCSVの取得
G空間情報センターから、加古川市_公用車走行データ（走行履歴）_2016_1をダウンロードします。
zipを解凍後のファイルから、probe_kaisen197_2016.csv.csvを使用します。
ファイル名が珍妙なので、probe_kaisen197_2016.csvにリネームしておきましょう。
中身を見ると、99,561レコードあるようです。

kaisen197,2016-01-04 15:25:40,34.748560,134.837897
...

のようなフィールドになっています。（確か2018年にダウンロードした時は年月日と時刻が分かれていたが、、、）

id,pdate,latitude,longitude
kaisen197,2016-01-04 15:25:40,34.748560,134.837897
...

のように1行目にヘッダをつけておきます。

加古川市のサンプル範囲の経緯度範囲は、(134.7641,34.6797,134.9534,34.861)くらいです。
プローブの走行軌跡データは日本測地系になっているので、あとでプログラム内部で世界測地系に変換します。

座標系の詳細は下記をご参照ください。
https://epsg.io/4301

#マップマッチングするマップの取得と加工
wgetコマンドを使用してPlanet OSM から Overpass API を使用してOSMデータをダウンロード必要な経緯度範囲のデータをダウンロードします。
加古川市のサンプル範囲の経緯度範囲は、(134.7641,34.6797,134.9534,34.861)なので、

wget -q --no-check-certificate --header=\"Content-Type: application/osm3s+xml\" https://overpass-api.de/api/map?bbox=134.7641,34.6797,134.9534,34.861 -O kakogawa.osm

で取得します。

##OSMデータのpostgresqlインポート
osm2pgroutingを使用して、OSMデータをpostgreSQLにインポートします。
実行時にpostgreSQLのbinディレクトリにインストールされたmapconfig.xmlファイルを指定します。

osm2pgrouting.exe -f kakogawa.osm -c <ファイルのあるディレクトリ>/mapconfig.xml --prefix kakogawa_ -d (データベース名) -U (ユーザー名) -W (パスワード) -h (ホスト名) --clean

正しく終了すると、kakogawa_ のついたテーブルがいくつかできます。
このうちkakogawa_waysを使用してマップマッチングを実行します。
osm2pgroutingについては下記を参照ください。
https://github.com/pgRouting/osm2pgrouting/wiki/Documentation-for-osm2pgrouting-v2.2

#プローブ走行軌跡 CSVデータのpostgreSQLインポート
ダウンロードしたプローブデータのCSVには連番がないので連番フィールドを追加します。
またプローブ解析に使用する経緯度ジオメトリを格納するためwayフィールドも追加します。

ダウンロードしたcsvファイルのあるフォルダで、psqlコマンドでpostgresqlに接続し、コマンドラインで操作を行います。

psql -U postgres mapmatching

###csvインポート先テーブル作成

CREATE TABLE public.probe_kaisen197_2016 (
	csv_id bigserial,
	id varchar NOT NULL,
	pdate timestamp NOT NULL,
	latitude double precision NOT NULL,
	longitude double precision NOT NULL,
	way geometry(POINT, 4326)  NULL,
	CONSTRAINT probe_kaisen197_2016_pkey PRIMARY KEY (csv_id)
);
CREATE INDEX probe_kaisen197_2016_way_idx ON probe_0122 USING gist (way);

###プローブデータコピー

\copy probe_kaisen197_2016(id,pdate,latitude,longitude) from 'probe_kaisen197_2016.csv' csv header

###ジオメトリ更新
プローブデータはEPSG:4301なので、wayで使用する座標系EPSG:4326に変換します。

UPDATE public.probe_kaisen197_2016 SET way=st_transform(ST_SetSRID(ST_Point(longitude,latitude),4301),4326) ;

##mapmatchコード作成
mapillary/map_maching/exampleのmap_matcher.pyをベースにコード作成します。

map_matcher.pyは、
経緯度だけが並んだCSVデータをstdinから読み込んで、配列に格納、
osm2pgroutingでインポートした道路データテーブルに対して、

• query_edges_in_sequence_bbox
  サーチ円の分拡張されたバウンディングボックス内のシーケンスのすべての道路エッジをクエリーする

• build_road_network
  エッジリストの双方向道路グラフデータを構築する

• query_candidates
  サーチ円内に存在するシーケンスデータの各々について計測データの候補をクエリーする
  クエリー結果をcandidateクラスに格納する

という処理を行っています。

コード変更するポイントは、

• 配列格納のデータをcsvからインポートしたテーブルに変更
• 配列からサブクエリーを作っていた処理をテーブルベースのクエリーに変更
• candidateクラスにマッチングした経緯度、ポイントの時刻、エッジデータかどうかを格納する

というところです。

###ソースコード変更
map_matcher.pyをmy_mapmacher.pyにコピーします。

準備したデータなどから使用するパラメータは下記の通りです。

データベース：mapmatching
テーブル：kakogawa_ways
プローブ：probe_kaisen197_2016

インポートライブラリを追加します。

my_mapmacher.py

import os
import sys
import datetime
import time
import locale
import psycopg2
import csv

メインの処理を下記のように変更します。
今回のプログラムでは引数の処理は行なっていないので、引数の評価部分は削除します。

my_mapmacher.py

def main(argv):

	# データベース接続パラメータ
	pguser='postgres'
	pgport='5432'
	pghost='localhost'
	pgdbname ='mapmatching'
	pgpassword='merryxmas'

	# postgresql://{username}:{password}@{hostname}:{port}/{database}
	dsn='postgresql://{0}:{1}@{2}:{3}/{4}'.format(pguser,pgpassword,pghost,pgport,pgdbname)

	# OSMデータダウンロード指定のファイル名から作成予定のOSMテーブル名を生成
	osmtbl='kakogawa_ways'

	# プローブCSVファイルをアップロードする
	csvtbl='probe_kaisen197_2016'

	#アプトプットファイル
	outputcsv='result_' + csvtbl + '.csv'

	# プローブテーブルを使用してマップマッチングを実行する
	# 下記のパラメータは必要に応じて変更しましょう
	# SEARCH_RADIUS :道路のサーチバッファ半径
	# MAX_ROUTE_DISTANCE ：サーチ候補の道路リンクの最大延長

	search_radius=30
	max_route_distance=2000

	start=time.time()
	conn = psycopg2.connect(dsn)
	candidates = map_match(conn, osmtbl,csvtbl, search_radius, max_route_distance)
	conn.close()

	process_time = time.time() - start

	print( 'process_time;',process_time )

	# 候補データに各エッジの最初と最後の識別フラグを追加する
	flg=0
	cb=None

	for candidate in candidates:
		candidate.edgeflg=(0 if flg == candidate.edge.id else 1)
		flg=candidate.edge.id
		if cb is not None :
			if candidate.edgeflg == 1 and cb.edgeflg==0 :
				cb.edgeflg=2
		cb=candidate

    print('process time min:{0}'.format(process_time/60))

	with open( outputcsv, "w" ) as f:
		f.write(u'mid,ptime,mlon,mlat,clon,clat,cid,cloc,cdist,edgeflg\n')
		for candidate in candidates:
			a= \
			'{0},'.format(candidate.measurement.id) +\
			'{0},'.format(candidate.ptime)+\
			'{0:.6f},{1:.6f},'.format(*map(float, (candidate.measurement.lon, candidate.measurement.lat))) +\
			'{0:.6f},{1:.6f},'.format(*map(float, (candidate.lon, candidate.lat)))+\
			'{0},'.format(candidate.edge.id) +\
			'{0:.2f},'.format(candidate.location) +\
			'{0:.2f},'.format(candidate.distance) +\
			'{0}\n'.format(candidate.edgeflg)
			f.write(a)
		f.close()

	return 0

パラメータsequenceはcsvテーブル名を渡し、map_match()の下位関数でテーブルを処理します。

マップマッチングの結果の返り値class Candidateにマッチングした経緯度、ポイントの時刻、エッジデータかどうかを示すフラグを追加します。

my_mapmacher.py

class Candidate(mm.Candidate):
    def __init__(self, measurement, edge, location, distance):
        super(Candidate, self).__init__(measurement=measurement, edge=edge, location=location, distance=distance)
        self.lon = None
        self.lat = None
        self.mlon=None	#マッチしたポイントの経度
        self.mlat=None	#マッチしたポイントの緯度
        self.ptime= None	#マッチデータの時刻
        self.edgeflg=None	#エッジかどうかを示すフラグ、osmデータのリンク番号を特定するために使用する

次に、関数query_edges_in_sequence_bbox() は、sequenceパラメータがプローブのデータテーブルになっていますので、配列からSQLを生成していた部分を下記のように書き換えます。

my_mapmacher.py

def query_edges_in_sequence_bbox(conn, road_table_name, sequence, search_radius):
    
    """
	サーチ円の分拡張されたバウンディングボックス内のシーケンスのすべての道路エッジをクエリーする
    Query all road edges within the bounding box of the sequence
    expanded by search_radius.
    """
    if not sequence:
        return

    stmt = '''
    -- NOTE the length unit is in km
    SELECT edge.gid, edge.source, edge.target, edge.length * 1000, edge.length * 1000
    FROM {road_table_name} AS edge
         CROSS JOIN (SELECT ST_Extent(ST_MakePoint(ST_X({sequence_name}.way), ST_Y({sequence_name}.way)))::geometry AS extent FROM {sequence_name}) AS extent
    WHERE edge.the_geom && ST_Envelope(ST_Buffer(extent.extent::geography, {search_radius})::geometry)
    '''.format(road_table_name=road_table_name,sequence_name=sequence,search_radius=search_radius)

    cur = conn.cursor()
    cur.execute(stmt)

    for gid, source, target, cost, reverse_cost in cur.fetchall():
        edge = Edge(id=gid,
                    start_node=source,
                    end_node=target,
                    cost=cost,
                    reverse_cost=reverse_cost)
        yield edge

    cur.close()

関数query_candidates()も同様にプローブのデータテーブルになっていますので、配列からSQLを生成していた部分を下記のように書き換えます。

my_mapmacher.py

def query_candidates(conn, road_table_name, sequence, search_radius):
    """
    サーチ円内に存在するシーケンスデータの各々の計測データの候補をクエリーする
    Query candidates of each measurement in a sequence within
    search_radius.
    """

    stmt = '''
        WITH 
	    --- WITH sequence AS (subquery here),
	    seq AS (SELECT *,
	                   ST_SetSRID(ST_MakePoint(ST_X({sequence_name}.way), ST_Y({sequence_name}.way)), 4326) AS geom,
	                   ST_SetSRID(ST_MakePoint(ST_X({sequence_name}.way), ST_Y({sequence_name}.way)), 4326)::geography AS geog
	        FROM {sequence_name})
	    
	    SELECT seq.csv_id, ST_X(seq.way) as lon, ST_Y(seq.way) as lat, seq.ptime,
	           --- Edge information
	           edge.gid, edge.source, edge.target,
	           edge.length, edge.length,
	
	           --- Location, a float between 0 and 1 representing the location of the closest point on the edge to the measurement.
	           ST_LineLocatePoint(edge.the_geom, seq.geom) AS location,
	
	           --- Distance in meters from the measurement to its candidate's location
	           ST_Distance(seq.geog, edge.the_geom::geography) AS distance,
	
	           --- Candidate's location (a position along the edge)
	           ST_X(ST_ClosestPoint(edge.the_geom, seq.geom)) AS clon,
	           ST_Y(ST_ClosestPoint(edge.the_geom, seq.geom)) AS clat
	
	    FROM seq CROSS JOIN {road_table_name} AS edge
	    WHERE edge.the_geom && ST_Envelope(ST_Buffer(seq.geog, {search_radius})::geometry)
	          AND ST_DWithin(seq.geog, edge.the_geom::geography, {search_radius})
    '''.format(road_table_name=road_table_name,sequence_name=sequence,search_radius=search_radius)

    cur = conn.cursor()
    cur.execute(stmt)

    for mid, mlon, mlat, mdt, \
        eid, source, target, cost, reverse_cost, \
        location, distance, \
        clon, clat in cur:

        measurement = Measurement(id=mid, lon=mlon, lat=mlat)

        edge = Edge(id=eid, start_node=source, end_node=target, cost=cost, reverse_cost=reverse_cost)

        assert 0 <= location <= 1
        candidate = Candidate(measurement=measurement, edge=edge, location=location, distance=distance)

        # Coordinate along the edge (not needed by MM but might be
        # useful info to users)
        candidate.lon = clon    # マッチングポイント X
        candidate.lat = clat    # マッチングポイント Y
        candidate.mlon = mlon   # プローブポイント X
        candidate.mlat = mlat   # プローブポイント Y
        candidate.ptime = mdt    # プローブ日付(TIMESTAMP)
        candidate.edgeflg = 0

        yield candidate

    cur.close()

#プログラムの実行と出力の確認
my_mapmacher.pyを実行します。

pytoh3 my_mapmacher.py

私のPCの仮想環境ではprobe_kaisen197_2016のデータで8分くらいでマッチングが終了します。

result_probe_kaisen197_2016.csv

mid,ptime,mlon,mlat,clon,clat,cid,cloc,cdist,edgeflg
1,2016-01-22 09:12:53,134.830824,34.750382,134.830894,34.750255,24122,0.35,15.33,1
2,2016-01-22 09:12:54,134.830584,34.750102,134.830591,34.750090,24122,0.77,1.51,0
3,2016-01-22 09:12:55,134.830495,34.750122,134.830530,34.750056,24122,0.85,7.90,2
4,2016-01-22 09:12:56,134.830285,34.749842,134.830253,34.749899,24121,0.96,6.84,1
5,2016-01-22 09:12:57,134.830174,34.749802,134.830151,34.749842,24121,0.93,4.78,0
6,2016-01-22 09:12:58,134.830054,34.749732,134.830030,34.749774,24121,0.91,5.02,0
.....

のようなデータが出力されます。
csvファイルのうち、2016-1-22のデータを一部抽出しています。

edgeflgは、1=開始ポイント、0=リンク内ポイント、2=終了ポイントになっています。
解析結果について、シーケンスと時刻が概ね連続しているcidが同じレコードブロック内で、edgeflg=2のptimeからedgeflg=1のptimeを引いた値がリンク所要時間になります。
出力サンプルでは、24122のリンク所要時間は2秒です。

#解析結果のチェック
出力されたマップマッチング結果をQGISで評価してみましょう。
result_probe_kaisen197_2016.csv の結果から2016-1-22のデータを抽出して
2016-1-22 9:12:53 から 2016-1-22 16:45:31の2266データをQGIS上にポイント表示します。

全体として、概ねうまくマッチングができています。
細かく確認しましょう。

##リンクがうまくマッチしている例
09:14:31 に道路リンクに開始マッチし、09:14:38 に終了マッチしています。
09:14:30、09:14:39がほぼリンク上であったはずですが、位置が別リンクに近かったため、別リンク番号がマッチしています。

##プローブデータが荒れたため別のリンクにマッチしている例
09:15:03から09:15:26のデータはGPSデータが荒れたため、一部のデータが併走する高架上の道路にマッチされています。

##OSMで道路リンクがないため、データがマッチングされていない例
16:31のデータの一部は、OSMデータ上に道路リンクがないためにマッチングされていません。

#まとめ
ご紹介したMapillaryのmap_matchingライブラリを利用したpython3とpostGISのマップマッチングは、カスタマイズがしやすく、実行スピードもかなり速いです。

プローブデータはGPSのHzに応じて一定時間に出力されるポイント数が増えます。
1Hzで1日分だと3600*24=86400レコード、もしも100Hzだとこの100倍のポイント数を処理することになり、マッチング対象の道路リンクデータのレコード数も考慮すると処理にかなりの時間がかかります。
試してはいませんが、100Hzのデータで同じ範囲のマップマッチングを行った場合、800分くらい必要？
2017年の加古川のプローブデータが95ファイルで1.45GBですから100Hzのデータと同じ程度のレコード数としてこれを全て登録してマップマッチングするとどのくらいかかるかわかるでしょう。

#補足：map_matchingライブラリの修正

Mapillaryのmap_matchingはpython2系用なので、python3用にすこし修正が必要です。

map_matching.py

@7,8,9
from . import shortest_path
from . import viterbi_path
from . import road_routing

road_routing.py

@4
from . import shortest_path as sp

viterbi_path.py

@80
 def __next__(self):
@214
try:
    next_state = next(states)
except StopIteration:
    return

この記事のソースコードは下記にあります。
https://github.com/picaosgeo/mymapmatching