GeoHash和Google S2

生活场景

日常生活中，我们常常会有需求查找附近的xxx，比如出租车，停车场等等，实现方式多种多样，常见的有

MongoDB
PostgreSQL+PostGIS
ElasticSearch
Mysql外接正方形
Mysql+GeoHash
Redis+GeoHash

基本原理

GeoHash是一种地址编码方法，它能把一个二维的经纬度编码成一个字符串。

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

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

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

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

算法过程

以经纬度：(119.984686, 30.229212)为进行算法说明。

对纬度分区编码

对纬度30.229212进行逼近编码 (地球纬度区间是[-90,90]）：

区间[-90,90]进行二分为[-90,0),[0,90]，称为左右区间，可以确定30.229212属于右区间[0,90]，给标记为1
接着将区间[0,90]进行二分为 [0,45),[45,90]，可以确定30.229212属于左区间 [0,45)，给标记为0
递归上述过程30.229212总是属于某个区间[a,b]。随着每次迭代区间[a,b]总在缩小，并越来越逼近30.229212
如果给定的纬度（30.229212）属于左区间，则记录0，如果属于右区间则记录1，序列的长度跟给定的区间划分次数有关，如下图:

纬度范围

划分区间0(左)

划分区间1(右)

30.229212

(-90, 90)

(-90, 0)

(0, 90)

(0, 45)

(45, 90)

(0, 45)

(0, 22.5)

(22.5, 45)

(22.5, 33.75)

(33.75, 45)

(22.5, 33.75)

(22.5, 28.125)

(28.125, 33.75)

(28.125, 30.9375)

(30.9375, 33.75)

(28.125, 30.9375)

(28.125, 29.53125)

(29.53125, 30.9375)

(29.53125, 30.234375)

(30.234375, 30.9375)

(29.53125, 30.234375)

(29.53125, 29.8828125)

(29.8828125, 30.234375)

(29.8828125, 30.05859375)

(30.05859375, 30.234375)

(30.05859375, 30.146484375)

(30.146484375, 30.234375)

(30.146484375, 30.1904296875)

(30.1904296875, 30.234375)

(30.1904296875, 30.2124023438)

(30.2124023438, 30.234375)

(30.2124023438, 30.2233886719)

(30.2233886719, 30.234375)

(30.2233886719, 30.228881836)

(30.228881836, 30.234375)

对经度分区编码

同理，地球经度区间是[-180,180]，可以对经度119.984686进行编码。通过上述计算

经度范围

划分区间0(左)

划分区间1(右)

119.984686

(-180, 180)

(-180, 0)

(0, 180)

(0, 90)

(90, 180)

(90, 135)

(135, 180)

(90, 135)

(90, 112.5)

(112.5, 135)

(112.5, 123.75)

(123.75, 135)

(112.5, 123.75)

(112.5, 118.125)

(118.125, 123.75)

(118.125, 120.9375)

(120.9375, 123.75)

(118.125, 120.9375)

(118.125, 119.53125)

(119.53125, 120.9375)

(119.53125, 120.234375)

(120.234375, 120.9375)

(119.53125, 120.234375)

(119.53125, 119.8828125)

(119.8828125, 120.234375)

(119.8828125, 120.05859375)

(120.05859375, 120.234375)

(119.8828125, 120.05859375)

(119.8828125, 119.970703125)

(119.970703125, 120.05859375)

(119.970703125, 120.0146484375)

(120.0146484375, 120.05859375)

(119.970703125, 120.0146484375)

(119.970703125, 119.9926757812)

(119.9926757812, 120.0146484375)

(119.970703125, 119.9926757812)

(119.970703125, 119.9816894531)

(119.9816894531, 119.9926757812)

通过上面的计算，纬度产生的编码是 101010101111111，经度产生的编码是110101010101001

合并经纬度编码计算

偶数位放经度，奇数位放纬度，下标是从0开始计算的。

经度 1 1 0 1 0 1 0 1 0 1 0 1 0 0 1
纬度  1 0 1 0 1 0 1 0 1 1 1 1 1 1 1
—————————————————————————————————————
结果 111001100110011001110111010111

将11100 11001 10011 00111 01110 10111转成十进制，按照每5位转换成十进制，分别是28, 25, 19, 7, 14, 23。最后用0-9、b-z（去掉a, i, l, o）这32个字母进行base32编码为：wtm7fr。

地图上标记的范围：

SELECT * FROM customer WHERE geohash LIKE 'wtm7fr%'

GeoHash缺陷

曲线突变

大部分而言，编码相似的距离也相近，但是部分编码存在突变性，编码相邻但距离却相差很远，比如本文图2中的0111与1000，编码是相邻的，但距离相差很大。

边界问题

比如现在是定位点1，根据geohash算法定位到定位点1所在的矩形块，返回附近的点就会得到点B。但是我们发现实际情况是，点A虽然不在定位点1所在的那个矩形区域，但是点A显然定位点更近。这就是geohash的边界问题，就是如果我们刚好处在矩形的边界处，那么离我们最近的点不一定是和我们在同一个矩形框中的点，很可能是旁边的矩形框中的某些点离我们更近。

解决方案：

目前比较通行的做法就是我们不仅获取当前我们所在的矩形区域，还获取周围8个矩形块中的点。那么怎样定位周围8个点呢？关键就是需要获取周围8个点的经纬度，那我们已经知道自己的经纬度，只需要用自己的经纬度减去最小划分单位的经纬度就行。因为我们知道经纬度的范围,有知道需要划分的次数，所以很容易就能计算出最小划分单位的经纬度。不理解的可以参考后面的代码实现部分。

精度对照表

应用场景

半径范围内热点位置数据列表（附近人、商店、美食等等）
范围内点位置数据GIS地图聚合（地图海量点数据聚合展示）
网格化运营数据统计

Google S2

原理

将经纬度S(lng, lat)使用坐标表示f(x,y,z)

将球面投影到外切正方形上，可以将f(x, y, z) => g(face, u, v),face是正方形的六个面，u，v对应的是六个面中的一个面上的x，y坐标。

上一步我们把球面上的球面矩形投影到正方形的某个面上，形成的形状类似于矩形，但是由于球面上角度的不同，最终会导致即使是投影到同一个面上，每个矩形的面积也不大相同。

通过一个二次变化 g(face,u,v) => h(face,s,t)。得到修正后的坐标

精度表

功能

经纬度转cellId
获取任意形状内的S2块
根据cellId获取其所在区域的4个点
判断2个形状的包含关系

@Test
    @DisplayName("经纬度转CellId")
    void lngLat2CellId() {
        double lat = 30.229212;
        double lng = 119.984686;
        int currentlevel = 15;
        GhPoint point = GeoTransformUtils.bd09_To_gps84(new GhPoint(lng, lat));

        S2LatLng s2LatLng = S2LatLng.fromDegrees(point.getLat(), point.getLng());
        S2CellId cellId = S2CellId.fromLatLng(s2LatLng).parent(currentlevel);
        log.info("cellId = {}", cellId);
    }

    @Test
    @DisplayName("经纬度=>多边形")
    void lnglatToPolygon() {
        double lat = 30.229212;
        double lng = 119.984686;
        int currentlevel = 15;
        GhPoint point = GeoTransformUtils.bd09_To_gps84(new GhPoint(lng, lat));

        S2LatLng s2LatLng = S2LatLng.fromDegrees(point.getLat(), point.getLng());
        S2CellId cellId = S2CellId.fromLatLng(s2LatLng).parent(14);
        S2Cell cell = new S2Cell(cellId);
        S2LatLng center = new S2LatLng(cell.getCenter());
        S2LatLng v_0 = new S2LatLng(cell.getVertex(0));
        S2LatLng v_1 = new S2LatLng(cell.getVertex(1));
        S2LatLng v_2 = new S2LatLng(cell.getVertex(2));
        S2LatLng v_3 = new S2LatLng(cell.getVertex(3));
    }

    @Test
    @DisplayName("获取所有的cell")
    void getAllCell() {
        double lat = 30.229212;
        double lng = 119.984686;
        int currentlevel = 14;
        GhPoint point = GeoTransformUtils.bd09_To_gps84(new GhPoint(lng, lat));
        S2LatLng s2LatLng = S2LatLng.fromDegrees(point.getLat(), point.getLng());
        S2CellId cellId = S2CellId.fromLatLng(s2LatLng).parent(14);
        S2RegionCoverer coverer = S2RegionCoverer.builder().setMaxLevel(15).setMinLevel(15).build();
        ArrayList<S2CellId> s2CellIds = Lists.newArrayList();
        coverer.getInteriorCovering(new S2Cell(cellId), s2CellIds);
    }

    @Test
    @DisplayName("是否包含")
    void polygonRelation() {
        double lat = 30.229212;
        double lng = 119.984686;
        int currentlevel = 14;
        GhPoint point = GeoTransformUtils.bd09_To_gps84(new GhPoint(lng, lat));
        S2LatLng s2LatLng = S2LatLng.fromDegrees(point.getLat(), point.getLng());
        S2CellId cellId = S2CellId.fromLatLng(s2LatLng).parent(currentlevel);
        S2Polygon s2Polygon = new S2Polygon(new S2Cell(cellId));
        GhPoint point_other = GeoTransformUtils.bd09_To_gps84(new GhPoint(119.982228, 30.228535));
        S2LatLng s2LatLng_other = S2LatLng.fromDegrees(point_other.getLat(), point_other.getLng());
        S2CellId cellId_other = S2CellId.fromLatLng(s2LatLng_other).parent(9);
        S2Polygon s2Polygon_other = new S2Polygon(new S2Cell(cellId_other));
        boolean contains = s2Polygon_other.contains(s2Polygon);
        System.out.println("是否包含：" + contains);
    }

参考

文档

s2文档

演示工具

百度地图添加覆盖物网址

https://lbsyun.baidu.com/jsdemo.htm#dAddClearOverlay

本地项目地址

http://localhost/static/geohash_click.html

上一页杂记

最后更新于1年前