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

:::note info
これは [MIERUNE AdventCalendar 2022](https://qiita.com/advent-calendar/2022/mierune) 15日目の記事です。
昨日は @dayjournal さんによる [QGISのプロセッシングを色々とためしてみた](https://qiita.com/dayjournal/items/dc2009cf6c30aa7f978c) でした。
:::

# 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万件
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/168196/a0e85525-81e0-28e1-8e59-d4ec12ec623a.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/168196/e25fed5a-e2e4-79ca-bd0e-2c9c38222c0a.png)


- 日本の市区町村Polygon
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/168196/b8957a84-6cae-f7d2-1d0f-1861424f1f79.png)



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

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

↓ 結果イメージ
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/168196/2e813060-62e4-547c-9079-d6a3f212ff93.png)



## 3.1. 緯度経度カラムのままで
```sql: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時間でいけました... :confounded:
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/168196/09234e72-a47a-2739-d42c-f246ec0b77df.png)


## 3.2. Pointのgeoｍetryカラムで

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

```sql: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 !
```sql: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分… :confounded:
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/168196/4ccc4e20-1010-b4aa-e385-494c13226f55.png)



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

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

```sql:add_index
CREATE INDEX gps_spatial_index ON scratch.gps_strava USING GIST (geom);
```

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

→結果 : 3分59秒！
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/168196/6eb68816-b820-40fd-1dbb-65cc92e4a0c8.png)

## 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
```sql: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秒
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/168196/cd318b21-0b0c-7549-4363-6a4617940b2d.png)

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

```sql: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)))')
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/168196/7e87bd6d-7d6a-56b1-7381-7bc462f69013.png)


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

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


```python: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秒で完了できました！:metal_tone2:
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/168196/bdf10e82-00de-0ce4-4d57-452efb76ce4f.png)

#### 3.4.4. GIST Indexの影響
- 処理の行動を見ると、対象市区町村にGPSデータが無い場合すぐにスキップされて、データがたくさん入ってると難病かかります
- GIST Indexを抜けると対象市区町村にGPSデータが無くても処理時間が約2秒かかりまして、無駄になります↓
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/168196/09e0395e-405f-a4a8-4629-a36c514038de.png)


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

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


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


## おまけ ： ST_INTERSECTSでどうになるのか
ST_WITHIN ではなくST_INTERSECTSで実行したらスピードがあんまり変わらない

```sql: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秒

```sql: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秒

以上 :rocket:

:::note info
明日は@@xinmiao1995さんによる[PLATEAUの3D都市モデルでMIERUNEマステープの作り方](https://qiita.com/xinmiao1995/private/6cb6aad8b397526c0b17)です！お楽しみにー
:::