最近有项目需要使用 geohash 来生成经纬度之间的相临关系。语言选择使用正在学习的 go,库使用 mmcloughlin/geohash 。
为了更好学习 go 语言,并且稍微深入的理解 geohash 原理。花了些时间来研究 mmcloughlin 的实现,将算法和语言的实现记录在这里备忘。
经纬相当于二维的一个坐标,将位置信息表示出来。实际上,因为精度问题,经纬度指向的是一个区域,而非一个点。
经纬度表示有什么问题呢?我想到的有两个问题:
Geohash 的思路是将经纬度二维信息编织(interleave)成一个一维信息。并且与二维数据信息一样,有严格的顺序性,这样只需要按顺序就能获取点或点的周边区域。
经纬度分别是 [-180,180]
和 [-90,90]
的度数。这种保存不利于转换成一维信息。所以需要将两个维度表示成标准形式。
标准形式是转换为对应于整个轴线范围的归一化数值。
可以这样简单的理解:
代码中实现比例直接,以 python 说明
lng = (lng + 90.0) / 180.0
print lng * ( 2 ** 32 )
经度和纬度已经被编码成 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),
}
}