聊聊GeoHash

版权声明 本站原创文章 由 萌叔 发表
转载请注明 萌叔 | https://vearne.cc

本文基于 Redis 3.2.0

前言

鲁班门前弄大斧，笔者今天要聊一下geohash。
redis/elasticsearch/mongoDB 中地理位置的搜索，比如某个位置，附近1公里内的POI点，类似这样的功能都是基于geohash算法实现的。

首先看看维基百科的定义

Geohash is a public domain geocode system invented in 2008 by Gustavo Niemeyer[1], which encodes a geographic location into a short string of letters and digits. It is a hierarchical spatial data structure which subdivides space into buckets of grid shape, which is one of the many applications of what is known as a Z-order curve, and generally space-filling curves.

翻译

Geohash是Gustavo Niemeyer于2008年发明的公共领域地理编码系统，将地理位置编码为一小段字母和数字。 它是一种分层的空间数据结构，可将空间细分为网格状的桶。

本文在描述geohash的过程会穿插介绍geohash在redis中的具体实现，以便于读者更好的理解

1. 曼托卡投影 Mercator projection

见参考资料1

墨卡托投影，是正轴等角圆柱投影。由荷兰地图学家墨卡托(G.Mercator)于1569年创立。假想一个与地轴方向一致的圆柱切或割于地球，按等角条件，将经纬网投影到圆柱面上，将圆柱面展为平面后，即得本投影。

图1

简单而言就是将3维的球面坐标系转换为2维的平面坐标。墨卡托投影导致了二极地区(南极、北极)被放大，越接近高纬度地区，形变越严重。

2. geohash的理念

整个二维平面(维度:[-90, 90],经度:[-180, 180]) 被划分成若干的小方格。

图2
图1所示的方格可以进一步切分，形成更多的方格如图2。

显然地球上的每个区域都能够被方格所表示。大的区域对应大的方格，小的区域对应小方格。为了方便传输对每个方格进行base32编码。
geohash编码的特点是

• geohash编码越长所表示的区域越小,越精确
  5位的编码能表示10平方千米范围的矩形区域，而6位编码能表示更精细的区域（约0.34平方千米）
  
  图3 精度速查表

• 相邻近的区域有公共前缀
  公共前缀越长，2个区域挨得越近

比如华东地区可以用 "wt" 表示
而上海和南京则可以表示为 "wtw3"和"wtsq"
上海徐汇区的美罗城可以表示为 "wtw37q"

图4

读者可以用参考资料8的链接，将经纬度转换成geohash编码

图5

3. geohash编码的实现

geohash编码实际是将二维平面坐标转换为了一维线性坐标。回顾一下，从三维球面坐标系，到二维平面坐标，再到一维线性坐标，不断再降维度。

(longitude, latitude) -> geohash

3.1 把经度和维度值转换成整型

redis中实现

#define GEO_STEP_MAX 26

long = ((longitude - (-180)) / (180 - (-180))) * (1 << GEO_STEP_MAX)
lat = ((latitude - (-85))/ (85- (-85))) * (1 << GEO_STEP_MAX)

以上代码将经纬度从浮点数转换成了整数，各占26bits。
提示

• 为什么是52，根据精度速查表，只要有52个bits，误差已经可以降低到0.5米以下，所以不需要在更多bits了
• 这里维度的范围不是[-90, 90], 因为高纬度地区几乎没有人居住，所以这部分区域实际是被忽略了
  实际这也是一种巧妙优化手段

比如考虑中国大概位于东经70度到140度之间，北纬0度到60度之间，我们可以用更少的bits来表征一个POI，而且精度不会下降。但是通过这种方法得到的不是标准的geoHash值，如果与其它系统交互，需要进行转换

3.2 转换为一维线性坐标

为实现将二维平面坐标转换为了一维线性坐标，geohash引入了Z阶空间填充曲线，见参考资料9

图6

把long和lat都转换成2进制，然后按照偶数位放经度，奇数位放纬度的规则重新组成bit数组
注意: 这里的偶数位是从第0位开始。 问一个问题，为什么第0位要是精度，而不是维度。
提示考虑下精度和维度的可取值范围

![

Z阶曲线在局部是有序的，在区域的边界上突变问题, 在每个 Z 字母的拐角，都有可能出现突变
希尔伯特曲线 能够解决Z阶曲线的突变问题，见参考资料9

4. 求附近的POI

图7
如上图所示，红色的点是定位点，黄色的点是它附近的POI点，求在它附近(R距离)内的POI点。定位点在区域E，区域A、B、C、D、F、G、H、I是区域E的邻近区域。

4.1 确定定位点所在区域(box)的大小

基本的要求是区域(方格)的边长大于2*R

/* You must *ONLY* estimate steps when you are encoding.
 * If you are decoding, always decode to GEO_STEP_MAX (26). */
uint8_t geohashEstimateStepsByRadius(double range_meters, double lat) {
    if (range_meters == 0) return 26;
    int step = 1;
    while (range_meters < MERCATOR_MAX) {
        range_meters *= 2;
        step++;
    }
    step -= 2; /* Make sure range is included in the worst case. */
    /* Wider range torwards the poles... Note: it is possible to do better
     * than this approximation by computing the distance between meridians
     * at this latitude, but this does the trick for now. */
    if (lat > 67 || lat < -67) step--;
    if (lat > 80 || lat < -80) step--;

    /* Frame to valid range. */
    if (step < 1) step = 1;
    if (step > 26) step = 25;
    return step;
}

4.2 确定实际需要检索的邻近区域(box)

从图7可以看出，区域A、D、G、H、I不可能会包含我们所要的POI点，因此只需要搜索区域B、C、E、F。

对区域B、C、E、F的搜索，直接利用Z阶空间曲线的局部有序性，可以直接在zset中进行检索。（geoHash作为score值），时间复杂是O(log2N)

提问: 经度180(-180度)的邻近区域怎么算？
见图二，由于180(-180度)大致穿过的白令海峡附近，经线穿过的区域几乎都是海洋, 也就是说很少可能有POI点，因此Redis 3.2.0的实现没有考虑这个问题

4.3 计算附近的点

对4.2 得到的POI再计算一次球面距离(见参考资料10)，以确保POI点都在R距离内

后记

虽然希尔伯特曲线能够解决Z阶曲线的突变问题，但是它的实现较为复杂。因此目前redis的geoHash实现仍然使用Z阶空间曲线。

2021年4月7日
elasticsearch 5.x以后，一维数值型字段，二维的geo信息以及多维数值型向量，索引都使用的是Block k-d tree，不再是GeoHash。

参考资料

1. Mercator projection
2. Geohash
3. Base32
4. 地理空间距离计算优化
5. 深入浅出空间索引：为什么需要空间索引
6. GeoHash核心原理解析
7. 高效的多维空间点索引算法 — Geohash 和 Google S2
8. geohash.org
9. Z-order curve
10. Haversine formula

---

请我喝瓶饮料