引言

你是否好奇, 定位系统是如何快速判断你在哪个城市的.

为什么有的想吃的外卖无法送达你的小区, 隔壁的小区又可以送达.

为什么你的儿童手表或者宠物追踪环, 有超出范围预警, 这一切是如何做到的呢 ?

这都离不开地理围栏的作用.

本文尽量通俗的从数学原理讲起, 然后从工程优化和算法优化两个层面讲解是如何提高判断效率的. 最后会介绍我们用地理围栏做了什么.

什么是地理围栏

地理围栏（Geo-fencing）是LBS的一种新应用，就是用一个虚拟的栅栏围出一个虚拟地理边界。当手机进入、离开某个特定地理区域，或在该区域内活动时，手机可以接收自动通知和警告。有了地理围栏技术，位置社交网站就可以帮助用户在进入某一地区时自动登记。与LBS以某一区域为中心覆盖周边3公里或5公里半径的定位方式不同，地理围栏更侧重于对区域边界的界定，不再是以某点为圆心向外等距离画圆，而是准确勾勒出小区、写字楼等特定坐标的实际形状、区域及面积。

接下来看下如何根据地理围栏判定定位用户位置。

地理围栏定位判断算法

每个城市作为一条记录，包含了该城市的地理围栏数据(抽象成一个多边形区域)，对于给定点（纬度，经度）定位所属城市。问题就可以抽象为判断一个点是否属于一个多边形区域内部。

如何判断一个点是否在多边形内？

判断一个点是否在多边形内是处理空间数据时经常面对的需，例如GIS中的点选功能、根据多边形边界筛选出位于多边形内的点、筛选不在多边形内的点等，判定一个点是否在多边形内有几种不同的方法

射线法：从判断点向某个统一方向作射线，依交点个数的奇偶判断；
转角法：按照多边形顶点逆时针顺序，根据顶点和判断点连线的方向正负（设定角度逆时针为正）求和判断；
夹角和法：求判断点与所有边的夹角和，等于360度则在多边形内部；
面积和法：求判断点与多边形边组成的三角形面积和，等于多边形面积则点在多边形内部；

其中面积和法涉及多个面积的计算，比较复杂；夹角和法以及转角法用到角度计算，会涉及反三角函数，计算开销比较大；射线法主要是循环多边形的每条边进行求交运算，但大部分边可以通过简单坐标比对直接排除，相对是比较好的一种方法，下面主要是分析如何利用射线法进行定位判定。

射线法的具体实现

射线法的定义

射线法其实就是从判断点开始，向左边或者右边(下面我们讲解的都是围绕向右做射线的情况包括代码和图示)水平的做一条射线，判断该射线穿过多边形的次数，如果为奇数则认为该点在区域内部，如果为偶数则认为该点在区域外部。该方法对于复合多边形也能正确的判断。

举个例子：比如现在我们知道一个点的经纬度信息，以及某个省的地理围栏，接下来我们如何判断这个点是否在这个省内？

我们需要把每个省的地理围栏数据连线组成的多边形想象成在一个二维平面内，依次判断从判定点沿x轴向右引出的射线与这个多边形的每条边是否相交。具体可参考下图简单介绍了射线法的定义，接下来让我们看下具体如何实现相交点个数的判定。

射线法的实现

在这之前我们先来看下在平面坐标轴中的多边形和判定点的相交情况，如下图所示：从图中可以看出对于标注中的特殊的边我们可以首先排除掉，不参与相交判断。具体可以分为如下几种情况：

边与射线平行或者重合：这时边的两个端点的Y坐标相等。如图中的 1、5 两条边
边在射线上边或下边代表：这时边的两个端点的Y坐标均大于或小于判定点的Y坐标。如图中的 2、8 两条边
边在射线左边：这时边的两个端点的X坐标均小于判定点的Y坐标。如图中的 10、11 两条边
边的上端点或者下端点与射线相交：这时边的起始点中的一个点的Y左标等于判定点的Y左标，另一个点的Y坐标大于判定点Y坐标(小于判定点的视为相交)。如图中的 4 这条边

我们首先排除了上面几种特殊情况，那么对于3、7、6、9几条边我们如何判定某个边与射线的交点是在判定点的X坐标的左边还是右边呢？

具体的判断需要用到相似三角形的知识，首先来看下什么是相似三角形：

