( ꒪⌓꒪) ゆるよろ日記

( ゚∀゚)o彡°オパーイ!オパーイ! ( ;゚皿゚)ノシΣ フィンギィィーーッ!!!

GeoHashのdecodeのアルゴリズムの解説します & ScalaのGeoHashライブラリを作ってみました(仮)

GeoHash(http://en.wikipedia.org/wiki/Geohash)は、緯度経度を文字列のハッシュで表現する仕様です。


GeoHashにより表現された緯度経度の情報は、一つの文字列で緯度と経度という2次元の情報に加えて精度も表すことができるという特徴を持っています。


例えば、どうでしょうバカの聖地である北海道札幌市の平岸高台公園は、北緯43.025東経141.377ですが、これをGeoHashで表現すると、"xpssc0"となります。

f:id:yuroyoro:20100115122758p:image


この"xpssc0"というGeoHash表現は、「北緯43.0224609375から43.0279541015625の間で、東経141.3720703125から141.383056640625の矩形範囲」であり、座標はこの矩形範囲の中心点になります。


@masuidrive blogさんの緯度経度を文字列で表すGeoHash - @masuidrive blogを読んでこの仕様を知り、面白そうだと思ってScala版のライブラリを作ってみました。GeoHashの解説もこちらの記事が非常に参考になります。

GeoHashをデコードするアルゴリズム

Wikipediaにも書いてありますが、文字列で表現されたGeoHashを実際の緯度経度にデコードするアルゴリズムを解説します。

1.Base32によるデコード

まず、GeoHashの文字列は、"BASE32"というエンコードが行われた結果の文字列です。
この"BASE32"は、0から32までの数値を、0から9までの数字と"a","i","l","o"を除いたアルファベットに変換するエンコードです。


具体的な変換テーブルは下記の通りです。

数値 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
Base32 0 1 2 3 4 5 6 7 8 9 b c d e f g

数値 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
Base32 h j k m n p q r s t u v w x y z


で、先ほどの"xpssc0"というGeoHashをBase32によりデコードすると、以下のような結果になります。

GeoHash(Base32) x p s s c 0
数値 31 21 24 24 11 0


次に、Base32によって変換された0から31までの数値を5bitのbit列の表現に変換します。

GeoHash(Base32) x p s s c 0
数値 31 21 24 24 11 0
bit列 11111 10101 11000 11000 01011 00000
2.デコードしたbit列から緯度と経度を分離させる

1でデコードした結果、"xpssc0"というGeoHashは、"11111 10101 11000 11000 01011 00000"というbit列の表現に変換されました。
このbit列から、緯度に対応するbitと経度に対応するbitを分離させます。


bit列の左から0を始点としてみて、偶数bitが経度、奇数bitが緯度に対応するbit列です。

bit番号 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
bit列 1 1 1 1 1 1 0 1 0 1 1 1 0 0 0 1 1 0 0 0
bit番号 20 21 22 23 24 25 26 27 28 29
bit列 0 1 0 1 1 0 0 0 0 0

色をつけたところが奇数なので緯度、色がないところは偶数なので緯度です。
ここから取り出してみると、緯度経度は以下のようなbit列として分離できます。

経度(偶数bit) 111001001000100
緯度(奇数bit) 101111010011000
3.緯度経度それぞれのbit列から範囲を算出する

さて、ここまで緯度と経度それぞれのbit列に変換できました。ここで、このbit列から具体的な緯度経度の数値を算出します。


算出の方法ですが、bit列に応じて範囲を2分割しながら絞り込む形をとります。
緯度の場合で説明します。


緯度の場合は、絞り込みの範囲が"-90度から90度"の範囲で開始します。


まず、この-90度から90度の範囲を2分割します。2分割すると、"-90から0"と"0から90"の二つの範囲になります。


与えられた緯度のbit列が、この2つの範囲のどちらに属するかは、bit列の1bit目が表しています。
1bit目は"1"ですので、大きい方の範囲である"0から90"の範囲の緯度であることが絞り込まれます。


次に、"0から90"の範囲を2分割します。2分割すると、"0から45"と"45から90"の二つの範囲になります。
2bit目は"0"ですので、小さい方の範囲である"0から45"の範囲の緯度であることが絞り込まれます。


この要領で、bit列の全てのbitに対して絞り込みを繰り返します。次のようなイメージです。

bit max mid min val
1 -90.000 0.000 90.000 45.000
0 0.000 45.000 90.000 22.500
1 0.000 22.500 45.000 33.750
1 22.500 33.750 45.000 39.375
1 33.750 39.375 45.000 42.188
1 39.375 42.188 45.000 43.594
0 42.188 43.594 45.000 42.891
1 42.188 42.891 43.594 43.242
0 42.891 43.242 43.594 43.066
0 42.891 43.066 43.242 42.979
1 42.891 42.979 43.066 43.022
1 42.979 43.022 43.066 43.044
0 43.022 43.044 43.066 43.033
0 43.022 43.033 43.044 43.028
0 43.022 43.028 43.033 43.025

最終的に、緯度は北緯43.025と算出されました。範囲は、"北緯43.022から北緯43.028"であることがわかります。


同じ要領で、経度のbit列に対しても"-180から180"を開始範囲として絞り込むこと経度が算出できます。

GeoHashをエンコード/デコードするScalaのライブラリ

で、書いてみました。


結構やっつけでエラーチェックとかしてないです。

object GeoHash{
  val base32 = "0123456789bcdefghjkmnpqrstuvwxyz"
  
  def encode( longtitude:Double, latitude:Double, range:Double) = {
    def asBit( d:Double,max:Double,min:Double ) : List[Boolean] = {
      val mid = ( max + min ) / 2
      var res = ( d >= mid )
      val (m,n) = if( res ) (max,mid) else (mid,min)
      if( Math.abs( m - n ) <= range ) 
        res :: Nil
      else
        res :: asBit( d ,m ,n )
    }
    
    def asBase32( xs:List[Boolean] ):String = {
      base32( (0.0 /: xs.reverse.zipWithIndex){ 
        case( n,(b,i) ) => if( b ) n + Math.pow( 2, i ) else n 
        } toInt
      ).toString
    }
    
    def encodeBits( xs:List[Boolean] ):String = xs match{
      case Nil => ""
      case _ =>
        val (h,t) = xs.splitAt(5)
        asBase32(h) + encodeBits( t )
    }
    
    val bits = asBit(longtitude , 180,-180 ).
      zipAll( asBit( latitude, 90, -90 ), false, false ).
      flatMap{ case (a,b) => a::b::Nil }
      
    encodeBits( bits ).reverse.dropWhile( '0' == ).reverse.mkString
  }

  def decode( geohash:String ):((Double,Double),(Double,Double)) = {
    def toBitList( s:String ) = s.toLowerCase.flatMap{ 
      c => ("00000" + base32.indexOf(c).toBinaryString ).
             reverse.take(5).reverse.map('1' == ) } toList
 
    def split( l:List[Boolean] ):(List[Boolean],List[Boolean]) ={
      l match{
        case Nil => (Nil,Nil)
        case x::Nil => ( x::Nil,Nil)
        case x::y::zs => val (xs,ys) =split( zs );( x::xs,y::ys)
      }
    }
 
    def dehash( xs:List[Boolean], min:Double, max:Double ):(Double,Double) = {
      ((min,max) /: xs ){ 
        case ((min,max) ,b) => 
          if( b )( (min + max )/2 , max )
          else ( min,(min + max )/ 2 )
       }
    }
    
    val ( xs ,ys ) = split( toBitList( geohash ) ) 
    ( dehash( xs, -180,180 ) , dehash( ys ,-90,90) )
  }
}