Postgisを使って緯度経度から２点間の距離を計算した

#背景
インターン先のサービス提携先の各店舗の商圏内にいるユーザーIDとその店舗を算出して個別のユーザーアプローチを取ることになったので、やってみました。
最初はmetabaseからユーザーと店舗の位置情報テーブルを個別に二つcsvファイルで用意してpythonで読み込んで計算したんですが、postgresqlを採用してたのでpostgisを使ってmetabase内で計算することにしました。
metabaseは社内で利用している無料のBIツールです。
pythonで計算した時の参考記事はこちら

#やりたいこと
各ユーザーと1km以内にある店舗の組み合わせを全て出したい

#データ
・faclility_id,lat,longの３列＊店舗数のfacilitiesテーブル
・user_id,lat,longの３列＊ユーザー数のuserテーブル

#コード

distance

SELECT
    dist.user_id
    dist.facility_id
    dist.distance
FROM
    (SELECT
        users.user_id AS user_id,
        facilities.facility_id AS facility_id,
        st_distance_sphere(st_makepoint(users.long, users.lat),
                           st_makepoint(facilities.long, facilities.lat)) AS distance,
    FROM users
    CROSS JOIN facilities) AS dist
WHERE dist.distance < 1000;

CROSSJOINでユーザーレコードと店舗レコードの全組み合わせのレコードを作り、その全部に対してpostgisのst_distance_sphere関数を使うことで全組み合わせの距離を算出しています。

#やってみて
この関数では完全な球形に対して計算しているので回転楕円体である地球により近似するためにはst_distance_spheroid関数を使うとより精密になりますが、今回の計算で算出する距離は長くとも100km程度で、欲しいのは1km以内かどうかなのでこれで十分かと思います。

ちなみに、pythonでやってみた時との誤差は数m程度だったので、この誤差が上記の厳密さに起因しているのかと思いました。

参考にした記事
postgresql-PostGisの距離計算