一 使用场景     

 目前基于某个位置查附近的人，附近的商家等等，查出来的结果添加距离，或者查附近多大范围内的人或者商家，然后按距离排序已经是IT界一个很通用的功能了。

二 距离计算搜索(百万点集以下)

2.1 球的定义

2.2 两点之间距离计算

     

2.2.1两点间直线距离

    欧式距离计算法: 欧式距离（Euclidean Distance）是欧几里得空间中两点之间的直线距离，广泛应用于二维、三维及更高维度的空间计算12。其核心思想基于勾股定理，通过平方差的累加与平方根运算实现距离度量.

       

2.2.2余弦定理算球面距离

• 球面距离

        

如果AOB三点共线的时候我们可能会有N个这样的圆，如果AOB三点不共线的时候我们要算的是劣圆弧长。弧度数的概念实际就是用弧长除以半径称作单位半径弧长。

2.2.2.1认识地球

      

 本初子午线以东叫东经，以西叫西经。赤道以北叫北纬，赤道以南叫南纬。

2.2.2.2 正余弦定理

        正弦定理

       

• 余弦定理

       

2.2.2.3 同经度弧长计算

         

2.2.2.4 同维度弧长计算

           

2.2.2.5 经纬度都不同

          

   如上图所示，如何求出A1A2的弧长.也就是我们说的球面的劣弧距离。也就是我们地球上任何物体之间的球面距离。

          

 4）AB平方的求解过程

              

5)化解

               

2.2.3 Haversine模型算球面距离(哈弗赛)

2.2.4 Vincenty模型算球面距离

2.3 距离计算搜索落地

    现在几乎所有的O2O应用中都会存在“按范围搜素、离我最近、显示距离”等等基于位置的交互，那这样的功能是怎么实现的呢？本文提供的实现方式，适用于所有数据库。

为了方便下面说明，先给出一个初始表结构

CREATE TABLE `customer` (
    `id` INT(11) UNSIGNED NOT NULL AUTO_INCREMENT COMMENT '自增主键',
    `name` VARCHAR(5) NOT NULL COMMENT '名称',
    `log` DOUBLE(9,6) NOT NULL COMMENT '经度',
    `lat` DOUBLE(8,6) NOT NULL COMMENT '纬度',
    PRIMARY KEY (`id`)
)
COMMENT='商户表'
CHARSET=utf8mb4
ENGINE=InnoDB

实现过程主要分为四步： 

1. 搜索 
在数据库中搜索出接近指定范围内的商户，如：搜索出1公里范围内的。 
2. 过滤 
搜索出来的结果可能会存在超过1公里的，需要再次过滤。如果对精度没有严格要求，可以跳过。 
3. 排序 
距离由近到远排序。如果不需要，可以跳过。 
4. 分页 
如果需要2、3步，才需要对分页特殊处理。如果不需要，可以在第1步直接SQL分页。

第1步数据库完成，后3步应用程序完成。

2.3.1 搜索

<dependency>
    <groupId>org.locationtech.spatial4j</groupId>
    <artifactId>spatial4j</artifactId>
    <version>0.8</version>
</dependency>

2.3.1.1 普通列表距离计算

1. Law of Cosines(余弦定理)

     本计算式理论模型为余玄定理,平面计算通过三角函数,误差较大,底层原理就是采用基本的球面距离计算,底层实现就是采用的余弦定理。

double lat1 = DistanceUtils.toRadians(0.0);
double lon1 = DistanceUtils.toRadians(0.0);
double lat2 = DistanceUtils.toRadians(4.0);
double lon2 = DistanceUtils.toRadians(3.0);
//余玄定理，平面计算DistanceUtils.distLawOfCosinesRAD
double distLawOfCosinesRad = DistanceUtils.distLawOfCosinesRAD(lat1,lon1,lat2,lon2);
//弧度转换为角度
double distLawOfCosinesDegrees = DistanceUtils.toDegrees(distLawOfCosinesRad);
System.out.println("距离(度):"+distLawOfCosinesDegrees);

2 Haversine球面模型
本计算式理论模型为椭球模型，计算时会多次迭代，理论计算精度较高,入参返回

//椭球模型
DistanceUtils.distVincentyRAD
DoubledistVincentyDegree = DistanceUtils.toDegrees(DistanceUtils.distVincentyRAD(lat1,lon1,lat2,lon2));
System.out.println("椭球模型,距离(度):"+distVincentyDegree);

