射线法进行地理空间坐标计算

移动端应用开发过程中经常需要计算地理空间关系，例如：快递员/外卖员是否已抵达小区附近，办公打卡签到时是否在公司周边，共享单车是否停放在指定区域或禁停区等等。

前面我们介绍过如何使用MySQL进行地理坐标计算，不过这个方案存在几个局限性：

所有地理数据必须入库存储，有时候我们可能只是单纯想内存计算一下；
MySQL不支持曲线计算，只能通过线段拟近，存在一定误差；
对于航空空域等计算场景，当地理坐标跨越了换日线时，无法直接用几何图形表示地理关系。

上面第三点是作者在做实际项目时遇到的问题，因此在这里介绍另一种计算地理空间关系的方案，以及个人在实践过程中踩过的一些坑。

一、地理坐标系统

地理空间计算一般有两种坐标系统：地理坐标系以及投影坐标系。

地理坐标系就是按球面坐标计算地理关系，坐标单位是经纬度，通过经纬线将参考椭球划分格网以度量其表面特定位置。当我们以东/西经和南/北纬指代位置时，我们就在使用地理坐标系。

投影坐标系是平面坐标系，参考平面是水平面，坐标单位是米、千米等。投影坐标系的基础是地理坐标系，它定义了如何将球面上的一个点映射到平面坐标上。根据投影方法的不同，各类投影坐标的平面图形也不尽相同，如下图所示：

圆锥投影
圆柱投影
平面投影
极方位投影

因此不管是地理坐标系还是投影坐标系，都可以通过一定的转换方式将其转化为平面坐标计算。

二、场景分解

地理空间计算主要有几种常见类型：相交、包含、最短距离等等，这里为简化计算，选取点是否被平面包含这一例子进行分析。

简单地理空间计算

计算点是否被特定地理平面所包含，一个简单的方法是利用图形API库，例如Java里的2D API。在内存中构建对应图形对象，使用关系函数即可求得结果。不过这里的前提是需要将地理数据转换为图形对象，对于数据量较小的场景可以快速处理，例如上面办公打卡签到，只需要对公司大楼图形坐标进行转换即可。

复杂计算场景

考虑另一个场景，飞机航行过程中，如何判断当前处于哪些空域范围内。国际上将各个地理空间范围规划为不同的空域，并赋予不同的管理级别，例如我国民航可用空域面积319.53万平方公里，占比32%，其余为军用或受控制空域。实时掌握飞机当前所处空域，能够帮助规避潜在风险。与上面例子不同的是，这里我们无法事先知道有哪些空域需要参与计算，并且空域数量庞大、形状不规则（如下图所示），如果使用轮询的方式套搬上面的方法，会导致计算量过大，无法保证实时性。

这里我们对问题进行拆解，第一步对空域数据进行检索，排除非相关部分；第二步对检索得到的空域进行坐标计算，判断是否为包含关系。

三、地理数据检索

通常使用R树索引对地理数据进行检索，这是一种非精确的索引，利用图形最小包围盒MBR进行快速过滤。

由于R树索引的特性，构建索引时图形粒度越小，检索的效率就会越高。如果将空域作为一个整体构建索引，在检索时会返回较多不相关数据，导致过滤率较低。因此我们转换另一个思路，对构成空域的各个线段进行索引，提升过滤率。

四、射线法计算空间关系

索引细化后，我们无法直接判断点与原空域的空间关系，这里我们改用几何光学里的射线法进行空间计算。
射线法的原理是：从目标点向无穷远处发射一条直线，如果直线与多边形的交点个数为奇数，则该点位于多边形的内部。

结合R树检索与射线法计算过程后，整体算法调整如下：

从目标点（x, y）到水平无穷远处（∞, y）构建一条直线A；
利用R树索引检索所有与直线A相交的空域线段；
使用射线法判断各个空域相交线段的数量，如果为奇数，则目标点位于该空域内。

以上算法的优点是最大限度的利用了R树索引，索引得到的结果即为第三步射线法计算所需要的全集。在实际使用过程中该算法表现也符合预期，例如空域数量10W时，对应空域线段数量在200W左右，经过R树索引过滤后相关线段只有200+。

五、特殊处理

射线法特殊场景

使用射线法判断相交时，有几类特殊场景需要处理，具体如下：

射线与线段交点刚好为两线段的交接处；

此时射线同时穿过两线段，但实际交点只有一个。因此在进行相交判断时，需要按左开右闭图形进行计算，同时构成图形的线段数据需要按照有向图排列。
交点在两线段交接处，且两线段在射线同侧；

此时虽然射线与线段有交点，但实际并未发生多边形内外区域穿越，这类交点需要排除不纳入统计。
射线起点在线段上；

相交的一种特殊场景，可以按正常相交进行计算。
射线经过平行线；

同第二种场景，由于未发生内外区域穿越，这类交点同样需要排除。

地理坐标特殊场景

对于地理坐标计算时，当出现坐标跨越换日线时，此时再用平面几何的方式进行计算会出现异常。例如计算点坐标（170, 0）水平射线与线段(170, 10), (-170, -10)是否相交：

在地理坐标系里，坐标之间始终按照劣弧方式进行连接，因此(170, 10), (-170, -10)表示的是跨越换日线的线段，而平面坐标系里则不存在这一限制。此时可以对图形坐标进行水平扩展，例如(170, 10), (-170, -10)转换为(170, 10), (190, -10)。

曲线的计算

上面的计算里将几何图形简化为了直线连接，不过实际计算过程中往往会遇到圆弧连接、大圆线连接等场景，此时需要自行计算其最小包围矩形并加入到R树索引，或采用线段拟近的方式进行拆解。

参考资料：

地理坐标系与投影坐标系辨析
投影类型
Wiki: Point in polygon

阮伟聪