5
$\begingroup$

I'm having difficulty converting a geohash string to latitude and longitude with Mathematica and wondered if I'm missing a built in function that already has the functionality available?

(I'm still feeling my way with Mathematica and am not familiar with many of its capabilities or built in functionality so any advice or tips would be greatly welcomed)

$\endgroup$
3
$\begingroup$

Building upon FiddleMeRagged's solution, here's a revision that returns answers that match geohash.org:

narrowChoice[lower_, upper_, choice_] := 
  Module[{centre}, centre = (upper + lower)/2.0;
  If[choice == 0, {lower, centre}, {centre, upper}]];

relevantDigits[v_, err_] := Module[{digits = Round[-Log10[err] - 1]},
  NumberForm[v, {digits + RealDigits[v][[2]], digits}]];

decodeGeohash[s_] := 
 Module[{assoc32, bits, b, latitude, longitude, lowerLat, upperLat, 
   lowerLong, upperLong, actualLat, actualLong}, 
  assoc32 = <|"0" -> 0, "1" -> 1, "2" -> 2, "3" -> 3, "4" -> 4, 
    "5" -> 5, "6" -> 6, "7" -> 7, "8" -> 8, "9" -> 9, "b" -> 10, 
    "c" -> 11, "d" -> 12, "e" -> 13, "f" -> 14, "g" -> 15, "h" -> 16, 
    "j" -> 17, "k" -> 18, "m" -> 19, "n" -> 20, "p" -> 21, "q" -> 22, 
    "r" -> 23, "s" -> 24, "t" -> 25, "u" -> 26, "v" -> 27, "w" -> 28, 
    "x" -> 29, "y" -> 30, "z" -> 31|>;
  (*separate out the bits*)
  bits = Flatten[
    IntegerDigits[Lookup[assoc32, Characters[s]], 2, 5]];
  latitude = bits[[2 ;; ;; 2]];
  longitude = bits[[;; ;; 2]];
  lowerLat = -90; upperLat = 90;
  lowerLong = -180; upperLong = 180;
  (*for all the bits in latitude*)
  Do[{lowerLat, upperLat} = narrowChoice[lowerLat, upperLat, b], {b, 
    latitude}];
  (*for all the bits in longitude*)
  Do[{lowerLong, upperLong} = 
    narrowChoice[lowerLong, upperLong, b], {b, longitude}];
  actualLat = 
   relevantDigits[(lowerLat + upperLat)/2.0, 90/2^Length[latitude]];
  actualLong = 
   relevantDigits[(lowerLong + upperLong)/2.0, 
    180/2^Length[longitude]];
  {actualLat, actualLong}]

The (lowerLat+upperLat)/2.0 and (lowerLong+upperLong)/2.0 values are exact, but precision depends on the number of bits that were used to encode the values.

The number of digits to the right of the decimal point depends on the latitude and longitude errors. Each bit of latitude or longitude divides the error in half, so the latitude error is 90/2^(number of latitude bits) and the longitude error is 180/2^(number of longitude bits). relevantDigits[] accepts an exact value and error, and uses NumberForm[] to format the answer.

These answers agree nicely with geohash.org:

decodeGeohash/@{"9q8znb12j1", "mh7w", "0123e", "ezs42","t3b9m","c2b25ps", "80021bgm", "k484ht99h2", "8buh2w4pnt"}
(*{{37.79298,-122.3967},{-20.,50.},{-82.7,-179.5},{42.6,-5.6},{10.1,57.2},{49.26,-123.26},{0.005,-179.567},{-30.55555,0.2000},{5.00001,-140.6000}}*)
$\endgroup$
4
$\begingroup$

Here's a pretty compact implementation:

lrules = Thread[DeleteCases[CharacterRange["a", "z"], "a" | "i" | "l" | "o"] -> 
                CharacterRange["a", "v"]];

SetAttributes[geohash, Listable];
geohash[s_String] := Module[{bits, lal, lat, lol, lon},
        bits = IntegerDigits[FromDigits[StringReplace[s, lrules], 32],
                             2, 5 StringLength[s]];
        {lol, lal} = Flatten[Partition[bits, 2, 2, 1, {}], {2}];
        {lat, lon} = Mean /@ MapThread[Fold[ReplacePart[#1, #2 -> Mean[#1]] &, ##] &,
                                       {{{-90, 90}, {-180, 180}}, Mod[{lal, lol}, 2, 1]}];
        GeoPosition[N[MapThread[Round[#, 10^Ceiling[RealExponent[#2]]] &,
                                {{lat, lon}, {90, 180}/2^{Length[lal], Length[lol]}}]]]]

Using creidhne's example:

geohash[{"9q8znb12j1", "mh7w", "0123e", "ezs42", "t3b9m", "c2b25ps",
         "80021bgm", "k484ht99h2", "8buh2w4pnt"}]
   {GeoPosition[{37.793, -122.397}], GeoPosition[{-20., 50.}],
    GeoPosition[{-82.7, -179.5}], GeoPosition[{42.6, -5.6}],
    GeoPosition[{10.1, 57.2}], GeoPosition[{49.261, -123.26}], 
    GeoPosition[{0.0051, -179.567}], GeoPosition[{-30.5556, 0.2}],
    GeoPosition[{5.00001, -140.6}]}
$\endgroup$
  • $\begingroup$ This is a big improvement. However, there's a rounding problem here if the results are to match geohash.org. For example, "9q8znb12j1" should round to 5 and 4 places respectively for latitude and longitude. Instead, the result is rounded to 5 places for both. Replacing 10^Ceiling[RealExponent[#2]] with 10^(Round[RealExponent[#2]] + 1) seems to fix the problem. Note that even with this change, "t3b9m" returns {10.1,57.20000000000001}, but that difference gets lost in the noise. $\endgroup$ – creidhne Oct 3 '17 at 2:02
3
$\begingroup$

The GeohashUtils class from the Spatial4j library can do this, and querying a java library from Mathematica is pretty painless.

Download the jar file to some location on your hard drive, here I put it in ~/Downloads/spatial/.

<< JLink`

AddToClassPath["~/Downloads/spatial"];

We need to directly load 2 classes,

LoadJavaClass["org.locationtech.spatial4j.io.GeohashUtils"];
LoadJavaClass["org.locationtech.spatial4j.context.SpatialContext"];

The functions we need are static, so it isn't necessary to create new objects.

geohashToGeoPosition[hash_String] := 
 Module[{point = 
    GeohashUtils`decode["u4pruydqqvj", SpatialContext`GEO]},
  GeoPosition[{point@getY[], point@getX[]}]
  ]

And it can be tested via

geohashToGeoPosition@"dp1hpq9p"
(* GeoPosition[{40.1151, -88.2736}] *)

GeoListPlot[geohashToGeoPosition@"dp1hpq9p", GeoZoomLevel -> 5]

Mathematica graphics

You could generate geohashes with the same class,

geoPositionToGeohash[GeoPosition[{lat_, lon_, h_: 0}]] := GeohashUtils`encodeLatLon[lat, lon]


geoPositionToGeohash@
 GeoPosition[{40.1151180267334`, -88.27360153198242`}]
geohashToGeoPosition@%
(* "dp1hpq9p7zzz" *)
(* GeoPosition[{40.1151, -88.2736}] *)

But you may need to tweak this to deal with GeoPosition objects with specified datums.

$\endgroup$
  • $\begingroup$ Many thanks, this works perfectly though I'm surprised there's not a native function at the moment. $\endgroup$ – FractalDoctor Sep 28 '17 at 8:22
2
$\begingroup$

So, I had the time to look into this further and took the opportunity to try to write a Mathematica function to decode a geohash given the description of the algorithm on Wikipedia.

    (* decoding geohashes *)

narrowChoice[lower_, upper_, choice_] := 
  Module[{centre},
   centre = (upper + lower) / 2.0;
   If[choice == 0, {lower, centre}, {centre, upper}]
   ];

    decodeGeohash[s_] :=
     Module[{assoc32},
      assoc32 = <| "0" -> 0, "1" -> 1, "2" -> 2, "3" -> 3, "4" -> 4, 
        "5" -> 5, "6" -> 6, "7" -> 7, "8" -> 8, "9" -> 9, "b" -> 10, 
        "c" -> 11, "d" -> 12, "e" -> 13, "f" -> 14, "g" -> 15, "h" -> 16, 
        "j" -> 17, "k" -> 18, "m" -> 19, "n" -> 20, "p" -> 21, "q" -> 22, 
        "r" -> 23, "s" -> 24, "t" -> 25, "u" -> 26, "v" -> 27, "w" -> 28, 
        "x" -> 29, "y" -> 30, "z" -> 31|>;
      bits = Flatten[
        IntegerDigits[Lookup[assoc32, Characters[s]], 2, 5]];
      (* separate out the bits *)
      latitude = bits[[2 ;; ;; 2]];
      longitude = bits[[;; ;; 2]];
      (* for all the bits in latitude *)
      lowerLat = -90;    
      upperLat = 90;
      lowerLong = -180; upperLong = 180;
      Do[{lowerLat, upperLat}     = 
        narrowChoice[lowerLat,   upperLat,  
         latitude [[element]]],  {element, 1, Length[latitude], 1}];
      Do[{lowerLong, upperLong} = 
        narrowChoice[lowerLong, upperLong, 
         longitude[[element]]], {element, 1, Length[longitude], 1}];
      actualLat = (lowerLat + upperLat) /2.0;
      actualLong = (lowerLong + upperLong)/2.0;
      {actualLat, actualLong}
      ]

Taking the geodes of the San Francisco Embarcadero subway as 9q8znb12j1 I get an answer of {37.793, -122.397} compared to geohash.org which gives {37.79298 -122.3967}.

I'm probably taking liberties taking the average of my final result, and I'm sure there's a more functional way than with the two Do loops (perhaps a way with Map, but I couldn't figure out how to feed back in the previous results). Any comments or improvements are welcome.

$\endgroup$
  • $\begingroup$ You can likely use Fold or Nest instead of Do $\endgroup$ – b3m2a1 Oct 1 '17 at 18:20
2
$\begingroup$

Here's a little Compiled version of @J.M.'s code:

geohashDecode[s_String] :=

  GeoPosition@geohashDecodeCore@ToCharacterCode[s];
geohashDecode~SetAttributes~Listable

geohashDecodeCore =
  Compile[{{chars, _Integer, 1}},
   Module[{base32digits, base32int, bits, lal, lat, lol, lon},
    base32digits =
     Map[
      Which[
        48 <= # <= 57, # - 48,
        98 <= # <= 104, # - 88,
        106 <= # <= 107, # - 89,
        109 <= # <= 110, # - 90,
        112 <= # <= 122, # - 91,
        True, -1
        ] &,
      chars
      ];
    base32int = 
     Apply[Plus, 
      MapIndexed[#*32^(#2[[1]] - 1) &, Reverse@base32digits]];
    bits = IntegerDigits[base32int, 2, 5 Length[chars]];
    lol = bits[[Range[1, Length[bits], 2]]];
    lal = bits[[Range[2, Length[bits], 2]]];
    lat = {-90., 90.};
    Scan[
     Function[Set[lat[[#]], Apply[Plus, lat]/2.]],
     Mod[lal, 2, 1]
     ];
    lon = {-180., 180.};
    Scan[Set[lon[[#]], Apply[Plus, lon]/2.] &,
     Mod[lol, 2, 1]
     ];
    MapThread[
     Round,
     {
      {Apply[Plus, lat]/2, Apply[Plus, lon]/2},
      10.^Ceiling@Log10[{90., 180.}/2^{Length[lal], Length[lol]}]
      }
     ]
    ]
   ];

It's not really necessary here, since J.M.'s code is so fast already:

geohashDecJM@{"9q8znb12j1", "mh7w", "0123e", "ezs42", "t3b9m", 
    "c2b25ps", "80021bgm", "k484ht99h2", "8buh2w4pnt"} // 
  RepeatedTiming // First
geohashDecode[{"9q8znb12j1", "mh7w", "0123e", "ezs42", "t3b9m", 
    "c2b25ps", "80021bgm", "k484ht99h2", "8buh2w4pnt"}] // 
  RepeatedTiming // First

0.0034

0.00045

But maybe if you have a huge number of hashes it'll be helpful. In any case it's kinda fun.

$\endgroup$
  • $\begingroup$ It's a bit better to use Scan[] instead of Map[], since you don't really care about the output of Set[] for lat and lon. $\endgroup$ – J. M. will be back soon Oct 1 '17 at 20:40
  • $\begingroup$ @J.M. Oh true, I had Map from an old way of going about it. I'll tweak that. $\endgroup$ – b3m2a1 Oct 1 '17 at 20:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.