三角分别相等，三边成比例的两个三角形叫做相似三角形。相似三角形的面积比等于边长比的平方。

我们需要运用相似三角形的性质，求出射线与边的交点的X坐标，然后判断交点的X坐标是否大于判定点的X坐标，大于则相交，交点个数加1，否则不相交。

计算公示如下图一，同时我们以 9 这条边为例判断其是否与判定点向右的射线相交如下图二。

如下是基于python的是代码实现

判定边与射线是否相交的实现方法如下：

# 射线与边是否有交点
def isRayIntersectsSegment(poi, sPoi, ePoi):  # [x,y] [lng,lat]
    # 排除 与射线平行、重合、是一个点的情况
    if sPoi[1] == ePoi[1]:
        return False
    # 排除 线段在射线上边
    if sPoi[1] > poi[1] and ePoi[1] > poi[1]:
        return False
    # 排除 线段在射线下边
    if sPoi[1] < poi[1] and ePoi[1] < poi[1]:
        return False
    # 排除 交点为下端点
    if sPoi[1] == poi[1] and ePoi[1] > poi[1]:
        return False
    # 排除 交点为上端点
    if ePoi[1] == poi[1] and sPoi[1] > poi[1]:
        return False
    # 排除 线段在射线左边
    if sPoi[0] < poi[0] and ePoi[0] < poi[0]:
        return False
    # 求交，相似多边形性质
    xseg = ePoi[0] + (poi[1] - ePoi[1]) / (sPoi[1] - ePoi[1]) * (sPoi[0] - ePoi[0])
    # 交点在射线起点的左侧
    if xseg < poi[0]:
        return False
    #  排除上述情况之后真正有效的点
    return True

针对多多边形MultiPolygon的最终判断代码如下：

# MultiPolygon
# 输入：点，多边形三维数组 
# poly =  [[[116.637525000110358, 41.060467999820851], [116.713636000435855, 40.910223999305245]]] 三维数组
def poiInPoly(poi, poly, tolerance=0.0001):
    # 先判断点是否在外包矩形内
    if not isPoiWithinBox(poi, [[0, 0], [180, 90]], tolerance):
        return False
    # 交点个数
    count = 0
    # 逐个二维数组进行判断
    for epoly in poly:
        # [0,len-1]
        for i in range(len(epoly) - 1):
            s_poi = epoly[i]
            e_poi = epoly[i + 1]
            if isRayIntersectsSegment(poi, s_poi, e_poi):
                # 有交点加1
                count += 1
    return count % 2 != 0

城市边界为geojson格式，关于geojson的介绍可以参考这篇文章：GEOJSON标准格式。

这样我们就完成了点是否在多边形的判断。当多边形数目较少时，我们可以依次遍历每一个多边形（暴力遍历法），然后用射线法进行判断，这样效率也很高。但是当多边形数目较多时，响应时间会比较长。比如我们的需求是需要根据用户上报的经纬度信息判断用户所在的行政区域精确到乡镇级别，每天大概上报40万条定位数据。目前全国目前大概有44821个乡镇级行政区，也就是有4万多个地理围栏，那么对每个定位使用射线发判断的时间复杂度就是O(n)。40万条定位数据分别判定属于这4万多个地理围栏中的哪一个的时间复杂度就是O(n2)，显然这样的时间复杂度效率是比较低的。接下来介绍一种优化方法，对给定的所有地理围栏数据建立RTree索引

判定算法的优化

基于R树的优化

暴力遍历法效率低下的主要原因是与每一个多边形都进行了射线法判断，如果能减少射线法的调用次数性能就能提升，因此我们可以通过粗筛的方法快速找到符合条件的少量多边形，然后对粗筛后的多边形使用射线法判断，这样射线法的执行次数会大大减少，执行效率也能大大提高。

如何进行筛选呢？

对于一维数据我们常使用索引的方法，比如通过B树(B树是平衡多路查找树,如下图一)索引找到某一个范围区间段，然后对此范围区间段进行遍历查找,关于B树可以参考这篇文章：B树、B+树、B*树；
对于二维空间数据常使用空间索引的方法，比如通过R(Rectangle)树(R树是一棵用来存储高维数据的平衡树，如下图二)找到范围区间内的多边形，然后对此范围内的多边形进行精确判断； R树主要运用了空间分割的理念。采用了一种称为MBR(Minimal Bounding Rectangle)的方法，可以将它成为“最小外接矩形”。从叶子结点开始用矩形（rectangle）将空间框起来，结点越往上，框住的空间就越大，以此对空间进行分割。 R树无论是叶子结点还是非叶子结点，它们都对应着一个矩形。树形结构上层的结点所对应的矩形能够完全覆盖它的孩子结点所对应的矩形。根结点也唯一对应一个矩形，而这个矩形是可以覆盖所有我们拥有的数据信息在空间中代表的点的。

