Gaucheでgeohash できた

Gaucheでgeohashする一連の投稿だけど、やっとライブラリが完成した。

モジュール化もやってみたよ。
(use geohash)
で使える。

関数はふたつ。

geohash-encode LATITUDE LONGITUDE (:precision 9)
  => geohash文字列

geohash-decode GEOHASH
  => LATITUDE LONGITUDE LAT下限 LAT上限 LON下限 LON上限

デコードのときはgeohashの文字列の長さで有効桁が長くなるので、返した座標値のほかにその座標値がどれくらいの範囲まで正確なのかの値も返してる。
多値を返せるって便利だなあ。

#!/usr/bin/env gosh
;; Copyright (c) 2008 KOGA Kazuo
;; MIT License
;;
;; Geohash
;; see http://en.wikipedia.org/wiki/Geohash
;;

(define-module geohash
  (use srfi-1)
  (use gauche.sequence)
  (use gauche.collection)
  (use gauche.uvector)
  (export geohash-encode
          geohash-decode))
(select-module geohash)

(define (medium x y) (/ (+ x y) 2))

(define (interval-choice which min max)
  (if (= 0 which)
      (values min (medium min max))
      (values (medium min max) max)))

(define (interval-fold code min max)
  (fold2 interval-choice min max
         code))

(define (interval-fold-fold codes min max)
  (fold2 interval-fold min max
         codes))

(define (interval-bit n min max)
  (let1 mid (medium min max)
    (if (< n mid)
        (values 0 min mid)
        (values 1 mid max))))

(define base32-base
  #(#\0 #\1 #\2 #\3 #\4 #\5 #\6 #\7
    #\8 #\9 #\b #\c #\d #\e #\f #\g
    #\h #\j #\k #\m #\n #\p #\q #\r
    #\s #\t #\u #\v #\w #\x #\y #\z))

(define base32-bits
  #(#u8(0 0 0 0 0) #u8(0 0 0 0 1) #u8(0 0 0 1 0) #u8(0 0 0 1 1)
       #u8(0 0 1 0 0) #u8(0 0 1 0 1) #u8(0 0 1 1 0) #u8(0 0 1 1 1)
       #u8(0 1 0 0 0) #u8(0 1 0 0 1) #u8(0 1 0 1 0) #u8(0 1 0 1 1)
       #u8(0 1 1 0 0) #u8(0 1 1 0 1) #u8(0 1 1 1 0) #u8(0 1 1 1 1)
       #u8(1 0 0 0 0) #u8(1 0 0 0 1) #u8(1 0 0 1 0) #u8(1 0 0 1 1)
       #u8(1 0 1 0 0) #u8(1 0 1 0 1) #u8(1 0 1 1 0) #u8(1 0 1 1 1)
       #u8(1 1 0 0 0) #u8(1 1 0 0 1) #u8(1 1 0 1 0) #u8(1 1 0 1 1)
       #u8(1 1 1 0 0) #u8(1 1 1 0 1) #u8(1 1 1 1 0) #u8(1 1 1 1 1)
       ))

(define (pos-even?)
  (let1 even #f
    (lambda (_)
      (set! even (not even))
      even)))

(define (decode-base32 s)
  (map (lambda (c)
         (ref base32-bits
              (find-index (cut eqv? c <>) base32-base)))
       s))

(define (geohash-decode str)
  "geohash => latitude longitude latitude-min latitude-max longitude-min longitude-max"
  (let ((bits (decode-base32 str))
        (pos (pos-even?)))
    (let loop ((b bits)
               (lon-min -180.0)
               (lon-max 180.0)
               (lat-min -90.0)
               (lat-max 90.0))
      (if (null? b)
          (values (medium lat-min lat-max)
                  (medium lon-min lon-max)
                  lat-min lat-max
                  lon-min lon-max)
          (receive (lon lat)
              (partition pos (car b))
            (receive (lomin lomax)
                (interval-fold lon lon-min lon-max)
              (receive (lamin lamax)
                  (interval-fold lat lat-min lat-max)
                (loop (cdr b)
                      lomin lomax lamin lamax))))))))

(define (bit->number b0 b1 b2 b3 b4)
  (+ (* b0 16)
     (* b1 8)
     (* b2 4)
     (* b3 2)
     b4))

(define (encode longitude latitude precision)
  (define (req-bits n)
    (receive (q r) (quotient&remainder n 2)
      (+ (* n 2) q r)))
  (define (morton-bits num min max count)
    (let loop ((i 0)
               (min min)
               (max max)
               (ret '()))
      (if (< i count)
          (receive (bit min max) (interval-bit num min max)
            (loop (+ i 1) min max (cons bit ret)))
          (reverse! ret))))
  (define (composer a b)
    (values (bit->number (ref a 0) (ref b 0)
                         (ref a 1) (ref b 1)
                         (ref a 2))
            (drop a 3) (drop b 2)))

  (let* ((need-bits (req-bits precision))
         (lon-bits (morton-bits longitude -180 180 need-bits))
         (lat-bits (morton-bits latitude -90 90 need-bits)))
    (let loop ((lon lon-bits)
               (lat lat-bits)
               (even #t)
               (ret '()))
      (if (null? lon)
          (map-to <string>
                  (lambda (e)
                    (ref base32-base e))
                  (reverse! ret))
          (receive (n a b)
              (if even
                  (composer lon lat)
                  (composer lat lon))
            (if even
                (loop a b (not even) (cons n ret))
                (loop b a (not even) (cons n ret))))))))

(define (geohash-encode latitude longitude . opts)
  "latitude longitude => geohash"
  (unless (and (<= -90 latitude 90)
               (<= -180 longitude 180))
    (error "geohash-encode: bad domain: latitude or longitude"))
  (let-keywords opts ((precision 9))
    (encode longitude latitude precision)))

(provide "geohash")