はじめに

緯度・経度があるデータに、その位置に関する情報を付与したいことがあります。
国土数値情報　ダウンロードサービスには、色んな情報が無償で公開されています。それを使って、最寄り駅とその距離を求めます。

PostGISのインストール

PostGISを使います。 PGSQLのサイトからだと便利です。
（以前はPowerShell用の手順もあったのですが、WSLを使いましょうになってきました。。。）

シェイプファイルからインポート

駅別乗降客数データをダウンロードします。展開すると、*.shpなどのファイルがあります。ここに移動して、シェイプファイルからインポート用のSQLに変換します。

$ shp2pgsql -W cp932 -D -I -s 4326 S12-18_NumberOfPassengers.shp station_passenger > station_passenger.sql

$ psql -U postgres sample < station_passenger.sql

距離をもとめる関数を使うために、geography型のカラムを追加します。

psql> ALTER TABLE station_passenger ADD COLUMN geo geography(MultiLineString,4326);
psql> UPDATE station_passenger SET geo = geom::geography;
psql> CREATE INDEX idx_station_passenger_geo ON station_passenger USING gist (geo);

最寄り駅とその距離を算出

緯度・経度があるデータを、下記のようなテーブルに作ってpostgisにいれます。

create table sample (
    id integer primary key,
    longitude double precision,
    latitude double precision,
    geo geography(Point, 4324)
)

緯度・経度から、geography型のカラムに値をいれ、インデックスをつけます。

psql> create index sample_geo on sample using gist(geo);
psql> update sample set geo = st_setsrid(st_point(longitude, latitude), 4326)::geography;

cross join と、遅延評価のlateralを使います。
イメージ的には、sampleテーブルの1行ごとに、最短距離のレコードを問い合わせます。
geographyのカラムにインデックスが貼っていないと、めちゃめちゃ遅くなります。

sql

select a.id, station_name, passengers_per_day, distance from sample a 
CROSS JOIN LATERAL (
SELECT 
S12_001 as station_name, 
S12_033 as passengers_per_day, 
a.geo <-> b.geo as distance 
FROM station_passenger b ORDER BY a.geo <-> b.geo LIMIT 1) AS closest_pt

copy文で出力するには、下記のようになります。

copy (select a.id, S12_001, S12_033, dist from rpay_geo a CROSS JOIN LATERAL (SELECT S12_001, S12_033, a.geo <-> b.geo as dist FROM station_passenger b ORDER BY a.geo <-> b.geo LIMIT 1) AS closest_pt) TO '/tmp/nearest_station_passenger.csv'(FORMAT csv, FORCE_QUOTE *);

参考リンク

https://qiita.com/boiledorange73/items/f2d8388e88d534066bc0
https://stackoverflow.com/questions/48551479/postgis-minimum-distance-between-two-large-sets-of-points/48554474
https://qiita.com/junkor-1011/items/896f2071546b5dc1c3e5