Login

Factual Blog /

How Geohashes Work

We use geohashes all the time at Factual, so one day I got curious and read through the canonical Java implementation to figure out exactly how they work. It was an enlightening read, but along the way I encountered some unfortunate bits of coding horror that got me wondering about the fastest way to encode and decode geohashes.

The geohash format
A geohash is a series of bits that repeatedly bisects a search space of latitudes and longitudes. The first bit bisects the longitude, the next one latitude, the next longitude, etc. As a result, some geohashes specify square-ish (in spherical coordinates) regions, while others specify 2:1 rectangles. Because geohashes are fundamentally binary data, they’re often written in base-32. The alphabet is a little unusual in that it omits some letters that might be confused with numbers.

Like floating-point numbers, geohashes contain some error depending on how many bits they have. A geohash library often returns the range of latitude and longitude values you get when decoding, since geohash error is typically much larger than the underlying floating-point error. You can’t technically decode a geohash to a specific point, but you can usually approximate by taking the midpoint of the error box. Wikipedia has a great example that illustrates how this works.

Typical encoding strategy
Geohashing is basically storing the intermediate results of two binary searches. The basic idea looks something like this in code:

long encodeGeohash(double lat, double lng, int bits) {
  double minLat = -90,  maxLat = 90;
  double minLng = -180, maxLng = 180;
  long result = 0;
  for (int i = 0; i < bits; ++i)
	if (i % 2 == 0) {                 // even bit: bisect longitude
	  double midpoint = (minLng + maxLng) / 2;
	  if (lng < midpoint) {
		result <<= 1;                 // push a zero bit
		maxLng = midpoint;            // shrink range downwards
	  } else {
		result = result << 1 | 1;     // push a one bit
		minLng = midpoint;            // shrink range upwards
	  }
	} else {                          // odd bit: bisect latitude
	  double midpoint = (minLat + maxLat) / 2;
	  if (lat < midpoint) {
		result <<= 1;                 // push a zero bit
		maxLat = midpoint;            // shrink range downwards
	  } else {
		result = result << 1 | 1;     // push a one bit
		minLat = midpoint;            // shrink range upwards
	  }
	}
  return result;
}

Decoding is similar; instead of bisecting based on midpoint comparisons, you just go backwards by setting either the min or max value to the midpoint based on the bits from the geohash.

Optimizing the algorithm
The code above is already pretty fast. We aren’t allocating any memory or making any method calls, so it’s a handful of primitive operations and conditions per bit we want to encode. The only way to make it substantially faster is to get rid of the loop; this involves finding a sublinear way to do both the bisections and the interleaving.

The latitude and longitude are in some sense bisected already since they’re specified as doubles. The only problem is that the mantissa bisects the range between two powers of two, not the limits of latitude or longitude. But we can change that easily enough with an add and a divide:

double adjustedLat = (lat + 90) / 180.0;
double adjustedLng = (lng + 180) / 360.0;

Now our numbers are both in the range [0, 1]. We can convert them to bits either by using Double.doubleToLongBits, or by multiplying and doing a cast:

long latBits = (long) (adjustedLat * 0x80000000l) & 0x7fffffffl;
long lngBits = (long) (adjustedLng * 0x80000000l) & 0x7fffffffl;

At this point we have the two binary search results, each stored in the low 31 bits of the respective longs. (Technically 30 is all we need since most geohashes end up written in base-32.) Now we need to interleave them, also without looping.

It turns out that this is a well-studied problem with some very cool solutions. My favorite one is from the indispensable Stanford bithacks website:

static const unsigned int B[] = {0x55555555, 0x33333333, 0x0F0F0F0F, 0x00FF00FF};
static const unsigned int S[] = {1, 2, 4, 8};
unsigned int x; // Interleave lower 16 bits of x and y, so the bits of x
unsigned int y; // are in the even positions and bits from y in the odd;
unsigned int z; // z gets the resulting 32-bit Morton Number.  
				// x and y must initially be less than 65536.
