geohash算法

前言

用户附近位置如何计算?

常见的很多app功能都有附近位置的计算,那是如何实现的?是通过两点经纬度距离计算?想快速了解附近100米的超市?遍历每个超市拿经纬度计算距离显然不合适。

经纬度与物理距离介绍
经纬度是经度与纬度的合称组成一个坐标系统,称为地理坐标系统,它是一种利用三度空间的球面来定义地球上的空间的球面坐标系统,能够标示地球上的任何一个位置。

geohash算法.assets/57eb711c1d9cded841228b29b821eeba.png

在一定误差范围内,通常情况下,经纬线和米的换算为:经度或者纬度0.00001度,约等于1米。以下表格列出更细致的换算关系:

在纬度相等的情况下 在经度相等的情况下
经度每隔0.00001度,距离相差约1米;每隔0.0001度,距离相差约10米;每隔0.001度,距离相差约100米;每隔0.01度,距离相差约1000米;每隔0.1度,距离相差约10000米。 纬度每隔0.00001度,距离相差约1.1米;每隔0.0001度,距离相差约11米;每隔0.001度,距离相差约111米;每隔0.01度,距离相差约1113米;每隔0.1度,距离相差约11132米。

Geohash算法

GeoHash是空间索引的一种方式,其基本原理是将地球理解为一个二维平面,通过把二维的空间经纬度数据编码为一个字符串,可以把平面递归分解成更小的子块,每个子块在一定经纬度范围内拥有相同的编码。

以GeoHash方式建立空间索引,可以提高对空间poi数据进行经纬度检索的效率。

编码规则为:先将纬度范围(-90, 90)平分成两个区间(-90, 0)和(0, 90),如果目标维度位于前一个区间,则编码为0,否则编码为1,然后根据目标纬度所落的区间再平均分成两个区间进行编码,以此类推,直到精度满足要求,经度也用同样的算法,对(-180, 180)依次细分,然后合并经度和纬度的编码,奇数位放纬度,偶数位放经度,组成一串新的二进制编码,按照Base32进行编码。

以当前所在办公区【两江国际】的位置坐标为例, 经纬度为(104.059684,30.559545)

第一步:将经纬度转换为二进制

序号 纬度范围 划分区间0 划分区间1 30.559545所属区间
1 (-90, 90) (-90, 0.0) (0.0, 90) 1
2 (0.0, 90) (0.0, 45.0) (45.0, 90) 0
3 (0.0, 45.0) (0.0, 22.5) (22.5, 45.0) 1
4 (22.5, 45.0) (22.5, 33.75) (33.75, 45.0) 0
5 (22.5, 33.75) (22.5, 28.125) (28.125, 33.75) 1
6 (28.125, 33.75) (28.125, 30.9375) (30.9375, 33.75) 0
7 (28.125, 30.9375) (28.125, 29.53125) (29.53125, 30.9375) 1
8 (29.53125, 30.9375) (29.53125, 30.234375) (30.234375, 30.9375) 1
9 (30.234375, 30.9375) (30.234375, 30.5859375) (30.5859375, 30.9375) 0
10 (30.234375, 30.5859375) (30.234375, 30.41015625) (30.41015625, 30.5859375) 1
11 (30.41015625, 30.5859375) (30.41015625, 30.498046875) (30.498046875, 30.5859375) 1
12 (30.498046875, 30.5859375) (30.498046875, 30.541992188) (30.541992188, 30.5859375) 1
13 (30.541992188, 30.5859375) (30.541992188, 30.563964844) (30.563964844, 30.5859375) 0
14 (30.541992188, 30.563964844) (30.541992188, 30.552978516) (30.552978516, 30.563964844) 1
15 (30.552978516, 30.563964844) (30.552978516, 30.55847168) (30.55847168, 30.563964844) 1

最后得到维度的二进制编码为:101010110111011, 用同样的方式可以得到精度(104.059684)的二进制编码:110010011111111。

第二步:将经纬度的二进制编码合并