什么最小外接矩形？

最小外接矩形是指以二维坐标表示的若干二维形状（例如点、直线、多边形）的最大范围，即以给定的二维形状各顶点中的最大横坐标、最小横坐标、最大纵坐标、最小纵坐标定下边界的矩形。如下图这样的一个矩形包含给定的二维形状，且边与坐标轴平行。 行政区域判断

了解了什么是R树以及最小外包矩形的概念，那么我们如何基于R树进行行政区域查询的优化呢？

例如现在我们知道一个点的经纬度，并且现在我们知道中国、六大行政区、各省、各市、各县区的地理围栏，判定它是的定位位置精确到县区级。

首先先建立中国、六大行政区、各省、各市、各县以及各乡镇的最小外接矩形
把中国作为整个R树，首先判断是否在中国最小外接矩形内，并以中国作为R树的根结点
然后判断是在那个大区，各大区作为第一层的结点
判定属于华北大区内，接着判断是属于那个华北大区的那个省，华北大区各省作为第二层结点。
判定属于河南省内，接着判断是属于那个河南省的那个市，河南省各市作为第三层结点。
判定属于郑州市内，最后依次遍历郑州市的各区，判断判定节点是在那个区内比如在金水区。

具体如下图所示这样我们就完成了基于R树的索引，R树平均查询复杂度为O(Log(n))，n为多边形个数，优与O(n)的时间复杂度。

多边形边数较多的情况

大多数应用的地理围栏多边形都比较简单，但有时也会遇到一些特别复杂的多边形，比如单个多边形的边数就超过十几万条，这时候对此复杂多边形执行一次射线法也非常耗时（因为射线法时间复杂度为O(n)，n为多边形边数）。

如何提高对复杂多边形执行射线法的计算效率呢？同样使用R树索引，对边数较多（如超过1万）的多边形的边再单独进行R树索引，首先对多边形的每条边构建最小外包矩形，然后在这些最小外包矩形基础上构建R树索引，这样射线法求交点的时候首先通过R树判断射线是否与外包矩形相交，最后对R树粗筛后的边进行精确求交判断，时间复杂度从O(n)降到O(Log(n))，可大大提高了计算效率。

例子

以上就完成了基于地理围栏的定位判断，下面是利用此算法判断一批经纬度所在的行政区域是否正确的例子：

import requests
import json
import time
import pandas

region_api = 'https://maf.meituan.com/region'
region_list = []
root_url = ''


def getRegionInfo():
    reg_list = []
    url = ''
    params = {"key": "c01bc331-ce26", "region": 100000, "subDistrict": 3, "regionType": 2,
              "returnGeometry": 'true'}
    response = requests.get(url=region_api, params=params)
    resData = json.loads(response.text)
    if int(resData['status']) == 200:
        reg_list = resData['result']['admin_list'][0]['children']
        url = resData['result']['rootUrl']
    return {"reg_list": reg_list, "root_url": url}


# 点是否在外包矩形内
# sbox=[[x1,y1],[x2,y2]]
def isPoiWithinBox(poi, sbox, toler=0.0001):
    if sbox[0][0] < poi[0] < sbox[1][0] and sbox[0][1] < poi[1] < sbox[1][1]:
        return True
    if toler > 0:
        pass
    return False


# 射线与边是否有交点
# [x,y] [lng,lat]
def isRayIntersectsSegment(poi, sPoi, ePoi):
    # 排除 与射线平行、重合、是一个点的情况
    if sPoi[1] == ePoi[1]:
        return False
    # 排除 线段在射线上边
    if sPoi[1] > poi[1] and ePoi[1] > poi[1]:
        return False
    # 排除 线段在射线下边
    if sPoi[1] < poi[1] and ePoi[1] < poi[1]:
        return False
    # 排除 交点为下端点
    if sPoi[1] == poi[1] and ePoi[1] > poi[1]:
        return False
    # 排除 交点为上端点
    if ePoi[1] == poi[1] and sPoi[1] > poi[1]:
        return False
    # 排除 线段在射线左边
    if sPoi[0] < poi[0] and ePoi[0] < poi[0]:
        return False
    # 求交
    # xseg = ePoi[0] - (ePoi[0] - sPoi[0]) * (ePoi[1] - poi[1]) / (ePoi[1] - sPoi[1])
    xseg = ePoi[0] + (poi[1] - ePoi[1]) / (sPoi[1] - ePoi[1]) * (sPoi[0] - ePoi[0])
    # 交点在射线起点的左侧
    if xseg < poi[0]:
        return False
    #  排除上述情况之后真正有效的点
    return True