3 .Vincenty椭球模型

  本计算式选取地球模型为球模型，以赤道半径为基准，故计算时纬度越高误差会越大，但胜在快速,入参返回(单位弧度)

//球模型DistanceUtils.distHaversineRAD
double distHaversineDegree = DistanceUtils.toDegrees(DistanceUtils.distHaversineRAD(lat1,lon1,lat2,lon2));
System.out.println("球模型,距离(度):"+distHaversineDegree);

 4 结果对比

余玄定理,距离(度):4.998536879937023

椭球模型,距离(度):4.998536879936981

球模型,距离(度):4.998536879936979

2.3.1.2列表范围查找 

     customer表中使用两个字段存储了经度和纬度，如果提前计算出经纬度的范围，然后在这两个字段上加上索引，那搜索性能会很不错。 
那怎么计算出经纬度的范围呢？已知条件是移动设备所在的经纬度，还有满足业务要求的半径，这很像初中的一道平面几何题：给定圆心坐标和半径，求该圆外切正方形四个顶点的坐标。而我们面对的是一个球体，可以使用spatial4j来计算。

// 移动设备经纬度
double lon = 116.312528, lat = 39.983733;
// 千米
int radius = 1;
 
SpatialContext geo = SpatialContext.GEO;
Rectangle rectangle = geo.getDistCalc().calcBoxByDistFromPt(
        geo.makePoint(lon, lat), radius * DistanceUtils.KM_TO_DEG, geo, null);
System.out.println(rectangle.getMinX() + "-" + rectangle.getMaxX());// 经度范围
System.out.println(rectangle.getMinY() + "-" + rectangle.getMaxY());// 纬度范围

计算出经纬度范围之后，SQL是这样：

SELECT id, name
FROM customer
WHERE (lon BETWEEN ? AND ?) AND (lat BETWEEN ? AND ?);
需要给lon、lat两个字段建立联合索引：
INDEX `idx_lon_lat` (`lon`, `lat`)

三 GeoHash 搜索(百万点集以上)

   为了提高数据量的检索，我们可能需要对地球表面的X公里范围内空间的点做索引，一提到索引，大家脑子里马上浮现出B树索引，因为大量的数据库（如MySQL、oracle、PostgreSQL等）都在使用B树。B树索引本质上是对索引字段进行排序，然后通过类似二分查找的方法进行快速查找，即它要求索引的字段是可排序的，一般而言，可排序的是一维字段，比如时间、年龄、薪水等等。但是对于空间上的一个点（二维，包括经度和纬度），如何排序呢？又如何索引呢？解决的方法很多，下文介绍一种方法来解决这一问题。

如果能通过某种方法将二维的点数据转换成一维的数据，那样不就可以继续使用B树索引了嘛。那这种方法真的存在嘛，答案是肯定的。目前很火的GeoHash算法就是运用了上述思想，下面我们就开始GeoHash之旅吧。

3.1 GeoHash 核心原理

                              

    GeoHash将二维的经纬度转换成字符串，比如下图展示了9个区域的GeoHash字符串，分别是WX4ER，WX4G2、WX4G3等等，每一个字符串代表了某一矩形区域。也就是说，这个矩形区域内所有的点（经纬度坐标）都共享相同的GeoHash字符串，这样既可以保护隐私（只表示大概区域位置而不是具体的点），又比较容易做缓存，比如左上角这个区域内的用户不断发送位置信息请求餐馆数据，由于这些用户的GeoHash字符串都是WX4ER，所以可以把WX4ER当作key，把该区域的餐馆信息当作value来进行缓存，而如果不使用GeoHash的话，由于区域内的用户传来的经纬度是各不相同的，很难做缓存。

       1) 字符串越长，表示的范围越精确。如图所示，5位的编码能表示10平方千米范围的矩形区域，而6位编码能表示更精细的区域（约0.34平方千米）字符串相似的表示距离相近（特殊情况后文阐述），这样可以利用字符串的前缀匹配来查询附近的目标位置信息。如上图，一个在城区，一个在郊区，城区的GeoHash字符串之间比较相似，郊区的字符串之间也比较相似，而城区和郊区的GeoHash字符串相似程度要低些。

       2) 通过上面的介绍我们知道了GeoHash就是一种将经纬度转换成字符串的方法，并且使得在大部分情况下，字符串前缀匹配越多的距离越近，只需要将所在位置经纬度转换成GeoHash字符串，通过前缀匹配的GeoHash字符串进行前缀匹配，匹配越多的距离越近。

        3) 由于GeoHash是将区域划分为一个个规则矩形，并对每个矩形进行编码，这样在查询附近的事务信息时会导致以下问题，比如距型WX4GF1和矩形WX4GD1 他们两个距离很近，但是由于我们只取了五位导致了我们在计算WX4GD1附近的距离的时候容易忽略掉WX4GF1中的距离,这个问题往往产生在边界处,那对于这个问题我们如何解决呢?

       4) 那就是我们查询时，除了使用定位点的GeoHash编码进行匹配外，还使用周围8个区域的GeoHash编码，这样可以避免这个问题。 比如我们找WX4GD2附近的区域的事物的时候，我们就需要扫描到它周边的8个相似区域的事务。

