March 23, 2017

Geohash 算法与实现备忘

说明

最近有项目需要使用 geohash 来生成经纬度之间的相临关系。语言选择使用正在学习的 go,库使用 mmcloughlin/geohash

为了更好学习 go 语言,并且稍微深入的理解 geohash 原理。花了些时间来研究 mmcloughlin 的实现,将算法和语言的实现记录在这里备忘。

算法

基本思想

经纬相当于二维的一个坐标,将位置信息表示出来。实际上,因为精度问题,经纬度指向的是一个区域,而非一个点。

经纬度表示有什么问题呢?我想到的有两个问题:

  1. 使用二维坐标来表示位置,在查询一个区域时比较复杂,需要进行两次过滤
  2. 找出中心点周边的 N 个邻近点,这类场景二维坐标处理困难

Geohash 的思路是将经纬度二维信息编织(interleave)成一个一维信息。并且与二维数据信息一样,有严格的顺序性,这样只需要按顺序就能获取点或点的周边区域。

经纬度编码

经纬度分别是 [-180,180][-90,90] 的度数。这种保存不利于转换成一维信息。所以需要将两个维度表示成标准形式。

标准形式是转换为对应于整个轴线范围的归一化数值。

可以这样简单的理解:

  1. 度数的范围形成一个一米长的尺子
  2. 点的位置对应在尺子上的一个刻度
  3. 计算出刻度对应在尺子上的相对比例
  4. 比例以二进制来表示,位数越多,则越精确说明点的位置
  5. 无论精度多高,总会产生一定误差,所以实际上描述的是一个中心点加上一个误差范围

代码中实现比例直接,以 python 说明

lng = (lng + 90.0) / 180.0
print lng * ( 2 ** 32 )

interleave

经度和纬度已经被编码成 32 位的 bit,在忽略地球是个扁球形,两个维度上的每个字节所表示的长度是相近的。将经度和纬度的字节顺序的交替的编织起来,生成一个 64 位的整数。

整数在显示或输入上不方便,编码成可见字符比较简单。geohash 使用 32 个可见字符(因为是 \(2^5\) ,容易处理 ),来进行编码。

每个字符表示 5 个 bit,12 个字符可以表示 60 个 bit。经度和纬度的位数都是 30 个字。

这里代码实现交织的方式值得参考,这里引用如下:

// Spread out the 32 bits of x into 64 bits, where the bits of x occupy even
// bit positions.
func spread(x uint32) uint64 {
	X := uint64(x)
	X = (X | (X << 16)) & 0x0000ffff0000ffff
	X = (X | (X << 8)) & 0x00ff00ff00ff00ff
	X = (X | (X << 4)) & 0x0f0f0f0f0f0f0f0f
	X = (X | (X << 2)) & 0x3333333333333333
	X = (X | (X << 1)) & 0x5555555555555555
	return X
}

// Interleave the bits of x and y. In the result, x and y occupy even and odd
// bitlevels, respectively.
func interleave(x, y uint32) uint64 {
	return spread(x) | (spread(y) << 1)
}

精度

以 5 个字符来计算,能表达 25 bits。经度在前,纬度在后,因此以 13bit 表示经度,12bit 表示纬度。

经度的误差为 \({0.5}^{14} * 360 = 0.02197265625\) ,同理纬度的误差为 \({0.5}^{13} * 180 = 0.02197265625\)

以地球周长约 40000km 约数计算(即每度距离 111km),公里误差是

\[\sqrt{(0.022*0.022) + (0.022*0.022)} * 40000 / 360 = 3.45km\]

相邻

由于 geohash 有一定的精度,实际上表示的是一个方形区域。但正方形有个问题,是上下左右边距离近于对角距离。因此,为了提取所有相邻区域,需要将东南西北及东南、东北、西南、西北八个相邻区域全部考虑进去

mmcloughlin/geohash 中的代码引用如下

// Neighbors returns a slice of geohash strings that correspond to the provided
// geohash's neighbors.
func Neighbors(hash string) []string {
	box := BoundingBox(hash)
	lat, lng := box.Center()
	latDelta := box.MaxLat - box.MinLat
	lngDelta := box.MaxLng - box.MinLng
	precision := uint(len(hash))
	return []string{
		// N
		EncodeWithPrecision(lat+latDelta, lng, precision),
		// NE,
		EncodeWithPrecision(lat+latDelta, lng+lngDelta, precision),
		// E,
		EncodeWithPrecision(lat, lng+lngDelta, precision),
		// SE,
		EncodeWithPrecision(lat-latDelta, lng+lngDelta, precision),
		// S,
		EncodeWithPrecision(lat-latDelta, lng, precision),
		// SW,
		EncodeWithPrecision(lat-latDelta, lng-lngDelta, precision),
		// W,
		EncodeWithPrecision(lat, lng-lngDelta, precision),
		// NW
		EncodeWithPrecision(lat+latDelta, lng-lngDelta, precision),
	}
}