# MultiPolygon
# 输入：点，多边形三维数组
# poly =  [[[116.637525000110358, 41.060467999820851], [116.713636000435855, 40.910223999305245]]] 三维数组
def poiInPoly(poi, poly, tolerance=0.0001):
    # 先判断点是否在外包矩形内
    if not isPoiWithinBox(poi, [[0, 0], [180, 90]], tolerance):
        return False
    # 交点个数
    count = 0
    # 逐个二维数组进行判断
    for epoly in poly:
        # [0,len-1]
        for i in range(len(epoly) - 1):
            s_poi = epoly[i]
            e_poi = epoly[i + 1]
            if isRayIntersectsSegment(poi, s_poi, e_poi):
                # 有交点加1
                count += 1
    return count % 2 != 0


def getRegions(poi, regions, region_res):
    for index, region in enumerate(regions):
        geofence_url = root_url + '/' + region['geoUrl']
        response = requests.get(geofence_url)
        try:
            resData = json.loads(response.text)
            coordinates = []
            for coordinateList in resData['coordinates']:
                coordinates.extend(coordinateList)
            isInPoly = poiInPoly(poi, coordinates)
            if isInPoly:
                region_res.append(region['name'])
                getRegions(poi, region['children'], region_res)
                break
        except:
            print(response.text)


def paseCityLocs():
    file_path = 'cityData1.csv'
    city_file = open(file_path, 'r')
    reader = pandas.read_csv(city_file, nrows=20000000)
    correctNum = 0
    num = 70
    for index, row in reader.iterrows():
        if index < num:
            # 获取定位行政区域
            # line 数据格式 115.382268,25.136305,江西省,赣州市,安远县
            print('原始数据', index, row.values)
            region_res = []
            getRegions([float(row['longitude']), float(row['latitude'])], region_list, region_res)
            print('通过射线法判断', region_res)
            try:
                # 排除没有县的情况
                if row['province'] == region_res[0] \
                        and (pandas.isnull(row['district']) or row['district'] == region_res[2]):
                    correctNum = correctNum + 1
                else:
                    print('不一致')
                print()
            except:
                print('error', region_res)

    print('总数:', num)
    print('一致数总计:', correctNum)
    print('一致率:', correctNum / num)


# 首先获取所有的行政区域到县级以及获取地理围栏的url
regionInfo = getRegionInfo()
region_list = regionInfo['reg_list']
root_url = regionInfo['root_url']
# 处理经纬度数据
start = time.time()
paseCityLocs()
end = time.time()
print('用时:', (end - start))

对200条数据进行测试，一致率99%。其中判断和原始数据不一致的最终确认(查询地址)后发现通过算法判断的是正确的，demo中获取行政区域及其地理围栏用的美团地图提供的省市区县编码API，最终运行结果如下图

JS实现点在多边形内的判断

function pointInPoly(point: number[], points: number[][]) {
  let isContain = false;
  let end = points.length - 1;
  for (let i = 0; i < points.length; i++) {
    const pointX = point[0];
    const pointY = point[1];
    const pointIX = points[i][0];
    const pointIY = points[i][1];
    const pointEndX = points[end][0];
    const pointEndY = points[end][1];
    if (pointY > pointIY && pointY < pointEndY || pointY < pointIY && pointY > pointEndY) {
      if ((pointY - pointIY) / (pointEndY - pointIY) * (pointEndX - pointIX) + pointIX < pointX) {
        isContain = !isContain;
      }
    }
    end = i;
  }
  return isContain;
}
pointInPoly([7, 5], [[2, 1], [8, 1], [8, 6], [6, 6], [6, 3], [4, 3], [4, 5], [2, 5]]);

算法篇：地理围栏算法解析

引言