Geohashを作成するコードを書いてみた


こんにちは

今日はいつもと趣向が違います。
このところやんごとなき理由により地理情報を扱うプログラムを書いています。
今ちょうどデータ変換中でちょっと時間が空いてしまったので、みなさんに共有できそうなコードを載せてみました。

それにたまにはちゃんとエンジニアだっていうことを証明するために技術的な話題もしないと。

タイトル通り、Geohashです。
知らなかったのですが、Wikipediaにある通り、色々と使い道がありそうです。
ジオハッシュ – Wikipedia

データベースにこのハッシュコードを入れて

select * from hoge where geohash like 'xxxxx%';

とやれば、通常の昇順インデックスをgeocodeに張ってあげるだけで近隣のオブジェクトをそれなりの速さで取ってこれるというなんとも便利なコードです。

うーん、使いたい!

というわけでコードを書いてみました。

といっても、Wikipediaに載っているアルゴリズムをそのまま実装するのは微妙なのでインチキしてみました。
要は、半分の半分の、、、というところで2進数表記のことだと気付いたので、それなら緯度と経度をそれぞれの値域の幅、それぞれ180と360で割ってしまって、そうすると-0.5〜0.5の間に収まるのであとは0.5を足せば、なんか良さそうだなと思います。
イメージ的には、緯度や経度が、

 0.10010010

とかなる感じですね。

あとは小数点以下をうまいこと取り出してごにょごにょすれば目的達成です。
pack/unpackを使ってビット列を取り出すことにします。

というわけで内部表現を調べることにします。
ここからは本当はアーキテクチャーによってビット表現が違うとは思うけど、割り切ってダークサイドの移植性0のコードを書きます。
いま使っているサーバーがIA-64なので、ビックエンデイアンでIEEE754形式なんだろうなと、割り切ります。

では早速、仕様をWikipediaに聞いてみます。
IEEE 754 – Wikipedia

なるほど、指数表現なんですね。浮動小数っていうぐらいだから当たり前か。
そうすると、あたまの数字に0の数によってビット列が左に動いてしまいます。例えば、

0.1001 x 2^0 = 1.001 x 2^-1
0.0001 x 2^0 = 1.000 x 2^-4

みたいな感じ。
それでは困ってしまうので、1を足してしまいます。そうすれば頭の重みが固定されます。例えば、

0.1001 -> 1.1001
0.0001 -> 1.0001

みたいな感じ。これなら頭が固定されるので、指数部が揃います。
というわけで、さきほど、180で割って0.5を足すと言いましたが、1.5を足すことにします。

あとは13ビット目から最後まで、52ビットを取ってくれば目的物がとれます。
あたまの1はケチ表記で省略されているそうなので全部持ってくることにします。

戻す時は、符号は0で、指数部が0なんだけどバイアスの1023を気にして、頭の部分に

001111111111

をつけてあげます。

動かしてみたところ、まあそれなりに動いているようです。
実際にはハッシュとして使いたいわけではないので、BASE64でエンコードしたほうが該当するエリアが正方形になるので使いやすいです。
あるいは、BASE16(なんてあるのか?)でエンコードして細かく精度をコントロールするのもアリかと思います。

最後にコードを載せておきます。
決して良いコードとは言えないけど、ノリでこういうのがかけちゃうのはRubyの良いところということで、お茶を濁しておきます。

module GeoHash
  # 緯度経度オブジェクト
  class LatLng

    attr_accessor :lat, :lng
    # コンストラクタ
    def initialize(lat, lng)
      @lat = lat
      @lng = lng
    end

    # ハッシュへ変換する
    def to_hash
      zip(f2b(lng / 360 + 1.5, f2b(lat / 180 + 1.5)).
        scan(/...../). # 5文字ずつ取ってくる
        map { |b| ('0b' + b).oct }. # 文字列を数値へ変換
        map { |i| CHAR_MAP[i] }.
        join
    end

    # ハッシュを緯度経度オブジェクトへ変換する
    def self.[](str)
      str += '0' if str.length % 2 == 1
      bin = unzip str.each_char.
        map { |c| CHAR_MAP.index c }.
        map { |i| '%05b' % i }.
        join

      LatLng.new (b2f(bin[1) - 1.5) * 180, (b2f(bin[0]) - 1.5) * 360
    end

  private
    # ハッシュに使用する文字
    CHAR_MAP = '0123456789bcdefghjkmnpqrstuvwxyz'

    # 二つのビット列を一つにまとめる
    def zip(s1, s2)
      s1.each_char.zip(s2.each_char).map(&:join).join
    end

    # ビット列を2つに分ける
    # 2文字づつへの分け方しかわからないので、転置行列をとる
    def unzip(s)
      s.each_char.each_slice(2).to_a.transpose.map(&:join)
    end

    # 浮動小数から二進数へ変換する
    def f2b(f)
      [f].pack('d').
         reverse. # ビッグエンディアンのみ動作
         unpack('B64')[0][12, 52] + '000' # IEEE754のみ動作
    end

    # 二進数から浮動小数へ変換する
    def b2f(b)
      [('001111111111' + b + ('0' * 52))[0, 64]]. # IEEE754のみ動作
        pack('B64').
        reverse. # ビッグエンディアンのみ動作
        unpack('d')[0]
    end
  end
end

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です