Geo::Hash->precision のよくわかる解説

PODを見ただけだと理解できなかったので、メモっておきます。

GeoHashについては、@masuidriveさんのエントリからリンクを辿ると色々わかると思います。

precision とはどういうメソッドか

PODによれば、適切なGeoHashの桁数を返してくれるということですが、これはどういうことでしょう。

Infer a suitable precision (number of character in hash) for a given lat, lon pair.

答えからいえば、渡した緯度経度の有効桁数を適切に表現できるGeoHashの桁数、となります。

例えば、経度を "100" と指定した場合に、 "100" と "101" という経度から同じGeoHashが生成されると困るわけです。この結果がきちんと分離されるために最低限必要なGeoHashの長さが、precisionの値となります。同様に、緯度を "45.00" で指定した場合は、経度 "45.00" と "45.01" によってそれぞれ違うGeoHashが生成されるために必要なGeoHashの長さ、となります。

実装を見てみる

precision の実装は以下です。マジックナンバーが散りばめられていて一目見ただけではよくわかりません。

# The number of bits necessary to represent the specified number of
# decimal digits
sub _d2b { int( shift() * 3.32192809488736 + 1 ) }

sub _bits_for_number {
    my $n = shift;
    return 0 unless $n =~ s/.*\.//;
    return _d2b( length $n );
}

sub precision {
    my ( $self, $lat, $lon ) = @_;
    my $lab = _bits_for_number( $lat ) + 8;
    my $lob = _bits_for_number( $lon ) + 9;
    return int( ( ( $lab > $lob ? $lab : $lob ) + 1 ) / 2.5 );
}

Geo/Hash.pm

順番にどのような処理が行われているのか考えてみましょう。まず、一番わかりやすいのは _bits_for_number です。これは、(最後に_d2bを噛ませている部分を除けば)文字列で表現された数字の小数点以下の有効桁数を求めています。"100"なら0桁、"100.00"なら2桁といった具合です。これをn桁とします。

さて、緯度と経度で二度計算を行いますが、内容は一緒なので緯度($lab)についてのみ考えましょう。緯度はGeoHashにおいては「緯度÷180」の形で表されますので、与えられた数値による最低の増分は "100" が与えられた場合は 1/180、 "100.00" が与えられた場合は 0.01/180 となります。一般には、 10^(-n)/180 となります。これをx桁の2進数で表せればいいので、「2^(-x) <= 10^(-n)/180」を満たす最小のxが、緯度をロストせずに表すために必要なビット数です。この式を2を底とする対数をとって整理すると、「x >= n * log_2(10) + log_2(180)」となります。

ここで、 log_2(10) = 3.32192809488736... であり、これがまず１つ目のマジックナンバーの正体です。また、log_2(180) = 7.49185309632967... となりますので与式は、「x >= n * 3.32192809488736 + 7.49185309632967」となります。ここで x が整数であることに着目して、「x >= ceil(n * 3.32192809488736 + 7.49185309632967)」がいえます。

現状のGeo::Hash の実装とこれを見比べてみると、ceilを右辺の1項と2項にそれぞれ展開して、「int( n * 3.32192809488736 + 1 )」と「8*1」とすれば尤もらしい式になります。ただし、この計算には些か不備があります。まず、ceilは分割できません*2。また、「ceil(x) ≠ int(x + 1)」です*3。

さて、ここまでで緯度と経度を表すのに必要なビット数 $lob と $lon が出ました。後は GeoHash が何桁あればこのビット数を表せるのかを求めれば、precision の値を求められます。GeoHash は32進数であり、1桁増える度に緯度と経度に使われるビット数は2ビット、または3ビット追加されます。この平均値が 2.5 なのでこの値で必要とされているビット数を割り、さらに整数値が欲しいので ceil してやれば、 precision() の最後の return に渡されている式が出てきます。

ただ、この実装もやはり不備があり、まず前述したように ceil は int では置き換えられません*4。後、桁数の算出に平均値の2.5を使っていますが、緯度と経度に分類し、きちんと計算しないと正しい値は出てきません*5。

まとめ

Geo::Hash->precision を使えば、渡した数字の有効桁数に適合する長さのGeoHashが自動的に生成されます。ただし、現行の実装には多少不備もありますので、「だいたいこの長さならいいよね」程度の精度だと思って利用すべきです。

なお、昨日Geo::Hash::XSにもこの関数を実装するパッチを送りましたが、現行ではあえてGeo::Hashと同じ実装にしています。

8/27 追記

間違ってるのを放置するのは忍びないのでパッチを送ってみました。

*1:経度の場合は log_2(360) == 8.49185309632968 を切り上げ。

*2:ceil(0.1+0.1) == 1 だが、ceil(0.1) + ceil(0.1) == 2

*3:ceil(1.0) == 1 だが、 int(1.0 + 1) == 2

*4:後、それならせめて「int( ( $lab > $lob ? $lab : $lob ) / 2.5 + 1 )」としないと、この実装だと経度が16bit必要な時に本当は7桁必要なのに結果は6桁となり、要件を満たさなくなってしまう。

*5:例えば、緯度が3ビットであれば必要なGeoHashは2桁。しかし、経度が3ビットであれば必要なGeoHashは1桁。