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

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

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

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

この"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) )
  }
}