空间索引 - GeoHash算法及其实现优化

h1,h2,h3,h4,h5,h6,p,blockquote { margin: 0; padding: 0 }
body { font-family: "Helvetica Neue", Helvetica, "Hiragino Sans GB", Arial, sans-serif; font-size: 13px; line-height: 18px; color: #737373; background-color: white; margin: 10px 13px 10px 13px }
table { margin: 10px 0 15px 0; border-collapse: collapse }
td,th { border: 1px solid #ddd; padding: 3px 10px }
th { padding: 5px 10px }
a { color: #0069d6 }
a:hover { color: #0050a3; text-decoration: none }
a img { border: none }
p { margin-bottom: 9px }
h1,h2,h3,h4,h5,h6 { color: #404040; line-height: 36px }
h1 { margin-bottom: 18px; font-size: 30px }
h2 { font-size: 24px }
h3 { font-size: 18px }
h4 { font-size: 16px }
h5 { font-size: 14px }
h6 { font-size: 13px }
hr { margin: 0 0 19px; border: 0; border-bottom: 1px solid #ccc }
blockquote { padding: 13px 13px 21px 15px; margin-bottom: 18px; font-family: georgia, serif; font-style: italic }
blockquote::before { content: "“"; font-size: 40px; margin-left: -10px; font-family: georgia, serif; color: #eee }
blockquote p { font-size: 14px; font-weight: 300; line-height: 18px; margin-bottom: 0; font-style: italic }
code,pre { font-family: Monaco, Andale Mono, Courier New, monospace }
code { background-color: #fee9cc; color: rgba(0, 0, 0, 0.75); padding: 1px 3px; font-size: 12px }
pre { display: block; padding: 14px; margin: 0 0 18px; line-height: 16px; font-size: 11px; border: 1px solid #d9d9d9; white-space: pre-wrap }
pre code { background-color: #fff; color: #737373; font-size: 11px; padding: 0 }
sup { font-size: 0.83em; vertical-align: super; line-height: 0 }
kbd { display: inline-block; padding: 3px 5px; font-size: 11px; line-height: 10px; color: #555; vertical-align: middle; background-color: #fcfcfc; border: solid 1px #ccc; border-bottom-color: #bbb }
* { }
code[class*="language-"],pre[class*="language-"] { color: black; background: none; font-family: Consolas, Monaco, "Andale Mono", "Ubuntu Mono", monospace; text-align: left; white-space: pre; word-spacing: normal; line-height: 1.5 }
pre[class*="language-"] { position: relative; margin: .5em 0; border-left: 10px solid #358ccb; background-color: #fdfdfd; background-image: linear-gradient(transparent 50%, rgba(69, 142, 209, 0.04) 50%); overflow: visible; padding: 0 }
code[class*="language"] { max-height: inherit; height: 100%; padding: 0 1em; display: block; overflow: auto }
:not(pre)>code[class*="language-"],pre[class*="language-"] { background-color: #fdfdfd; margin-bottom: 1em }
:not(pre)>code[class*="language-"] { position: relative; padding: .2em; color: #c92c2c; border: 1px solid rgba(0, 0, 0, 0.1); display: inline; white-space: normal }
pre[class*="language-"]::before,pre[class*="language-"]::after { content: ""; z-index: -2; display: block; position: absolute; bottom: 0.75em; left: 0.18em; width: 40%; height: 20%; max-height: 13em }
:not(pre)>code[class*="language-"]::after,pre[class*="language-"]::after { right: 0.75em; left: auto }
.token.comment,.token.block-comment,.token.prolog,.token.doctype,.token.cdata { color: #7D8B99 }
.token.punctuation { color: #5F6364 }
.token.property,.token.tag,.token.boolean,.token.number,.token.function-name,.token.constant,.token.symbol,.token.deleted { color: #c92c2c }
.token.selector,.token.attr-name,.token.string,.token.char,.token.function,.token.builtin,.token.inserted { color: #2f9c0a }
.token.operator,.token.entity,.token.url,.token.variable { color: #a67f59; background: rgba(255, 255, 255, 0.5) }
.token.atrule,.token.attr-value,.token.keyword,.token.class-name { color: #1990b8 }
.token.regex,.token.important { color: #e90 }
.language-css .token.string,.style .token.string { color: #a67f59; background: rgba(255, 255, 255, 0.5) }
.token.important { font-weight: normal }
.token.bold { font-weight: bold }
.token.italic { font-style: italic }
.token.entity { cursor: help }
.namespace { opacity: .7 }
.token.tab:not(:empty)::before,.token.cr::before,.token.lf::before { color: #e0d7d1 }
pre[class*="language-"].line-numbers { padding-left: 0 }
pre[class*="language-"].line-numbers code { padding-left: 3.8em }
pre[class*="language-"].line-numbers .line-numbers-rows { left: 0 }
pre[class*="language-"][data-line] { padding-top: 0; padding-bottom: 0; padding-left: 0 }
pre[data-line] code { position: relative; padding-left: 4em }
pre .line-highlight { margin-top: 0 }
pre.line-numbers { position: relative; padding-left: 3.8em; counter-reset: linenumber }
pre.line-numbers>code { position: relative }
.line-numbers .line-numbers-rows { position: absolute; top: 0; font-size: 100%; left: -3.8em; width: 3em; letter-spacing: -1px; border-right: 1px solid #999 }
.line-numbers-rows>span { display: block; counter-increment: linenumber }
.line-numbers-rows>span::before { content: counter(linenumber); color: #999; display: block; padding-right: 0.8em; text-align: right }

前言

上篇博客中提到了空间索引的用途和多种数据库对空间索引的支持情况,那么在应用层以下,好学的小伙伴应该会考虑空间索引的实现原理了。

目前空间索引的实现有 R树和其变种GIST树、四叉树、网格索引等。 网格索引不再多提,使用普通的hash表存储地点和风格之间的映射来实现。今天要介绍的GeoHash算法实现的空间索引,虽然是以B树实现,但我认为它也借用网格索引的一部分思想。


GeoHash

原理

GeoHash 算法的原理说起来是很简单的,如下图:

  1. 从横向上将整个方形纸分为左右两份,左侧部分为标记为 0, 右侧部分标记为 1
  2. 再将红点所在的部分划分为左右两块,再对红点位置做同样的标识,最后得出红点在横向上的标识为 10;
  3. 在纵向上对方形纸做同样的划分,左侧标识为0,右侧标识为 1,得出红点位置在纵向上的标识为 01;
  4. 将横向标识和纵向标识合并,规则为 纵向在奇数位,横向在偶数位 (也可纵横相反,但要在整个系统内保持一致),得出红点在方形纸上的标识为 1001;

只标记一个方格显得看不出什么规律,如果我们把这些都空格都标识后会发现 被划分在角落里的四个方格会有同样的前缀,如下图所示。

同样的前缀意味着可以使用 B树 索引查找有相同前缀的点作为附近的点,GeoHash 算法便是这些同样的前缀上面做文章。

墨卡托投影

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

墨卡托投影简单地说,就是可以 把整个地球平面作为一个正方形来处理,当然地球平面不是严格的正方形,此投影在两极附近的点会有误差,本文专注于原理,纠偏就不多提了(我也不懂,逃)。

实现

按照墨卡托投影的平面,我们可以按照上面划分方格纸的方式来将整个地球表面划分为各个小方格。

如(116.276349, 40.040875)这个点的经度划分:

  1. 经度在 [-180,0) 范围内的标识为0,经度范围在 [0, 180) 度的标识为 1;
  2. 继续划分,经度范围在 [0,90) 的标识为 0,经度范围在 [90,180) 的标识为 1;
  3. 这样,我们划分 20 次,方格的精度(见文末对照表)已达到 2m,得到经度的标识二进制串为11010010101011110111;
  4. 对纬度同样划分,得到纬度的标识二进制串为10111000111100100111;
  5. 我们对它组合,得到40位的二进制串11011 01110 00010 01110 11100 10111 01001 11111;
  6. 我们将这个二进制串使用 base32编码(原理同base64,可以见我的另一篇文章:WEB开发中的字符集和编码,位编码映射表见下),得到 GeoHash 编码为 3OCO4XJ7;

那么GeoHash编码前缀为 3OCO4XJ7的地理点就是离 (116.276349, 40.040875)两米内的点。如果我们把地理位置点和其GeoHash编码存入数据库的话,我们要查找 附近两米点的点,只需要限定条件 geo_code like ‘3OCO4XJ7%‘就行了;

边界点问题

可是最简版的 GeoHash 还有一个弱点,如下图:

如果每个方格的精度为 2km,那么我们直接按照前缀查询红点附近 2km 的点是查找不到离它很近的黑点的。

要解决这个问题,我们就需要所其周边八个方格也考虑上,将自身方格和周边八个方格内的点都遍历一次,再返回符合要求的点。那么如何知道周边方格的前缀呢?

仔细观察相邻方格,我们会发现两个小方格会在 经度或纬度的二进制码上相差1;我们通过 GeoHash 码反向解析出二进制码后,将其经度或纬度(或两者)的二进制码加一,再次组合为 GeoHash 码。


Redis的GEO函数

问题

我们常见的需求是查找 n米 范围内的点,那么 n米 与 GeoHash 码位数之间的映射如何实现呢?由于 GeoHash 码是由5位二进制码组成,每少一位,精度就会损失 2e(5/2)

方法当然有的,我们将二进制GeoHash码直接索引就可以,但很长的索引长度会导致 B树 索引查询效率会迅速下降。

方案

于是我们接着寻找解决方案,既然使用 base32 转换为 32进制码 会不好控制精度,保持二进制又导致索引长度过长,那么进制位数和索引长度有没有一个平衡呢?

另外 Redis 的 sorted set 支持 64位 的 double 类型的 score,我们把二进制的 GeoHash 码转为十进制放入 Redis 的 sorted set 中,不是可以实现 log(n)的查询效率了么。

说实话第一次看到 Redis 的 GEO 系列函数的时候我的内心是崩溃的,原来自己感觉极其良好的设计早已被人实现了(虽然这种情况经常出现)。。。

当然不能就这么算了,于是我使用PHP造了一遍轮子。。。

主要步骤如下:


代码实现

实现中我将 GeoHash 的最大精度设置为26位,此时它的距离精度为 0.3m。当然我们也可以充分利用 Redis 的 sorted set 的 score,设置精度为 32 位,刚好使用它的 double 类型。

数据入库:

将经纬度通过 GeoHash 算法获取到二进制 GeoHash 码,并将其转成十进制作为这个点的 score 存入 Redis 的 sorted set;

// GeoHash核心方法 传入float类型的度数和其对应的范围,经度和纬度公用方法
public function getBits($loc, $range, $level = self::LEVEL_MAX) {
    $bits = ‘‘;
    for ($i = 0; $i < $level; $i++) {
        $mid = ($range[‘min‘] + $range[‘max‘]) / 2;
        if ($loc < $mid) {
            $bits .= ‘0‘;
            $range = [‘min‘ => $range[‘min‘], ‘max‘ => $mid];
        } else {
            $bits .= ‘1‘;
            $range = [‘min‘ => $mid, ‘max‘ => $range[‘max‘]];
        }
    }

    return $bits;
}     

另外 php 的 bindec($bin_str) 方法能快速把二进制字符串转为十进制数字。

根据查询范围半径获取精度

上文说过,精度是由地图的划分次数决定的,划分次数多了,范围就小了,查询的出的数据就不全;划分次数少了,范围就会大了,我们对数据过滤时就会有过多的损耗。

private function getLevel($range_meter){
    $level = 0;
    $global = self::MERCATOR_LENGTH;
    while ($global > $range_meter) {
        $global /= 2;
        $level++;
    }

    return $level;
}   

上面代码的思想来自redis geo函数源码,真的很巧妙。

在墨卡托投影下,地球的表面可以作为一个正方形来看,它的边是地球周长中最长的一个。而学过初中地理的我们知道:“地球是一个两极稍扁,赤道略鼓的球体”,那么它最长的一个周长就是赤道周长了,于是我们得知墨卡托投影的长边为 2*PI*R=40075452.74M;

于是我们拿正方形的一个边来不停地进行二次划分,直到划分后的结果刚好比范围半径长,那么它构成的一个方块,便是我们需要的方格。

数据查询

数据查询时,我们需要获取中间方块的最小 score 值和其范围,最小 score 值很简单,直接将二进制位不足52位的在后面补0

此外,为了避免边界点问题,我们还需要把周围八个方格的 score 值范围也获取到。

我们在划分地图时,每多划分一次,会添加经度和纬度两个二进制位,在精度最高时,那么每一个方格的最大值和最小值之间差1。由此,我们通过下面的方法获取到一个方格的最大和最小 score 值之差。

private function getLevelRange($level) {
    $range = pow(2, 2 * (self::LEVEL_MAX - $level));

    return $range;
}

再由上面提过的边界点问题的解决方案,获取到周边八个方格的最小 score 值。

使用 Redis 的ZRANGEBYSCORE key hashInt hashInt+range命令将这九个方格内的点全部取到,再遍历九个方格,将距离不符合的数据过滤掉。


小结

花费了十多个小时,总算将 GeoHash 完全整体了一遍,完全理解 GeoHash 并没有想像中的那么简单。除了 GeoHash,四叉树和R树据说查询效率会更高,有时间再研究一下。

如果您觉得本文对您有帮助,可以点击下面的 推荐 支持一下我。博客一直在更新,欢迎 关注 。

参考:

GeoHash核心原理解析

Redis GEO 源码注释

GeoHash位数精度对照表(wiki百科):

GeoHash length lat bits lng bits lat error lng error km error
1 2 3 ±23 ±23 ±2500
2 5 5 ±2.8 ±5.6 ±630
3 7 8 ±0.70 ±0.70 ±78
4 10 10 ±0.087 ±0.18 ±20
5 12 13 ±0.022 ±0.022 ±2.4
6 15 15 ±0.0027 ±0.0055 ±0.61
7 17 18 ±0.00068 ±0.00068 ±0.076
8 20 20 ±0.000085 ±0.00017 ±0.019

base32 编码映射表:

Value Symbol Value Symbol Value Symbol Value Symbol
0 A 9 J 18 S 27 3
1 B 10 K 19 T 28 4
2 C 11 L 20 U 29 5
3 D 12 M 21 V 30 6
4 E 13 N 22 W 31 7
5 F 14 O 23 X    
6 G 15 P 24 Y    
7 H 16 Q 25 Z    
8 I 17 R 26 2    
时间: 05-14

空间索引 - GeoHash算法及其实现优化的相关文章

高斯模糊算法的全面优化过程分享(二)。

      相关链接: 高斯模糊算法的全面优化过程分享(一) 在高斯模糊算法的全面优化过程分享(一)一文中我们已经给出了一种相当高性能的高斯模糊过程,但是优化没有终点,经过上一个星期的发愤图强和测试,对该算法的效率提升又有了一个新的高度,这里把优化过程中的一些心得和收获用文字的形式记录下来. 第一个尝试   直接使用内联汇编替代intrinsics代码(无效) 我在某篇博客里看到说intrinsics语法虽然简化了SSE编程的难度,但是他无法直接控制XMM0-XMM7寄存器,很多指令中间都会用内

LBS地理位置距离计算方法之geohash算法

随着移动终端的普及,很多应用都基于LBS功能,附近的某某(餐馆.银行.妹纸等等).基础数据中,一般保存了目标位置的经纬度:利用用户提供的经纬度,进行对比,从而获得是否在附近.这里需要在设置出一个字段,是关于编码的字段,一会看下文哈…… 地理位置距离实现目标:查找附近多少公里内的人或者商家 比如:微信.陌陌.美团.基于O2O的一些APP这些应用或者移动网页都需要用到地理位置计算 目前来说:移动地理位置距离计算比较好的算法是geohash,特此整理分享. geohash有以下几个特点: 第一:geo

KMP串匹配算法解析与优化

朴素串匹配算法说明 串匹配算法最常用的情形是从一篇文档中查找指定文本.需要查找的文本叫做模式串,需要从中查找模式串的串暂且叫做查找串吧. 为了更好理解KMP算法,我们先这样看待一下朴素匹配算法吧.朴素串匹配算法是这样的,当模式串的某一位置失配时将失配位置的上一位置与查找串的该位置对齐再从头开始比较模式串的每一个位置.如下图所示. KMP串匹配算法解析 KMP串匹配算法是Knuth-Morris-Pratt算法的简称,KMP算法的思想就是当模式串的某一位置失配时,能不能将更前面的位置与查找串的该位

hiho一下 第四十八周 拓扑排序&#183;二【拓扑排序的应用 + 静态数组 + 拓扑排序算法的时间优化】

题目1 : 拓扑排序·二 时间限制:10000ms 单点时限:1000ms 内存限制:256MB 描述 小Hi和小Ho所在学校的校园网被黑客入侵并投放了病毒.这事在校内BBS上立刻引起了大家的讨论,当然小Hi和小Ho也参与到了其中.从大家各自了解的情况中,小Hi和小Ho整理得到了以下的信息: 校园网主干是由N个节点(编号1..N)组成,这些节点之间有一些单向的网路连接.若存在一条网路连接(u,v)链接了节点u和节点v,则节点u可以向节点v发送信息,但是节点v不能通过该链接向节点u发送信息. 在刚

蓝桥杯 算法提高 道路和航路 满分AC ,SPFA算法的SLF优化,测试数据还是比较水的,貌似没有重边

算法提高 道路和航路 时间限制:1.0s   内存限制:256.0MB 问题描述 农夫约翰正在针对一个新区域的牛奶配送合同进行研究.他打算分发牛奶到T个城镇(标号为1..T),这些城镇通过R条标号为(1..R)的道路和P条标号为(1..P)的航路相连. 每一条公路i或者航路i表示成连接城镇Ai(1<=A_i<=T)和Bi(1<=Bi<=T)代价为Ci.每一条公路,Ci的范围为0<=Ci<=10,000:由于奇怪的运营策略,每一条航路的Ci可能为负的,也就是-10,000

最短路问题(dijkstral 算法)(优化待续)

迪杰斯特拉算法是由荷兰计算机科学家狄克斯特拉于1959 年提出的,因此又叫狄克斯特拉算法.是从一个顶点到其余各顶点的最短路径算法,解决的是有向图中最短路径问题.迪杰斯特拉算法主要特点是以起始点为中心向外层层扩展,直到扩展到终点为止. 该算法复杂度为n^2 这里有一篇讲解的很清晰的文章:http://blog.chinaunix.net/uid-26548237-id-3834514.html 下面说说我个人的理解: 就以这张图为例: 要找出A点到其他点的最短路径,应该怎么找? 这个算法的思路是:

最短路问题(floyd算法)(优化待续)

问题描述: 最短路问题(short-path problem):若网络中的每条边都有一个数值(长度.成本.时间等),则找出两节点(通常是源节点和阱节点)之间总权和最小的路径就是最短路问题.最短路问题是网络理论解决的典型问题之一,可用来解决管路铺设.线路安装.厂区布局和设备更新等实际问题. 1.floyd算法 算法描述: Floyd算法又称为插点法,是一种用于寻找给定的加权图中多源点之间最短路径的算法.该算法名称以创始人之一.1978年图灵奖获得者.斯坦福大学计算机科学系教授罗伯特·弗洛伊德命名.

php通过geohash算法实现查找附近的商铺

geohash有以下几个特点: 首先,geohash用一个字符串表示经度和纬度两个坐标.利用geohash,只需在一列上应用索引即可. 其次,geohash表示的并不是一个点,而是一个矩形区域.比如编码wx4g0ec19,它表示的是一个矩形区域. 使用者可以发布地址编码,既能表明自己位于北海公园附近,又不至于暴露自己的精确坐标,有助于隐私保护. 第三,编码的前缀可以表示更大的区域.例如wx4g0ec1,它的前缀wx4g0e表示包含编码wx4g0ec1在内的更大范围. 这个特性可以用于附近地点搜索

快速排序算法实现及优化

#include  <iostream> #include  <cstdio> using namespace std; int a[100]; void quicksort(int left,int right) { int i=left,j=right,temp=a[left]; if(left>right) return; while(i!=j) { while(i<j&&a[j]>=temp) j--; while(i<j&&

《算法竞赛入门经典》之“算法设计与优化策略”

一.构造法 UVA 120 Stacks of Flapjacks Time Limit: 3000MS     64bit IO Format: %lld & %llu Submit Status uDebug Description Background Stacks and Queues are often considered the bread and butter of data structures and find use in architecture, parsing, op