_geohash.lisp

;;;; -*- Mode: Lisp; Syntax: Common-Lisp -*-

;;;; GEOHASH encoder/decoder written by @nfunato on 2020-04-18
;;;; Copyright (c) 2020 Nobuhiko Funato, released under the MIT license

;;;; TODO: add utility functions, e.g. proximity judgment etc.

;;;; change from gist 20200419 version

(defpackage #:geohash
  (:use #:cl)
  (:export #:geohash-p #:geohash-encode #:geohash-decode))

(in-package #:geohash)

;;;
;;; external APIs
;;;

(defconstant +max-gh-len+ 12)

(defun geohash-encode (lat lng nbyte)
  "Encode LATitude and LoNGitude into a geohash string. NBYTE means precision, the length of the answer geohash string. The accuracy of geohash encoding depends on \"IEEE754 double-floatness\" of LAT and LNG, so it should be assured by the caller of this encoder."
  (assert (and (typep lat 'double-float) (<=  -90 lat  90)))
  (assert (and (typep lng 'double-float) (<= -180 lng 180)))
  (assert (and (integerp nbyte) (<= 1 nbyte +max-gh-len+)))
  (let ((nbit (* 5 nbyte))
        (shsm (shsm-to-put)))
    (flet ((mid-fn (i lo hi)
             (if (< i nbit)
                 (let* ((mi   (mid-coord lo hi))
                        (ltm? (< (if (evenp i) lng lat) mi)))
                   (put-bit shsm (if ltm? 0 1))
                   (values mi ltm?)))))
      (multiple-value-bind (w e s n)
          (walk-quadtree #'mid-fn -180.0d0 180.0d0 -90.0d0 90.0d0)
        (values (extract-string shsm)
                (list w e s n))))))

(defun geohash-decode (geohash)
  "Decode GEOHASH string, and return corresponding latitude and longitude."
  (assert (geohash-p geohash))
  (let ((nbit (* 5 (length geohash)))
        (shsm (shsm-to-get geohash)))
    (flet ((mid-fn (i lo hi)
             (if (< i nbit)
                 (let ((mi (mid-coord lo hi))
                       (b  (get-bit shsm)))
                   (values mi (zerop b))))))
      (multiple-value-bind (w e s n)
          (walk-quadtree #'mid-fn -180.0d0 180.0d0 -90.0d0 90.0d0)
        (let ((lat (mid-coord s n))
              (lng (mid-coord w e)))
          (values (list lat lng)
                  (list w e s n)))))))

(defun geohash-p (h)
  (and (base32-string-p h)
       (<= (length h) +max-gh-len+)))

;;;
;;; common helper functions
;;;

(defun mid-coord (lo hi)
  (/ (+ lo hi) 2))

(defun walk-quadtree (mid-fn w e s n)
  (labels ((walk (i low1 high1 low2 high2)
             (multiple-value-bind (mid lower-p) (funcall mid-fn i low1 high1)
               (cond ((null mid) (if (evenp i)
                                     (values low1 high1 low2 high2)
                                     (values low2 high2 low1 high1)))
                     (lower-p    (walk (1+ i) low2 high2 low1 mid))
                     (t          (walk (1+ i) low2 high2 mid  high1))))))
    (walk 0 w e s n)))

;; wrappers to mimic "string-holder"
(defun string-holder ()
  (make-array +max-gh-len+ :element-type 'character :fill-pointer 0))
(defun string-holder-to-string (s)
  (copy-seq s)) ; try to discard fill-pointer (works at least on SBCL)

;; string-holder state machine
(defstruct shsm str idx byte bitcnt)
(defun shsm-to-get (s) (make-shsm :idx -1 :bitcnt 4 :str s))
(defun shsm-to-put  () (make-shsm :byte 0 :bitcnt 0 :str (string-holder)))
(defun get-bit (shsm)
  (with-slots (str idx byte bitcnt) shsm
    (when (= 5 (incf bitcnt))
      (setf byte (base32-to-bitvec (char str (incf idx)))
            bitcnt 0))
    (bit byte bitcnt)))
(defun put-bit (shsm bit)
  (with-slots (str byte bitcnt) shsm
    (setf byte (+ (* 2 byte) bit))
    (when (= 5 (incf bitcnt))
      ;; idx is not used in put-op for fill-pointer works instead of it
      (vector-push (int-to-base32 byte) str)
      (setf byte 0
            bitcnt 0))
    (values)))
(defun extract-string (shsm)
  (string-holder-to-string (shsm-str shsm)))

;;;
;;; base32 tools
;;;

(defvar *base32-chars* "0123456789bcdefghjkmnpqrstuvwxyz")

(defun int-to-base32 (d)
  (char *base32-chars* d))

(defvar +base32-alpha-code+     ; "ailo" are omitted
  ;; a   b  c  d  e  f  g  h  i   j  k  l   m  n  o   p
  #(nil 10 11 12 13 14 15 16 nil 17 18 nil 19 20 nil 21
  ;; q  r  s  t  u  v  w  x  y  z
    22 23 24 25 26 27 28 29 30 31))

(defun base32-char-p (c)        ; same with (position c *base32-chars*)
  "Tests whether char is a base32-char. If it is, its weight is returned as an integer; otherwise nil is returned."
  (or (digit-char-p c)
      (and (lower-case-p c)
           (svref +base32-alpha-code+ (- (char-code c) (char-code #\a))))))

(defun base32-string-p (s)
  (and (stringp s)
       (every #'base32-char-p s)))

(defun base32-to-int (c)
  (or (base32-char-p c)
      (error "base32-to-int")))

(defvar +base32-bitvecs+
  #(#*00000 #*00001 #*00010 #*00011 #*00100 #*00101 #*00110 #*00111
    #*01000 #*01001 #*01010 #*01011 #*01100 #*01101 #*01110 #*01111
    #*10000 #*10001 #*10010 #*10011 #*10100 #*10101 #*10110 #*10111
    #*11000 #*11001 #*11010 #*11011 #*11100 #*11101 #*11110 #*11111))

(defun base32-to-bitvec (c)
  (svref +base32-bitvecs+ (base32-to-int c)))

;;;
;;; some tests
;;;

#+:hoge
(progn

(defmacro gh-test (sexpr) `(format t "~s => ~s~%" ',sexpr ,sexpr))

(defun test1 ()
  ;; Cribbed from github.com/papachan/geo_neighbors for reference
  (gh-test (equal (geohash-encode 35.7101389d0 139.8108333d0 8) "xn77jkz4"))
  (gh-test (equal (geohash-encode 4.616335d0 -74.071275d0 6) "d2g662"))
  (gh-test (equal (geohash-encode 35.7101389d0 139.8108333d0 3) "xn7"))
  (gh-test (equal (geohash-encode 35.7101389d0 139.8108333d0 10) "xn77jkz4ss"))
  (gh-test (equal (length (geohash-encode 35.7101389d0 139.8108333d0 3)) 3))
  (gh-test (equal (length (geohash-encode 35.7101389d0 139.8108333d0 8)) 8))
  (gh-test (equal (length (geohash-encode 35.7101389d0 139.8108333d0 10)) 10))
  )

(defun test2 ()
  ;; from wikipedia
  (gh-test (equal (geohash-encode 57.64911d0 10.40744d0 11) "u4pruydqqvj"))

  ;; The followings are from this code -- just placed here for regression test
  (gh-test (equal (geohash-decode "u4pruydqqvj")
                  '(57.64911063015461d0 10.407439693808556d0)))
  (gh-test (equal (geohash-decode "u4pruydqquv")
                  '(57.6491092890501d0  10.407439693808556d0)))
  (gh-test (equal (geohash-decode "ezs42")
                  '(42.60498046875d0 -5.60302734375d0)))
  )

;; (let ((*read-default-float-format* 'double-float))
;;   (defun test2.5 ()
;;     ;; same as test2 except superficial literal format
;;     (break "hoge ~s: ~s ~s" *read-default-float-format* 57.64911 57.64911d0)
;;     (gh-test (equal (geohash-encode 57.64911 10.40744 11) "u4pruydqqvj"))
;;     (gh-test (equal (geohash-decode "u4pruydqqvj")
;;                     '(57.64911063015461 10.407439693808556)))
;;     (gh-test (equal (geohash-decode "u4pruydqquv")
;;                     '(57.6491092890501  10.407439693808556)))
;;     (gh-test (equal (geohash-decode "ezs42")
;;                     '(42.60498046875 -5.60302734375))))
;;   )

) ; progn

geohash-20200419.lisp

;;;; -*- Mode: Lisp; Syntax: Common-Lisp -*-

;;;; GEOHASH encoder/decoder written by @nfunato on 2020-04-18
;;;; Copyright (c) 2020 Nobuhiko Funato, released under the MIT license

;;;; TODO: add utility functions, e.g. proximity judgment etc.

(defpackage #:geohash
  (:use #:cl)
  (:export #:geohash-p #:geohash-encode #:geohash-decode))

(in-package #:geohash)

;;;
;;; external APIs
;;;

(defun geohash-encode (lat lng nbyte)
  "Encode LATitude and LoNGitude into a geohash string. NBYTE means precision, the length of the answer geohash string. The accuracy of geohash encoding depends on \"IEEE754 double-floatness\" of LAT and LNG, so it should be assured by the caller of this encoder."
  (assert (and (typep lat 'double-float) (<=  -90 lat  90)))
  (assert (and (typep lng 'double-float) (<= -180 lng 180)))
  (assert (and (integerp nbyte) (<= 1 nbyte 12)))
  (let ((nbit (* 5 nbyte))
        (gh-holder (make-string-holder))
        (acc 0)
        (bitcnt 0))
    (flet ((mid-fn (i lo hi)
             (if (< i nbit)
                 (multiple-value-prog1
                     (let* ((mi (mid-coord lo hi))
                            (ltm? (< (if (evenp i) lng lat) mi)))
                       (setq acc (+ (* 2 acc) (if ltm? 0 1)))
                       (values mi ltm?))
                   (when (= 5 (incf bitcnt))
                     (vector-push (int-to-base32 acc) gh-holder)
                     (setq acc 0
                           bitcnt 0))))))
      (multiple-value-bind (w e s n)
          (walk-quadtree #'mid-fn -180.0d0 180.0d0 -90.0d0 90.0d0)
        (values (string-holder-to-string gh-holder)
                (list w e s n))))))

(defun geohash-decode (geohash)
  "Decode GEOHASH string, and return corresponding latitude and longitude."
  (assert (geohash-p geohash))
  (let ((nbit (* 5 (length geohash)))
        (gh-idx -1)
        bits
        (bitcnt 4))
    (flet ((mid-fn (i lo hi)
             (if (< i nbit)
                 (progn
                   (when (= 5 (incf bitcnt))
                     (setq bits (base32-to-bitvec (aref geohash (incf gh-idx)))
                           bitcnt 0))
                   (let ((mi (mid-coord lo hi))
                         (b  (elt bits bitcnt)))
                     (values mi (zerop b)))))))
      (multiple-value-bind (w e s n)
          (walk-quadtree #'mid-fn -180.0d0 180.0d0 -90.0d0 90.0d0)
        (let ((lat (mid-coord s n))
              (lng (mid-coord w e)))
          (values (list lat lng)
                  (list w e s n)))))))

(defun geohash-p (h)
  (and (base32-string-p h)
       (<= (length h) 12)))

;;;
;;; common helper functions
;;;

(defun mid-coord (lo hi)
  (/ (+ lo hi) 2))

(defun walk-quadtree (mid-fn w e s n)
  (labels ((walk (i low1 high1 low2 high2)
             (multiple-value-bind (mid lower-p) (funcall mid-fn i low1 high1)
               (cond ((null mid) (if (evenp i)
                                     (values low1 high1 low2 high2)
                                     (values low2 high2 low1 high1)))
                     (lower-p    (walk (1+ i) low2 high2 low1 mid))
                     (t          (walk (1+ i) low2 high2 mid  high1))))))
    (walk 0 w e s n)))

;; wrappers to mimic "string-holder"
(defconstant +max-gh-len+ 12)
(defun make-string-holder ()
  (make-array +max-gh-len+ :element-type 'character :fill-pointer 0))
(defun string-holder-to-string (s)
  (copy-seq s)) ; try to discard fill-poiner (works at least on SBCL)

;;;
;;; base32 tools
;;;

(defvar *base32-chars* "0123456789bcdefghjkmnpqrstuvwxyz")

(defun int-to-base32 (d)
  (aref *base32-chars* d))

(defvar +base32-alpha-code+     ; "ailo" are omitted
  ;; a   b  c  d  e  f  g  h  i   j  k  l   m  n  o   p
  #(nil 10 11 12 13 14 15 16 nil 17 18 nil 19 20 nil 21
  ;; q  r  s  t  u  v  w  x  y  z
    22 23 24 25 26 27 28 29 30 31))

(defun base32-char-p (c)        ; same with (position c *base32-chars*)
  "Tests whether char is a base32-char. If it is, its weight is returned as an integer; otherwise nil is returned."
  (or (digit-char-p c)
      (and (lower-case-p c)
           (svref +base32-alpha-code+ (- (char-code c) (char-code #\a))))))

(defun base32-string-p (s)
  (and (stringp s)
       (every #'base32-char-p s)))

(defun base32-to-int (c)
  (or (base32-char-p c)
      (error "base32-to-int")))

(defvar +base32-bitvecs+
  #(#*00000 #*00001 #*00010 #*00011 #*00100 #*00101 #*00110 #*00111
    #*01000 #*01001 #*01010 #*01011 #*01100 #*01101 #*01110 #*01111
    #*10000 #*10001 #*10010 #*10011 #*10100 #*10101 #*10110 #*10111
    #*11000 #*11001 #*11010 #*11011 #*11100 #*11101 #*11110 #*11111))

(defun base32-to-bitvec (c)
  (svref +base32-bitvecs+ (base32-to-int c)))

;;;
;;; some tests
;;;

#+:hoge
(progn

(defmacro gh-test (sexpr) `(format t "~s => ~s~%" ',sexpr ,sexpr))

(defun test1 ()
  ;; Cribbed from github.com/papachan/geo_neighbors for reference
  (gh-test (equal (geohash-encode 35.7101389d0 139.8108333d0 8) "xn77jkz4"))
  (gh-test (equal (geohash-encode 4.616335d0 -74.071275d0 6) "d2g662"))
  (gh-test (equal (geohash-encode 35.7101389d0 139.8108333d0 3) "xn7"))
  (gh-test (equal (geohash-encode 35.7101389d0 139.8108333d0 10) "xn77jkz4ss"))
  (gh-test (equal (length (geohash-encode 35.7101389d0 139.8108333d0 3)) 3))
  (gh-test (equal (length (geohash-encode 35.7101389d0 139.8108333d0 8)) 8))
  (gh-test (equal (length (geohash-encode 35.7101389d0 139.8108333d0 10)) 10))
  )

(defun test2 ()
  ;; from wikipedia
  (gh-test (equal (geohash-encode 57.64911d0 10.40744d0 11) "u4pruydqqvj"))

  ;; The followings are from this code -- just placed here for regression test
  (gh-test (equal (geohash-decode "u4pruydqqvj")
                  '(57.64911063015461d0 10.407439693808556d0)))
  (gh-test (equal (geohash-decode "u4pruydqquv")
                  '(57.6491092890501d0  10.407439693808556d0)))
  (gh-test (equal (geohash-decode "ezs42")
                  '(42.60498046875d0 -5.60302734375d0)))
  )

;; (let ((*read-default-float-format* 'double-float))
;;   (defun test2.5 ()
;;     ;; same as test2 except superficial literal format
;;     (break "hoge ~s: ~s ~s" *read-default-float-format* 57.64911 57.64911d0)
;;     (gh-test (equal (geohash-encode 57.64911 10.40744 11) "u4pruydqqvj"))
;;     (gh-test (equal (geohash-decode "u4pruydqqvj")
;;                     '(57.64911063015461 10.407439693808556)))
;;     (gh-test (equal (geohash-decode "u4pruydqquv")
;;                     '(57.6491092890501  10.407439693808556)))
;;     (gh-test (equal (geohash-decode "ezs42")
;;                     '(42.60498046875 -5.60302734375))))
;;   )

) ; progn