[990]Geohash算法原理及实现

最近需要实现一个功能，查找车辆附近的加油站，如果车和加油站距离在200米以内，则查找成功。

加油站数量肯定不小，能否缩小查找范围，否则以遍历形式，效率肯定高不了。

Geohash算法就是将经纬度编码，将二维变一维，给地址位置分区的一种算法。

经纬度常识

• 经线是纵的，经度是横的，用于表示不同的经线，纬线是横的，纬度是纵的，用于表示不同的纬线，如下图
  

• 纬线：地球仪上的横线，lat，赤道是最大的纬线，从赤道开始分为北纬和南纬，都是0-90°，纬线是角度数值，并不是米；
• 经线：地球仪上的竖线，lng，子午线为0°，分为西经和东经，都是0-180°，经线也是角度数值；
• 经纬线和米的换算：经度或者纬度0.00001度，约等于1米，这个在GPS测算距离的时候可以体会到，GPS只要精确到小数点后五位，就是10米范围内的精度
• 经度0度的位置为本初子午线，在180度的位置转为西经，数字由大到小依次经过北美洲到达西欧.纬度0度的位置为赤道

基本原理

GeoHash是一种地址编码方法。他能够把二维的空间经纬度数据编码成一个字符串

我们知道，经度范围是东经180到西经180，纬度范围是南纬90到北纬90，我们设定西经为负，南纬为负，所以地球上的经度范围就是[-180， 180]，纬度范围就是[-90，90]。如果以本初子午线、赤道为界，地球可以分成4个部分。

如果纬度范围[-90°, 0°)用二进制0代表，（0°, 90°]用二进制1代表，经度范围[-180°, 0°)用二进制0代表，（0°, 180°]用二进制1代表，那么地球可以分成如下4个部分

如果在小块范围内递归对半划分呢？

可以看到，划分的区域更多了，也更精确了。geohash算法就是基于这种思想，划分的次数更多，区域更多，区域面积更小了。通过将经纬度编码，给地理位置分区

Geohash算法

Geohash算法一共有三步。

首先将经纬度变成二进制。

比如这样一个点（39.923201, 116.390705）
纬度的范围是（-90，90），其中间值为0。对于纬度39.923201，在区间（0，90）中，因此得到一个1；（0，90）区间的中间值为45度，纬度39.923201小于45，因此得到一个0，依次计算下去，即可得到纬度的二进制表示，如下表：

最后得到纬度的二进制表示为：

10111000110001111001

同理可以得到经度116.390705的二进制表示为：

11010010110001000100

第2步，就是将经纬度合并。

经度占偶数位，纬度占奇数位，注意，0也是偶数位。

11100 11101 00100 01111 00000 01101 01011 00001

第3步，按照Base32进行编码

Base32编码表（其中一种）：

Base32编码表的其中一种如下，是用0-9、b-z（去掉a, i, l, o）这32个字母进行编码。具体操作是先将上一步得到的合并后二进制转换为10进制数据，然后对应生成Base32码。需要注意的是，将5个二进制位转换成一个base32码。上例最终得到的值为

wx4g0ec1

Geohash比直接用经纬度的高效很多，而且使用者可以发布地址编码，既能表明自己位于北海公园附近，又不至于暴露自己的精确坐标，有助于隐私保护。

• GeoHash用一个字符串表示经度和纬度两个坐标。在数据库中可以实现在一列上应用索引（某些情况下无法在两列上同时应用索引）
• GeoHash表示的并不是一个点，而是一个矩形区域
• GeoHash编码的前缀可以表示更大的区域。例如wx4g0ec1，它的前缀wx4g0e表示包含编码wx4g0ec1在内的更大范围。 这个特性可以用于附近地点搜索

GeoHash的精度

以上表格来自:https://en.wikipedia.org/wiki/Geohash
从表格中可以看出，如果GeoHash值取7位长度，误差大概在76米左右；如果取8位，误差为19米左右。

编码越长，表示的范围越小，位置也越精确。因此我们就可以通过比较GeoHash匹配的位数来判断两个点之间的大概距离。

具体的计算方法

Latitude的范围是：-90 到 +90
Longitude的范围：-180 到 +180
地球参考球体的周长：40075016.68米

