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