geohash length	width	height
1	5,009.4km	4,992.6km
2	1,252.3km	624.1km
3	156.5km	156km
4	39.1km	19.5km
5	4.9km	4.9km
6	1.2km	609.4m
7	152.9m	152.4m
8	38.2m	19m
9	4.8m	4.8m
10	1.2m	59.5cm
11	14.9cm	14.9cm
12	3.7cm	1.9cm

     上表来自维基百科，可以看出，当geohash base32编码长度为8时，精度在19米左右，而当编码长度为9时，精度在5米左右，编码长度需要根据数据情况进行选择。 

3.2 GeoHash 算法

3.2.1 经纬度编码

     根据经纬度计算GeoHash二进制编码地球纬度区间是[-90,90]， 北海公园的纬度是39.928167，可以通过下面算法对纬度39.928167进行逼近编码:

1）区间[-90,90]进行二分为[-90,0),[0,90]，称为左右区间，可以确定39.928167属于右区间[0,90]，给标记为1；

2）接着将区间[0,90]进行二分为 [0,45),[45,90]，可以确定39.928167属于左区间 [0,45)，给标记为0；

3）递归上述过程39.928167总是属于某个区间[a,b]。随着每次迭代区间[a,b]总在缩小，并越来越逼近39.928167；

4）如果给定的纬度x（39.928167）属于左区间，则记录0，如果属于右区间则记录1，这样随着算法的进行会产生一个序列1011100，序列的长度跟给定的区间划分次数有关。

3.2.2编码转组码

　    通过上述计算，纬度产生的编码为10111 00011，经度产生的编码为11010 01011。偶数位放经度，奇数位放纬度，把2串编码组合生成新串：11100 11101 00100 01111。

　　最后使用0-9、b-z（去掉a, i, l, o）这32个字母进行base32编码，首先将11100 11101 00100 01111转成十进制，对应着28、29、4、15，十进制对应的编码就是wx4g。同理，将编码转换成经纬度的解码算法与之相反，具体不再赘述。

如图所示，如我们发现纬度放在偶数位，经度放在奇数位那为什么是这样的呢

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

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

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

     我们已经知道现有的GeoHash算法使用的是Peano空间填充曲线，这种曲线会产生突变，造成了编码虽然相似但距离可能相差很大的问题，比如上面比如0111与1000，编码是相邻的，但距离相差很大。那其实这就是我们所说的噪声点，在查询的时候我们是无法避免的，智能在我们搜索出所有的点以后，再做一个实际距离的计算，来去除噪声。

 3.3 GeoHash搜索落地

3.3.1 引入包

<dependency>
    <groupId>org.locationtech.spatial4j</groupId>
    <artifactId>spatial4j</artifactId>
    <version>0.8</version>
</dependency>

 3.3.2 建表

    在我们了解了 geohash。geohash算法能把二维的经纬度编码成一维的字符串，它的特点是越相近的经纬度编码后越相似，所以可以通过前缀like的方式去匹配周围的商户。 customer表要增加一个字段，来存储每个商户的geohash编码，并且建立索引。 