从偶数0开始,经度占偶数位,纬度占奇数位。

序号 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
编码 1 1 1 0 0 1 0 0 1 1 0 0 0 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1

第三步:将合并后的二进制数做Base32编码

按照每5位一组,分成6组,每组计算其对应的十进制数值,按照Base32进行编码。

Base32编码表的其中一种如下,是用0-9、b-z(去掉a, i, l, o)这32个字母进行编码.

geohash算法.assets/image-20230315022500048.png
1
2
11100 10011 00011 11011 11111 01111
28(w) 19(m) 3(3) 27(v) 31(z) 15(g)

最终得到的经纬度编码为:wm3vzg

如上文二进制编码的计算过程,如果递归的次数越大,则生成的二进制编码越长,因此生成的geohash编码越长,位置越精确。目前Geohash使用的精度说明如下:

geohash算法.assets/8eb57262b9cfb26dc94646465d72533c.png

GeoHash用一个字符串表示经度和纬度两个坐标, 比直接用经纬度的高效很多,而且使用者可以发布地址编码,既能表明自己位于某位置附近,又不至于暴露自己的精确坐标,有助于隐私保护。

编码过程中,通过二分范围匹配的方式来决定某个经纬坐标是编码为1还是0,因此某些邻近坐标的编码是相同的,因此GeoHash表示的并不是一个点,而是一个矩形区域。 GeoHash编码的前缀可以表示更大的区域。例如wm3vzg,它的前缀wm3vz表示包含编码wm3vzg在内的更大范围。 这个特性可以用于附近地点搜索。

网格与实际地理的精度

如果把某个区域或整个地图上的地理位置都按照Geohash编码,则会得到一个网格,编码递归粒度越细,网格的矩形区域越小,geohash编码的长度越大,则Geohash编码越精确。 不同的编码长度,生成的网格与实际地理的精度如下(Geohash字符串编码长度对应网格大小):

字符串长度 网格宽度 网格高度
1 5000Km 5000Km
2 1250Km 625Km
3 156Km 156Km
4 39.1Km 19.5Km
5 4.89Km 4.89Km
6 1.22Km 0.61Km
7 153m 153m
8 38.2m 19.1m
9 4.77m 4.77m
10 1.19m 0.596m

Geohash编码与网格

当前选取的编码长度为6,因此一个网格实际的地理差异在1.2公里与0.6公里,示例中两江国际对应的网格大致效果如图:

geohash算法.assets/b6c5b5530af147249e9a870f076f9978-1678818658611.png

使用geohash算法推算附近位置的弊端(边界问题)

geohash算法.assets/image-20230624011131681.png

邻近(8个)网格位置推算

结论
根据Geohash的编码规则将经纬度分解到二进制,结合地理常识,中心网格在南北(上下)方向上体现为纬度的变化,往北则维度的二进制加1,往南则维度的二进制减1,在东西(左右)方向上体现为经度的变化,往东则经度的二进制加1,往西则减1,可以计算出上下左右四个网格经纬度的二进制编码,再将加减得出的经纬度两两组合,计算出左上、左下、右上和右下四个网格的经纬度二进制编码,从而就可以根据Geohash的编码规则计算出周围八个网格的字符串。

正向推导

以Geohash编码长度为6为基础,网格的宽高与实际距离换算为:1.2Km*0.6Km.

参考上文提到的,在经度相同情况下,每隔0.001度,距离相差约111米。0.6Km换算为纬度为:0.005405405。

当前两江国际粗粒度的wgs84坐标(104.05503,30.562251), 纬度二进制编码:101010110111011,经度二进制编码:110010011111111, Geohash值为:wm3vzg

正北方向近邻的网格维度为增加一个网格的高度,即纬度增加0.005405405,为:30.562251 + 0.005405405 = 30.567656405, 转换为二进制编码后为(可用工具快速转换):101010110111100

正好是原纬度的二进制编码101010110111011 加1后的结果(101010110111011 + 000000000000001 = 101010110111100)

反向推导