geohash长度	Lat位数	Lng位数	Lat误差	Lng误差	km误差
1	2	3	±23	±23	±2500
2	5	5	± 2.8	±5.6	±630
3	7	8	± 0.70	± 0.7	±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.000086	±0.000172	±0.01911
9	22	23	±0.000021	±0.000021	±0.00478
10	25	25	±0.00000268	±0.00000536	±0.0005971
11	27	28	±0.00000067	±0.00000067	±0.0001492
12	30	30	±0.00000008	±0.00000017	±0.0000186

1.纬度相同，经度不同

在纬度相同的情况下：
经度每隔0.00001度，距离相差约1米；
每隔0.0001度，距离相差约10米；
每隔0.001度，距离相差约100米；
每隔0.01度，距离相差约1000米；
每隔0.1度，距离相差约10000米。

2.经度相同，纬度不同

纬度每隔0.00001度，距离相差约1.1米；
每隔0.0001度，距离相差约11米；
每隔0.001度，距离相差约111米；
每隔0.01度，距离相差约1113米；
每隔0.1度，距离相差约11132米。

字符串越长代表的精度越高 5位的编码能表示10平方千米，而6位的编码约0.34平方千米

问题

geohash算法有两个问题。首先是边缘问题。

如图，如果车在红点位置，区域内还有一个黄点。相邻区域内的绿点明显离红点更近。但因为黄点的编码和红点一样，最终找到的将是黄点。这就有问题了。

要解决这个问题，很简单，只要再查找周边8个区域内的点，看哪个离自己更近即可。

另外就是曲线突变问题。

本文第2张图片比较好地解释了这个问题。其中0111和1000两个编码非常相近，但它们的实际距离确很远。所以编码相近的两个单位，并不一定真实距离很近，这需要实际计算两个点的距离才行。

代码实现

geohash原理清楚后，代码实现就比较简单了。不过仍然有一个问题需要解决，就是如何计算周边的8个区域key值呢

假设我们计算的key值是6位，那么二进制位数就是 6*5 = 30位，所以经纬度分别是15位。我们以纬度为例，纬度会均分15次。这样我们很容易能够算出15次后，划分的最小单位是多少

  private void setMinLatLng() {
    minLat = MAXLAT - MINLAT;
    for (int i = 0; i < numbits; i++) {
        minLat /= 2.0;
    }
    minLng = MAXLNG - MINLNG;
    for (int i = 0; i < numbits; i++) {
        minLng /= 2.0;
    }
}

得到了最小单位，那么周边区域的经纬度也可以计算得到了。比如说左边区域的经度肯定是自身经度减去最小经度单位。纬度也可以通过加减，得到上下的纬度值，最终周围8个单位也可以计算得到。

可以到 http://geohash.co/ 进行geohash编码，以确定自己代码是否写错

整体代码如下所示：

public class GeoHash {
public static final double MINLAT = -90;
public static final double MAXLAT = 90;
public static final double MINLNG = -180;
public static final double MAXLNG = 180;

private static int numbits = 3 * 5; //经纬度单独编码长度

private static double minLat;
private static double minLng;

private 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' };

//定义编码映射关系
final static HashMap<Character, Integer> lookup = new HashMap<Character, Integer>();
//初始化编码映射内容
static {
    int i = 0;
    for (char c : digits)
        lookup.put(c, i++);
}

public GeoHash(){
    setMinLatLng();
}

public 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');
    }
    String code = base32(Long.parseLong(buffer.toString(), 2));
    //Log.i("okunu", "encode  lat = " + lat + "  lng = " + lon + "  code = " + code);
    return code;
}

public ArrayList<String> getArroundGeoHash(double lat, double lon){
    //Log.i("okunu", "getArroundGeoHash  lat = " + lat + "  lng = " + lon);
    ArrayList<String> list = new ArrayList<>();
    double uplat = lat + minLat;
    double downLat = lat - minLat;

    double leftlng = lon - minLng;
    double rightLng = lon + minLng;

    String leftUp = encode(uplat, leftlng);
    list.add(leftUp);

    String leftMid = encode(lat, leftlng);
    list.add(leftMid);

    String leftDown = encode(downLat, leftlng);
    list.add(leftDown);

    String midUp = encode(uplat, lon);
    list.add(midUp);

    String midMid = encode(lat, lon);
    list.add(midMid);

    String midDown = encode(downLat, lon);
    list.add(midDown);

    String rightUp = encode(uplat, rightLng);
    list.add(rightUp);

    String rightMid = encode(lat, rightLng);
    list.add(rightMid);

    String rightDown = encode(downLat, rightLng);
    list.add(rightDown);

    //Log.i("okunu", "getArroundGeoHash list = " + list.toString());
    return list;
}