CREATE TABLE `customer` (
    `id` INT(11) UNSIGNED NOT NULL AUTO_INCREMENT COMMENT '自增主键',
    `name` VARCHAR(5) NOT NULL COMMENT '名称' COLLATE 'latin1_swedish_ci',
    `logtitude` DOUBLE(9,6) NOT NULL COMMENT '经度',
    `latitude` DOUBLE(8,6) NOT NULL COMMENT '纬度',
    `geo_code` CHAR(12) NOT NULL COMMENT 'geohash编码',
    PRIMARY KEY (`id`),
    INDEX `idx_geo_code` (`geo_code`)
)
COMMENT='商户表'
CHARSET=utf8mb4
ENGINE=InnoDB
;

 3.3.3 设计搜索

        1) 搜索

      在新增或修改一个商户的时候，维护好geo_code，那geo_code怎么计算呢？spatial4j也提供了一个工具类,默认精度是12位。这个存储做好后，就可以通过geo_code去搜索了。

GeohashUtils.encodeLatLon(lat, lon)

   假设我们的需求是1公里范围内的商户，geo_code的长度设置为5就可以了，GeohashUtils.encodeLatLon(lat, lon, 5)。计算出移动设备经纬度的geo_code之后，SQL是这样：

      SELECT id, name FROM customer WHERE geo_code LIKE CONCAT(?, '%'); 这样会比区间查找快很多，并且得益于geo_code的相似性，可以对热点区域做缓存。但这个也会有一个问题，就是我们上面提到的边界距离计算的问题，容易算错，解决办法也如上讲解说到的，我们需要取周边8个 geo_code,那怎么计算出周围8个网格的geohash呢.

<dependency>
    <groupId>ch.hsr</groupId>
    <artifactId>geohash</artifactId>
    <version>1.3.0</version>
</dependency>

引入获取邻接的8个矩形
double lon = 116.312528, lat = 39.983733;
GeoHash geoHash = GeoHash.withCharacterPrecision(lat, lon, 10);
// 当前
System.out.println(geoHash.toBase32());
// N, NE, E, SE, S, SW, W, NW
GeoHash[] adjacent = geoHash.getAdjacent();
for (GeoHash hash : adjacent) {
    System.out.println(hash.toBase32());
}
最终我们的sql变成了这样：
SELECT id, name
FROM customer
WHERE geo_code LIKE CONCAT(?, '%')
OR geo_code LIKE CONCAT(?, '%')
OR geo_code LIKE CONCAT(?, '%')
OR geo_code LIKE CONCAT(?, '%')
OR geo_code LIKE CONCAT(?, '%')
OR geo_code LIKE CONCAT(?, '%')
OR geo_code LIKE CONCAT(?, '%')
OR geo_code LIKE CONCAT(?, '%')
OR geo_code LIKE CONCAT(?, '%');

        原来的1次查询变成了9次查询，性能肯定会下降，这里可以优化下。还用上面的需求场景，搜索1公里范围内的商户，从上面的表格知道，geo_code长度为5时，网格宽高是4.9KM，用9个geo_code查询时，范围太大了，所以可以将geo_code长度设置为6，即缩小了查询范围，也满足了需求。还可以继续优化，在存储geo_code时，只计算到6位，这样就可以将sql变成这样： SELECT id, name FROM customer WHERE geo_code IN (?, ?, ?, ?, ?, ?, ?, ?, ?); 这样将前缀匹配换成了直接匹配，速度会提升很多。

 2 ) 过滤

上面两种搜索方式，都不是精确搜索，只是尽量缩小搜索范围，提升响应速度。所以需要在应用程序中做过滤，把距离大于1公里的商户过滤掉。计算距离同样使用spatial4j。

// 移动设备经纬度
double lon1 = 116.3125333347639, lat1 = 39.98355521792821;
// 商户经纬度
double lon2 = 116.312528, lat2 = 39.983733;
 
SpatialContext geo = SpatialContext.GEO;
double distance = geo.calcDistance(geo.makePoint(lon1, lat1), geo.makePoint(lon2, lat2)) 
    * DistanceUtils.DEG_TO_KM;
System.out.println(distance);// KM

过滤代码就不写了，遍历一遍搜索结果即可。

3) 排序

    把这些计算出来的经纬度和1km距离范围内的点，通过redis的zset集合添加到集合中，距离作为score值进行添加。

4) 分页

    zrange借助redis进行内存分页查询。当搜索围从1km变成2km，或者5km的时候重复步奏1,2,3,4

四 地图围栏算法

五 地图空间索引

六 地图点聚合优化