この記事はYAMAPアドベントカレンダー2023に参加しています。
こんにちは、今年4月入社のKomatsuです。YAMAPではバックエンドを担当して開発しています。
はじめに
与えられた座標を使って、データベースにある座標の集合から最短距離の座標を検索したいと思います。
Nearest Neighbor Search(最近傍探索)と呼ばれるものを使います。
PostGISでは<->演算子とORDER BY句をあわせて使うことで、距離順の座標をGiSTインデックスを用いて高速に検索することが可能です。
参考: 29. Nearest-Neighbour Searching — Introduction to PostGIS
準備
今回は、国土交通省の国土計画局さんが公開されている郵便局データを使います。
Docker環境を使ったレポジトリを作りました。
セットアップ手順はREADMEを参照してください。
起動スクリプトについて補足します。こちらの記事を参考に作成しました。
#!/bin/sh
# Make working directory
mkdir /tmp/work && cd /tmp/work
# Create postgis extension
psql -U $POSTGRES_USER -d $POSTGRES_DB -c "CREATE EXTENSION postgis"
# Get data from website
wget https://nlftp.mlit.go.jp/ksj/gml/data/P30/P30-13/P30-13_40.zip
unzip P30-13_40.zip
# Convert to SQL
shp2pgsql -W cp932 -D -I -s 4612:4326 ./P30-13_40/P30-13_40.shp > ./init.sql
# Execute SQL to import data
psql -U $POSTGRES_USER -d $POSTGRES_DB -f ./init.sql
# Remove working directory
cd / && rm -rf /tmp/work
PostGISをインストールするとshp2pgsqlという、シェイプファイルをダンプするためのSQLに変換するコマンドがついてきます。
shp2pgsql -W cp932 -D -I -s 4612:4326 ./P30-13_40/P30-13_40.shp > ./init.sql
オプションはこちらの記事を参考にさせていただきましたが、一点、空間参照系を4612から4326に変換しています。pgAdminのGeometry Viewerで背景地図を表示させるためです。
実践
まずインポートしてきた郵便局のデータを見ます。
今回は福岡県内にある郵便局のデータを入れました。
geomカラムには、ジオメトリ型で郵便局の座標(緯度・経度)の情報が格納されています。
今回起点とする座標として、大濠公園の真ん中の島の座標(33.586057, 130.376449)を使います。
次にこの座標を使って、最も距離の近い郵便局の座標を算出します。
SELECT
p30_005 as name,
geom <-> ST_SetSRID(ST_Point(130.376449, 33.586057), 4326) as dist,
geom
FROM
"p30-13_40"
ORDER BY
dist
LIMIT
10
;
今回はわかりやすくLIMIT 10を指定しています。
この<->はST_Distanceと同じようにジオメトリ間の平面距離の最小値を返しますが、ST_Distanceと比べて高速です。ORDER BYを使うことで、Gistインデックスを用いた距離順のデータが返されます。
結果がこちら。
大濠公園の真ん中の島からは、大濠郵便局が一番近いということがわかりました。
またEXPLAIN ANALYZE句でp30-13_40_geom_idxというGiSTインデックスが使われていることが確認できます。
Limit (cost=0.14..7.45 rows=10 width=61) (actual time=0.062..0.077 rows=10 loops=1)
-> Index Scan using "p30-13_40_geom_idx" on "p30-13_40" (cost=0.14..595.14 rows=814 width=61) (actual time=0.061..0.075 rows=10 loops=1)
Order By: (geom <-> '0101000020E610000095F3C5DE0B4C6040C24B70EA03CB4040'::geometry)
Planning Time: 0.111 ms
Execution Time: 0.100 ms
mydb=# \d p30-13_40_geom_idx
Index "public.p30-13_40_geom_idx"
Column | Type | Key? | Definition
--------+--------+------+------------
geom | box2df | yes | geom
gist, for table "public.p30-13_40"
参考: <-> | Postgis.net
参考: ST_Distance
まとめ
- 算出するのは純粋な距離の差なので、道路も考慮できたら面白いかなあ。🤔
- GISの勉強を精進していきたいと思いました。🧐
出典
- 出典：国土交通省国土数値情報ダウンロードサイト(https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-P30.html)
Discussion