当前两江国际粗粒度的wgs84坐标(104.05503,30.562251), 纬度二进制编码:101010110111011,经度二进制编码:110010011111111, Geohash值为:wm3vzg

基于当前坐标的网格,正北方向近邻的网格N,其纬度二进制加1后为:101010110111100,经度不变,其Geohash值为:wm3vzu

通过http://geohash.co/ 反向转换其经纬坐标为:(104.0570068359375,30.56671142578125)

通过https://www.box3.cn/tools/lbs.html 查询2个坐标的实际位置,误差在531m(符合精度范围)

geohash算法.assets/d34fee6fe2a510de033a1c103a952c4d.png

邻近8个网格位置计算

Geohash编码:wm3vzs纬度二进制编码:101010110111100经度二进制编码:110010011111110公式:(Lat_bin + 1, Lon_bin - 1) Geohash编码:wm3vzu纬度二进制编码:101010110111100经度二进制编码:110100101101010公式:(Lat_bin + 1, Lon_bin) Geohash编码:wm6jbh纬度二进制编码:101010110111100经度二进制编码:110010100000000公式:(Lat_bin + 1, Lon_bin + 1)
Geohash编码:wm3vze纬度二进制编码:101010110111011经度二进制编码:110010011111110公式:(Lat_bin, Lon_bin - 1) Geohash编码:wm3vzg纬度二进制编码:101010110111011经度二进制编码:110010011111111公式:(Lat_bin, Lon_bin) Geohash编码:wm6jb5纬度二进制编码:101010110111011经度二进制编码:110010100000000公式:(Lat_bin, Lon_bin + 1)
Geohash编码:wm3vzd纬度二进制编码:101010110111010经度二进制编码:110010011111110公式:(Lat_bin - 1, Lon_bin - 1) Geohash编码:wm3vzf纬度二进制编码:101010110111010经度二进制编码:110010011111111公式:(Lat_bin - 1, Lon_bin) Geohash编码:wm6jb4纬度二进制编码:101010110111010经度二进制编码:110010100000000公式:(Lat_bin - 1, Lon_bin + 1)

附近3公里网格模型

geohash算法.assets/dc83cf2c5b09437938234726ca167a3e.png geohash算法.assets/image-20230315031553823.png

geohash进阶原理(Z阶曲线)

上文讲了GeoHash的计算步骤,仅仅说明是什么而没有说明为什么?为什么分别给经度和维度编码?为什么需要将经纬度两串编码交叉组合成一串编码?本节试图回答这一问题。

如下图所示,我们将二进制编码的结果填写到空间中,当将空间划分为四块时候,编码的顺序分别是左下角00,左上角01,右下脚10,右上角11,也就是类似于Z的曲线,当我们递归的将各个块分解成更小的子块时,编码的顺序是自相似的(分形),每一个子快也形成Z曲线,这种类型的曲线被称为 Peano 空间填充曲线。

这种类型的空间填充曲线的优点是将二维空间转换成一维曲线(事实上是分形维),对大部分而言,编码相似的距离也相近, 但Peano空间填充曲线最大的缺点就是突变性,有些编码相邻但距离却相差很远,比如0111与1000,编码是相邻的,但距离相差很大。  

geohash算法.assets/6f061d950a7b02080c4a5e7fd6dd51d4562cc8bd.png

除 Peano 空间填充曲线外,还有很多空间填充曲线,如图所示,其中效果公认较好是 Hilbert 空间填充曲线,相较于 Peano 曲线而言,Hilbert 曲线没有较大的突变。为什么 GeoHash 不选择 Hilbert 空间填充曲线呢?可能是 Peano曲线思路以及计算上比较简单吧,事实上,Peano 曲线就是一种四叉树线性编码方式。

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
import java.util.BitSet;
import java.util.HashMap;

