ユーザ関数で2点間距離算出

SQLにおいて与えられた2つの緯度経度より2点間距離を求める

サーバを整理していたら、1年以上前に開発したものが見つかったのでメモ。

指定された2つの緯度経度によりをSQLにて距離を求めたいという要望で
C言語にてpostgreSQLのユーザ関数を作成。
やりたい事は以下の通り。
以下のように緯度、経度を保持するテーブルがあって

                            Table "test"
   Column    |           Type           |   Modifiers   | Description
-------------+--------------------------+---------------+-------------
 test_id     | integer                  | not null      | テストID
 ido         | character varying(32)    |               | 緯度
 keido       | character varying(32)    |               | 経度

ある緯度・経度の地点より3km以内のtest_idを取得したい!
とか
ある緯度・経度の地点よりどのくらいの距離にあるの?
下記SQLのrange関数を自分で作っちゃえというもの。

select
 test_id 
from test 
where range('E141.16.55.6144', 'N43.6.21.8016', keido, ido) <= 3;

select 
 range('E141.16.55.6144', 'N43.6.21.8016', keido, ido) AS kyori
from
 test;

まずはC言語で、2つの緯度・経度の距離を算出するモジュールを作成。
定義だけ見て、自分で作成しつつちょっとびびる、かなりの数学。

/**
 * ヒュベニの公式を用いて経緯度から2点間距離を算出
 * @param   dLon1   第1地点経度
 * @param   dLat1   第1地点緯度
 * @param   dLon2   第2地点経度
 * @param   dLat2   第2地点緯度
 * @return  2点間距離[m]
 *
 * D = sqrt((M * dP) * (M * dP) + (N * cos(P) * dR) * (N * cos(P) * dR))
 * D  : 2点間の距離(m)
 * P  : 2点の平均緯度
 * dP : 2点の緯度差
 * dR : 2点の経度差
 * M  : 子午線曲率半径
 * N  : 卯酉線曲率半径
 *
 * M = 6334834 / sqrt((1 - 0.006674 * sin(P) * sin(P)) ^3)
 * N = 6377397 / sqrt(1 - 0.006674 * sin(P) * sin(P))
 *
 * Hubeny式については, 以下参照
 * http://www.kashmir3d.com/kash/manual/std_siki.htm
 */

でモジュール作って、今度はPostgreSQL用の関数を作成。
下記の場合は、SQLでコールされると常に100を返すサンプル

#include <postgres.h>
#include <fmgr.h>
/* プロトコルバージョン指定 */
PG_FUNCTION_INFO_V1(pg_range);

Datum
pg_range(PG_FUNCTION_ARGS){
    PG_RETURN_INT32(100);
}

コンパイルしてライブラリを作成

gcc -shared -Wl,-soname,pg_range.so -o pg_range.o range.o -lm -I${PG_HOME}/include pg_range.c

PostgreSQLに関数を定義

CREATE FUNCTION range(text, text, varchar, varchar)
RETURNS INTEGER
AS '/tmp/pg_range.so', 'pg_range'
LANGUAGE C
WITH (iscachable);

これで、SQLで独自に作成した range 関数を使用できます。
これにより、ある郵便番号より緯度・経度を取得して
その緯度・経度から10km圏内の店舗等の情報をDBから取得可能になりました。

複雑なSQLが必要な場合などには、独自でC言語で関数を作れば意外と便利。