x = (x | (x << S[3])) & B[3];
x = (x | (x << S[2])) & B[2];
x = (x | (x << S[1])) & B[1];
x = (x | (x << S[0])) & B[0];
y = (y | (y << S[3])) & B[3];
y = (y | (y << S[2])) & B[2];
y = (y | (y << S[1])) & B[1];
y = (y | (y << S[0])) & B[0];
z = x | (y << 1);

We need a 64-bit Morton number, so we’ll need one more stage. Here’s the code I ended up with to gap the bits:

static long widen(long low32) {
  low32 |= low32 << 16; low32 &= 0x0000ffff0000ffffl;
  low32 |= low32 << 8;  low32 &= 0x00ff00ff00ff00ffl;
  low32 |= low32 << 4;  low32 &= 0x0f0f0f0f0f0f0f0fl;
  low32 |= low32 << 2;  low32 &= 0x3333333333333333l;
  low32 |= low32 << 1;  low32 &= 0x5555555555555555l;
  return low32;
}

Then to finalize:

widen(lngBits) | widen(latBits) >>> 1;

This result can then be right-shifted to get the desired precision.

Calculating error bounds (decoding)
I mentioned earlier that geohash libraries typically return bounding boxes when decoding. Normally we’d have the required information available after the decoding loop, but the loopless version returns results without actually constructing the lat/lng ranges. We could trivially store these ranges in a lookup table since there are only 64 levels of precision, but there’s a fun algorithmic solution that’s almost as fast.

The latitude and longitude error are related: if #bits is even then longitude error = 2x latitude error, otherwise they’re equal. Given that, it suffices to calculate the latitude error first and then figure out the longitude error. Only every other bit impacts the latitude error, so there are 32 distinct possibilities. Specifically, the latitude error is 180 * 2-floor(#bits/2).

Math.pow() is the obvious solution, but it’s boring and is generally considered slow. Given that our exponent is always a 5-bit integer, we can use repeated squaring to very quickly come up with an exact answer:

double latError = 1;
for (int b = 0; b < 5; ++b) {
  latError *= latError;
  if ((bits & 2 << b) != 0)       // int bits is the geohash precision (0 to 63, inclusive)
	latError *= 0.5;
}
latError *= 90;

This can be optimized a bit further by unrolling and saving a couple of multiplies:

double error = 1.0;
	if ((bits & 32) != 0) error *= 0.25;
	if ((bits & 16) != 0) error *= 0.5; error *= error;
	if ((bits & 8)  != 0) error *= 0.5; error *= error;
	if ((bits & 4)  != 0) error *= 0.5; error *= error;
	if ((bits & 2)  != 0) error *= 0.5;

	double latError = error * 90;
	double lngError = error * ((bits & 1) != 0 ? 90 : 180);

Final notes: decoding base-32
A common strategy is to use a HashMap or similar to do reverse-lookups on characters. That’s how the reference implementation does it, for example. The problem is that Java generic data structures always work on boxed values, so every invocation is going to require allocating the boxed type, calculating its hashcode, doing a .equals() check, and then unboxing the result. Not only is this a lot of work at runtime, but the lookup structure itself is quite large due to the nontrivial per-object overhead imposed by the JVM.

Even though it seems wasteful, a primitive array is a much more efficient way to represent the reverse lookup table:

static final char[] BASE32 =
	  { '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' };

	static final byte[] BASE32_INV = new byte[(int) 'z' + 1];

	static {
	  for (int i = 0; i < BASE32.length; ++i)
	    BASE32_INV[(int) BASE32[i]] = (byte) i;
	}

Check out the code on Github
This code is available in the Factual beercode-open repository, which is a collection of libraries whose correctness is informally guaranteed by beer. And if Morton codes and repeated squaring are your thing, you might be interested in working with us.

- Spencer “keeping it primitive” Tipping, Software Engineer

We're hiring
See Job Openings