/**
* @Description TODO geohash算法实现工具类
*/
public class GeoHashUtil {
// todo 经纬度编码长度
private static int numbits = 6 * 5;
// todo 数字字母编码数组
final static char[] digits = {'0', '1', '2', '3', '4', '5', '6', '7',
'8', '9', 'b', 'c', 'd', 'e', 'f', 'g',
'h', 'j', 'k', 'm', 'n', 'p', 'q', 'r',
's', 't', 'u', 'v', 'w', 'x', 'y', 'z'};
// todo 存储编码字符循环存储的HashMap对象
final static HashMap<Character, Integer> lookup = new HashMap<>();

// todo 静态代码块,执行设置HashMap的key,value
static {
int i = 0;
for (char c : digits)
lookup.put(c, i++);
}

/**
* @desc todo 根据编码后的geohash字符串值,进行解码,得到经纬度数组
* @param geoHash
* @return [lat, lon]
*/
public static double[] decode(String geoHash) {
StringBuilder buffer = new StringBuilder();
for (char c : geoHash.toCharArray()) {
int i = lookup.get(c) + 32;
buffer.append(Integer.toString(i, 2).substring(1));
}

BitSet lonset = new BitSet();
BitSet latset = new BitSet();

// todo 经度,偶数位
int j = 0;
for (int i = 0; i < numbits * 2; i += 2) {
boolean isSet = false;
if (i < buffer.length())
isSet = buffer.charAt(i) == '1';
lonset.set(j++, isSet);
}

// todo 纬度,奇数位
j = 0;
for (int i = 1; i < numbits * 2; i += 2) {
boolean isSet = false;
if (i < buffer.length())
isSet = buffer.charAt(i) == '1';
latset.set(j++, isSet);
}

// todo 根据位编码、经度最小值、经度最大值计算出经度
double lon = decode(lonset, -180, 180);
// todo 根据位编码、纬度最小值、纬度最大值计算出经度
double lat = decode(latset, -90, 90);
// todo 返回纬度、经度数组
return new double[]{lat, lon};
}

/**
* @desc todo 编码方法,根据编码、维度[-90,90],经度[-180,180]
* @param bs
* @param floor
* @param ceiling
* @return
*/
private static double decode(BitSet bs, double floor, double ceiling) {
double mid = 0;
for (int i = 0; i < bs.length(); i++) {
mid = (floor + ceiling) / 2;
if (bs.get(i))
floor = mid;
else
ceiling = mid;
}
return mid;
}

/**
* @desc todo 解码方法,根据纬度、经度,返回32编码字符串
* @param lat 纬度
* @param lon 经度
* @return base32的字符串
*/
public static String encode(double lat, double lon) {
BitSet latbits = getBits(lat, -90, 90);
BitSet lonbits = getBits(lon, -180, 180);
StringBuilder buffer = new StringBuilder();
for (int i = 0; i < numbits; i++) {
buffer.append(lonbits.get(i) ? '1' : '0');
buffer.append(latbits.get(i) ? '1' : '0');
}
return base32(Long.parseLong(buffer.toString(), 2));
}

/**
* @desc todo 根据经纬度和范围,获取对应的二进制值
* @param d 经度 | 纬度
* @param floor 最小值
* @param ceiling 最大值
* @return 返回BitSet,java.util工具类,用于位移操作工具类
*/
private static BitSet getBits(double d, double floor, double ceiling) {
BitSet buffer = new BitSet(numbits);
for (int i = 0; i < numbits; i++) {
double mid = (floor + ceiling) / 2;
if (d >= mid) {
buffer.set(i);
floor = mid;
} else {
ceiling = mid;
}
}
return buffer;
}

/**
* @desc todo 将经纬度合并后二二进制进行指定32位编码
* @param i 被32编码的long值
* @return 32编码字符串
*/
private static String base32(long i) {
char[] buf = new char[65];
int charPos = 64;
boolean negative = (i < 0);
if (!negative)
i = -i;
while (i <= -32) {
buf[charPos--] = digits[(int) (-(i % 32))];
i /= 32;
}
buf[charPos] = digits[(int) (-i)];

if (negative)
buf[--charPos] = '-';
return new String(buf, charPos, (65 - charPos));
}
}