【postgreSQL】point型からの地域標準メッシュコードの取得

PostgreSQL Advent Calender 用のエントリです。
エントリー当初は、slony-I on Iaasの運用を書こうと思ってましたが、急遽、別のネタが思い浮かんだので、位置屋らしく本件でエントリします。

地域標準メッシュコードとは?

地域メッシュ - Wikipedia
日本地図の上に碁盤の目を張り巡らせて、1つ1つの枠にユニークなコードを割り振ました。ということですね。
以前、pythonで取得した記事となります。
【python】緯度経度から標準地域コードの取得

今回はpostgreSQLの上で取得してみようと思います。

準備

テーブルの作成

create table landmark ( name text, point point not null, mesh text);

データも登録しましょう

insert into landmark( name, point )values( '渋谷区役所', point(35.6640352,139.6982122) );
insert into landmark( name, point )values( '札幌市役所', point(43.0620958,141.3543763) );
insert into landmark( name, point )values( '那覇市役所', point(26.2285302,127.6891103) );

select * from landmark;
    name    |          point           | mesh
------------+--------------------------+------
 渋谷区役所 | (35.6640352,139.6982122) |
 札幌市役所 | (43.0620958,141.3543763) |
 那覇市役所 | (26.2285302,127.6891103) |

1次メッシュコードの取得

1次メッシュコードは、緯度 * 1.5の値を切り下げた値 || 経度 - 100の整数値となるので、

select cast( floor( point[0] * 1.5 ) as text)||cast( floor( point[1] - 100 ) as text) from landmark;
 ?column?
----------
 5339
 6441
 3927

簡単ですね。
せっかくなので、関数登録しておきましょう。

CREATE FUNCTION to_1st_mesh(point) RETURNS text AS '
declare
    geo alias for $1;
    left  text;
    right text;
begin
    left  = cast( floor(geo[0] * 1.5) as text );
    right = cast( floor(geo[1] - 100) as text ); 
    return left || right;
end;
' LANGUAGE 'plpgsql';

※ plpgsqlを使用しているので、有効にしてない方はcreate languageを行ってください。
同時に2次メッシュ取得時に1次メッシュの西南端の緯度、経度が必要となるので、併せて関数を作成しておきます。

CREATE FUNCTION to_1st_mesh_geo(point) RETURNS point AS '
declare
    geo alias for $1;
    dest_geo point;
begin
    dest_geo = point( floor(geo[0] * 1.5) / 1.5, floor(geo[1]) );
    return dest_geo;
end; ' LANGUAGE 'plpgsql';

2次メッシュコードの取得

2次メッシュコードは、1次メッシュ西南端から、緯度方向に5分(0.08333)刻みで進んだ位置 || 経度方向に7分30秒(0.125)進んだ位置となるので、先ほどの作成した1次メッシュ西南端取得関数を利用して、

select
    floor( (lat - parent[0]) / 0.083333 )||
    floor( (lon - parent[1]) / 0.125 )
from (
    select
        to_1st_mesh_geo(point) as parent,
        point[0] as lat,
        point[1] as lon
    from landmark ) t1;
 ?column?
----------
 35
 42
 25

となりますね。こちらも関数登録しましょう。
2次メッシュ取得

CREATE FUNCTION to_2nd_mesh(point) RETURNS text AS '
declare
    geo alias for $1;
    left  text;
    right text;
    parent_point point;
    parent_mesh text;
begin
    select to_1st_mesh(geo) into parent_mesh;
    select to_1st_mesh_geo(geo) into parent_point;
    left  = cast( floor( (geo[0] -  parent_point[0]) / 0.083333 ) as text );
    right = cast( floor( (geo[1] -  parent_point[1]) /  0.125 ) as text );
    return parent_mesh || left || right;
end; ' LANGUAGE 'plpgsql'; 

2次メッシュ西南端緯度、経度取得

CREATE FUNCTION to_2nd_mesh_geo(point) RETURNS point AS '
declare
    geo alias for $1;
    parent_geo point;
    current_mesh text;
    dest_geo point;
begin
    select to_1st_mesh_geo( geo ) into parent_geo;
    select to_2nd_mesh( geo ) into current_mesh;
    dest_geo = point(
	parent_geo[0] + cast( substr( current_mesh, 5, 1 ) as float ) * 0.083333,
	parent_geo[1] + cast( substr( current_mesh, 6, 1 ) as float ) * 0.125
    );
    return dest_geo;
end; ' LANGUAGE 'plpgsql';

3次メッシュコードの取得

3次メッシュコードは、2次メッシュ西南端から、緯度方向に30秒(0.008333)刻みで進んだ位置 || 経度方向に45秒(0.0125)進んだ位置となります。
こちらも2次メッシュ同様に、上記で作成した関数を利用して、

select
    floor( (lat - parent[0]) / 0.0083333 )||
    floor( (lon - parent[1]) / 0.0125 )
from (
    select
        to_2nd_mesh_geo(point) as parent,
        point[0] as lat,
        point[1] as lon
    from landmark ) t1;
 ?column?
----------
 95
 78
 75

create function

CREATE FUNCTION to_3rd_mesh(point) RETURNS text AS '
declare
    geo alias for $1;
    parent_geo point;
    left  text;
    right text;
    parent_mesh  text;
begin
    select to_2nd_mesh(geo) into parent_mesh;
    select to_2nd_mesh_geo(geo) into parent_geo;
    left  = cast(floor( (geo[0] - parent_geo[0]) / 0.008333 ) as text);
    right = cast(floor( (geo[1] - parent_geo[1]) / 0.0125 ) as text);
    return parent_mesh || left || right;
end;' LANGUAGE 'plpgsql';

仕上げ

3次メッシュを取得

select to_3rd_mesh(point) from landmark;
 to_3rd_mesh
-------------
 53393595
 64414278
 39272575

データ登録

update landmark set mesh = to_3rd_mesh(point);
select * from landmark;
    name    |          point           |   mesh
------------+--------------------------+----------
 渋谷区役所 | (35.6640352,139.6982122) | 53393595
 札幌市役所 | (43.0620958,141.3543763) | 64414278
 沖縄市役所 | (26.2285302,127.6891103) | 39272575

特に何か幸せになれるというTipsではありませんが、位置ポイントを集約したりするのには、有効な策の1つだと思っています。
example

select to_3rd_mesh(point),count(point) from landmark group by to_3rd_mesh(point);
 to_3rd_mesh | count
-------------+-------
 53393595    |     2
 64414278    |     1
 39272575    |     1

以上となりますが、上記、8系バージョンでデバッグを行ったので、バージョンによっては、、、
来年はアップグレードの体験記書けるといいなあと思っております。

では、明日のPostgreSQL Advent Calender@DQNEOさんです。宜しくお願いします。