これは MIERUNE AdventCalendar 2022 15日目の記事です。
昨日は @dayjournal さんによる QGISのプロセッシングを色々とためしてみたでした。

1. 概要

GPS log dataが基本的にlongitude,latitude,timeだけの情報がありますけど、市区町村に情報がよく足りない
Reverse geocoding (逆ジオコーディング)処理で各GPS logに市区町村を判定できる
ただ、GPS dataを集まるとVolumeが大きくなって処理が重たくなる

→ この記事はPostGISでどうやって逆ジオコーディングのパフォーマンスをOptimizeするのかを説明

2. Sample Data

PostGISに格納しているの：

Strava Running AppliのSessionを数件まとめてGPX fileを取得し、テスト用に重複して、約1000万件

日本の市区町村Polygon

3. PostGISで逆ジオコーディングQueries

基本的にPostGIS機能
ST_WITHIN({point geometry}, {polygon geometry})
を使って、PointとGeometryの書き方でパーフォーマンスを確認しましょう

↓ 結果イメージ

3.1. 緯度経度カラムのままで

Reverse_geocoding_method_1

UPDATE scratch.gps_strava g
SET city_code = c.city_code,
pref = c.pref_code,
city = c.city_name
FROM geo.japan_cities c
WHERE ST_WITHIN(ST_POINT(g.longitude,g.latitude)::geography::geometry,c.geom)

→結果 : 3時間でいけました...

3.2. Pointのgeoｍetryカラムで

先にGPS dataにGeometryカラムを追加して計算する

Reverse_geocoding_method_2

-- add geom
ALTER TABLE scratch.gps_strava ADD COLUMN "geom" geometry;

-- calculate geometry
UPDATE scratch.gps_strava 
SET geom = ST_POINT(longitude,latitude)::geography::geometry;

Let's try !

Reverse_geocoding_method_2

UPDATE scratch.gps_strava g
SET city_code = c.city_code,
pref = c.pref_code,
city = c.city_name
FROM geo.japan_cities c
WHERE ST_WITHIN(g.geom,c.geom)

→結果 : 4倍速い！けど、まだほぼ45分…

3.3. GPS dataのGeometryカラムに事前にGIST INDEXを付ける

GPSデータのテーブルにgeometryカラムにSpatial indexを追加したら？
参考 : https://www.postgresql.jp/document/8.3/html/textsearch-indexes.html

add_index

CREATE INDEX gps_spatial_index ON scratch.gps_strava USING GIST (geom);

そしてまた UPDATEやってみよう

→結果 : 3分59秒！

3.4. GPSデータにGIST Index付けたGeometryカラムと市区町村のEWKTテキスト形式で記入

1000万件のGPS dataで約4分でいけましたけど、データが1億件になったらまた時間かかりそうで、以下の改善方法をやってみる

1市区町村判定だけのQueryになりますけど、栃木県益子町でやってみよう
手動でDatabaseに益子のポリゴンをEWKT形式でコピーして、以下みたいにPasteする
ST_ASEWKT(geom)でgeometryをEWKT形式で取得できます。
比較しましょう

3.4.1. 1市区町村で比較テスト市区町村のGeometry column

Reverse_geocoding_method_341_mashiko_geom

UPDATE scratch.gps_strava g
SET city_code = c.city_code,
pref = c.pref_code,
city = c.city_name
FROM geo.mashiko c
WHERE ST_WITHIN(g.geom,c.geom);

→ 13秒

3.4.2. 1市区町村で比較テスト市区町村のGeometry column

Reverse_geocoding_method_342_mashiko

UPDATE scratch.gps_strava
SET city_code = '09342',
pref = '栃木県'
city = '益子町'
WHERE ST_WITHIN(geom,'SRID=4326;MULTIPOLYGON(((140.140443696 36.525078721,140.138973696 36.523677892,140.138593113 36.522226505,{座標が多すぎて一部だけ} 36.525828721,140.139953139 36.525421225,140.140443696 36.525078721)))')

→ 結果 : その市区町村で 7秒でいけました

3.4.3. 全市区町村にすると

→ 全市区町村にしたらQueryを1909回実行しないといけない...それはScriptで実行してみます。

reverse_geocoding_algorithm

import psycopg2
import psycopg2.extras
import time

start_time = time.time()

conn = psycopg2.connect(
        host="localhost",
        database="postgres",
        port="5432",
        user="bordoray",
        password="bordoray's password",
    )

cur = conn.cursor(cursor_factory=psycopg2.extras.DictCursor)

# get data from database
query = "SELECT city_code,pref_name,city_name,ST_ASEWKT(geom) as geometry FROM geo.japan_cities WHERE city_name IS NOT NULL ORDER BY city_code"
cur.execute(query)
res = cur.fetchall()

# update each cities
for city in res:
    sql = "UPDATE scratch.gps_strava "
    sql += "SET city_code = '"+city["city_code"]+"', "
    sql += "pref = '"+city["pref_name"]+"', "
    sql += "city = '"+city["city_name"]+"' "
    sql += "WHERE ST_Within(geom,'"+city["geometry"]+"') "

    cur.execute(sql)
    print(city["city_code"], city["pref_name"], city["city_name"])
print("--- %s seconds ---" % (time.time() - start_time))

→ 結果 : GPSにGIST indexが付けたら上記みたいなAlgorithmで約3分16秒で完了できました！

3.4.4. GIST Indexの影響

処理の行動を見ると、対象市区町村にGPSデータが無い場合すぐにスキップされて、データがたくさん入ってると難病かかります
GIST Indexを抜けると対象市区町村にGPSデータが無くても処理時間が約2秒かかりまして、無駄になります↓

4. まとめ：PostGISで逆ジオコーディングの UPDATE Queries評価まとめ

GPSデータ量が大きい時
-　GPS data にGIST INDEXを付けないと
-　各市区町村polygonでEWKT形式でテキストで記入したらもっと速くなりそう

Method	Indexなし	Index付き
① GPSデータの緯度経度カラムベース `ST_WITHIN(ST_POINT(gps.longitude,gps.latitude) ::geography::geometry,city.geom)`	3時間7分	-
② GPSデータのGeomカラムベースと市区町村のGeomカラム `ST_WITHIN(gps.geom,city.geom)`	44分48秒	3分59秒
③GPSデータのGeomカラムベースと市区町村のEWKT Geometry `ST_WITHIN(gps.geom,'SRID=4326...')`	1時間10分	3分16秒

おまけ： ST_INTERSECTSでどうになるのか

ST_WITHIN ではなくST_INTERSECTSで実行したらスピードがあんまり変わらない

ST_WITHIN

UPDATE scratch.gps_strava g
SET city_code = c.city_code,
pref = c.pref_code,
city = c.city_name
FROM geo.japan_cities c
WHERE ST_WITHIN(g.geom,c.geom)

Result : 239秒

ST_INTERSECTS

UPDATE scratch.gps_strava g
SET city_code = c.city_code,
pref = c.pref_code,
city = c.city_name
FROM geo.japan_cities c
WHERE ST_INTERSECTS(g.geom,c.geom)

Result : 251秒

以上

明日は@@xinmiao1995さんによるPLATEAUの3D都市モデルでMIERUNEマステープの作り方です！お楽しみにー

PostGISでGPS dataに逆ジオコーディング方法評価