Jun 20 2008

Randomised walks across your neighbourhood: geohashing

Published by Dougal at 5:30 pm under Computing, Humour, Programming

A few weeks ago Randall Munroe published XKCD 426: Geohashing and apparently invented a new sport…

The idea is quite simple. The world is divided into little degree-by-degree segments (by latitude and longitude). If you see these on a map they look mostly rectangular at the equator and gradually get more triangular at the poles, because the pole-wards side is shorter than the equator side. You can see a picture of the geographical segment around Edinburgh on the wiki.

If you imagine that each rectangle has a starting corner (the one nearest the equator and nearest the Greenwich meridian), which we’ll call (0.0, 0.0). We can identify any spot in your rectangle with a fractional offset from this point — like (0.456, 0.235).

If you want to know the major co-ordinates for your home then a good place to start would be this list for the major cities of the countries of the world.

I’ve put together a Haskell program to demonstrate the next stage of the procedure, though there are plenty of web-based tools to do the same thing. I just wanted to try out the new cabal-install package (akin to CPAN, gems etc for other languages).

module Main where
 
import Network.Curl
import System.Time
import System.Locale
import Data.Digest.MD5
import Data.Char
import Numeric
import Codec.Utils
import System.Environment
import Text.Printf
import Control.Monad
 
dowjonesurl = "http://irc.peeron.com/xkcd/map/data/%Y/%m/%d"
ymd = "%Y-%m-%d"
usage = "Usage: geohash <latitude> <longitude>"
 
getTimes adj = do
    time <- getClockTime
    today <- fmtcal ymd `liftM` toCalendarTime time
    let adjustment = noTimeDiff { tdDay = if adj then -1 else 0 }
    stockdate <- toCalendarTime $ addToClockTime adjustment time
    return (today, fmtcal dowjonesurl stockdate)
  where fmtcal = formatCalendarTime defaultTimeLocale
 
main = do
    args <- getArgs
    let [lat, lon] = case args of
                      [a,b] -> map read args
                      _     -> error usage
 
    (today, djia) <- getTimes (lat > -30)
 
    -- Get DJIA data, result not checked
    (_, avg) <- curlGetString djia []
 
    let inputstring = today ++ "-" ++ avg
        (latfrac, lonfrac) = calcFractions $ md5hash inputstring
    printf "%.6f %.6f\n" (lat .+ latfrac) (lon .+ lonfrac)
 
stringToOctets :: String -> [Octet]
stringToOctets = map (toEnum . ord)
 
calcFractions :: [Octet] -> (Double, Double)
calcFractions = both (hexfraction . concatMap nibbles) . splitAt 8
 
md5hash = hash . stringToOctets
hexfraction = foldr (\x s -> (x + s)/16) 0
 
nibbles :: Octet -> [Double]
nibbles o = let (d,r) = o `quotRem` 16 
            in [fromIntegral d, fromIntegral r]
 
-- Utils
(.+) :: Double -> Double -> Double
(.+) num part = let whole = fromIntegral $ truncate num
                in signum whole * (abs whole + abs part)
 
both f (a,b) = (f a, f b)

The result is just a simple latitude and longitude pair, which you can drop straight into Google Maps’ text box.

But as I said, the real purpose behind this was using cabal-install, which was marvellous. I pulled in the libcurl bindings for easy network access and the Crypto package for the MD5 hashes without any bother at all. Well done to all the hackers who’ve been working hard on Cabal and its associated tools, and on the contributors to Hackage.

Now it’s as simple as cabal install Crypto to download, build and install all dependencies for whatever package you want. Bravo!

Update: The original program was wrong as I did not think to test it against the reference implementation. I had just been using the decimal hash as the fractional part of the co-ordinates. In fact what I had to do was turn the hash into a hexadecimal fraction, ie each step to the right of the point is 1/16 of its neighbour. I’ve updated (and simplified) the program accordingly.

Trackback URI | Comments RSS

Leave a Reply