Fork me on GitHub

版权声明 本站原创文章 由 萌叔 发表
转载请注明 萌叔 | 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年创立。假想一个与地轴方向一致的圆柱切或割于地球,按等角条件,将经纬网投影到圆柱面上,将圆柱面展为平面后,即得本投影。

image_1dpf31rquc1819fd9ol1dlb1457e.png-433.8kB
图1

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

2. geohash的理念

整个二维平面(维度:[-90, 90],经度:[-180, 180]) 被划分成若干的小方格。
W020120920381438693315.jpg-42kB
图2
图1所示的方格可以进一步切分,形成更多的方格如图2。

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

  • geohash编码越长所表示的区域越小,越精确
    5位的编码能表示10平方千米范围的矩形区域,而6位编码能表示更精细的区域(约0.34平方千米)
    屏幕快照 2019-11-12 下午2.48.54.png-226.3kB
    图3 精度速查表

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

比如华东地区可以用 "wt" 表示
而上海和南京则可以表示为 "wtw3"和"wtsq"
上海徐汇区的美罗城可以表示为 "wtw37q"
56_5.png-251.5kB
图4

读者可以用参考资料8的链接,将经纬度转换成geohash编码
屏幕快照 2019-11-12 下午1.53.53.png-658.7kB
图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

600px-Four-level_Z.svg.png-27.8kB
图6

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

![image_1dpfbcfev5ut1ntt7hkrv5939.png-29kB

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

4. 求附近的POI

image_1dpfepi5s1dgivep16qk7p51al11a.png-78.7kB
图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

请我喝瓶饮料

微信支付码

3 对 “聊聊GeoHash”的想法;

    1. 确实写得有问题

          double lat_offset =
              (latitude - lat_range->min) / (lat_range->max - lat_range->min);
          double long_offset =
              (longitude - long_range->min) / (long_range->max - long_range->min);
      
          /* convert to fixed point based on the step size */
          lat_offset *= (1 << step);
          long_offset *= (1 << step);
          hash->bits = interleave64(lat_offset, long_offset);
      

发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据