//根据经纬度和范围，获取对应的二进制
private BitSet getBits(double lat, double floor, double ceiling) {
    BitSet buffer = new BitSet(numbits);
    for (int i = 0; i < numbits; i++) {
        double mid = (floor + ceiling) / 2;
        if (lat >= mid) {
            buffer.set(i);
            floor = mid;
        } else {
            ceiling = mid;
        }
    }
    return buffer;
}

//将经纬度合并后的二进制进行指定的32位编码
private 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));
}

private void setMinLatLng() {
    minLat = MAXLAT - MINLAT;
    for (int i = 0; i < numbits; i++) {
        minLat /= 2.0;
    }
    minLng = MAXLNG - MINLNG;
    for (int i = 0; i < numbits; i++) {
        minLng /= 2.0;
    }
}

//根据二进制和范围解码
private 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;
}

//对编码后的字符串解码
public 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();

    //偶数位，经度
    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);
    }

    //奇数位，纬度
    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);
    }

    double lon = decode(lonset, -180, 180);
    double lat = decode(latset, -90, 90);

    return new double[] {lat, lon};
}

public static void main(String[] args)  throws Exception{
    GeoHash geohash = new GeoHash();
//        String s = geohash.encode(40.222012, 116.248283);
//        System.out.println(s);
    geohash.getArroundGeoHash(40.222012, 116.248283);
//        double[] geo = geohash.decode(s);
//        System.out.println(geo[0]+" "+geo[1]);
}
}

geohash在mysql中的使用

mysql中有四个内置函数跟geohash相关。

函数	说明
ST_GeoHash()	产生geohash值
ST_LongFromGeoHash()	从geohash值返回经度
ST_LatFromGeoHash()	从geohash值返回纬度
ST_PointFromGeoHash()	将geohash值转换为POINT值

具体使用方法可参考：https://dev.mysql.com/doc/refman/8.0/en/spatial-geohash-functions.html#function_st-geohash

简单使用

#创建 point 表
create table `geom`(
`id` int not null auto_increment,
`geom` point not null,
primary key(`id`)
);
 
#插入数据
insert into `geom` (`geom`) VALUES
(st_geomFromText('Point(1 1)')),
(st_geomFromText('Point(2 2)'));
 
# geohash
SELECT st_geohash(`geom`,12) as geohash from `geom`;

结果

其中st_geomFromText()可以传入WKT类型的字符串生成geometry对象。

• st_geohash(point类型,geohash长度)
• st_geohash(精度,纬度,geohash长度)

sql使用

SELECT ST_AsText(ST_PointFromGeoHash('ws7gm',0));

select ST_PointFromGeoHash('ws7gm',6);

select ST_LatFromGeoHash('ws7gm'),ST_LongFromGeoHash('ws7gm');

select st_geohash(118,24.5,5);

select distinct ST_AsText(ST_PointFromGeoHash(a.geohash5,0)) from (
SELECT distinct ST_GeoHash(longitude,latitude,5) as geohash5 from poi_miniapp_total ) a

select distinct a.geohash5,ST_LatFromGeoHash(a.geohash5) as latitude,ST_LongFromGeoHash(a.geohash5) as longitude from (
SELECT distinct ST_GeoHash(longitude,latitude,5) as geohash5 from poi_miniapp_total ) a

参考：https://www.jianshu.com/p/2fd0cf12e5ba
https://www.cnblogs.com/feiquan/p/11380461.html
https://blog.csdn.net/ununie/article/details/100521397
https://blog.csdn.net/youhongaa/article/details/78816700

精度问题：
https://blog.csdn.net/ununie/article/details/96963642
https://segmentfault.com/a/1190000002513514
https://blog.csdn.net/bitcarmanlee/article/details/55824141
https://blog.csdn.net/qq_39341048/article/details/104910984