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

Bertha 。 2022-10-08 05:55 228阅读 0赞

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

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

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

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

## 一、地理坐标系统 ##

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

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

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

1.  圆锥投影  
    ![圆锥投影][20210620155144271.gif_pic_center]
2.  圆柱投影  
    ![圆柱投影][20210620155207895.gif_pic_center]
3.  平面投影  
    ![平面投影][20210620155231480.gif_pic_center]
4.  极方位投影  
    ![极方位投影][20210620155252927.gif_pic_center]

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

## 二、场景分解 ##

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

**简单地理空间计算**

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

![简单关系计算][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70]

**复杂计算场景**

考虑另一个场景，飞机航行过程中，如何判断当前处于哪些空域范围内。国际上将各个地理空间范围规划为不同的空域，并赋予不同的管理级别，例如我国民航可用空域面积319.53万平方公里，占比32%，其余为军用或受控制空域。实时掌握飞机当前所处空域，能够帮助规避潜在风险。与上面例子不同的是，这里我们无法事先知道有哪些空域需要参与计算，并且空域数量庞大、形状不规则（如下图所示），如果使用轮询的方式套搬上面的方法，会导致计算量过大，无法保证实时性。  
![广州周边空域概况][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70_pic_center]  
这里我们对问题进行拆解，第一步对空域数据进行检索，排除非相关部分；第二步对检索得到的空域进行坐标计算，判断是否为包含关系。

## 三、地理数据检索 ##

通常使用R树索引对地理数据进行检索，这是一种非精确的索引，利用图形最小包围盒MBR进行快速过滤。  
![R树索引原理][R]  
由于R树索引的特性，构建索引时图形粒度越小，检索的效率就会越高。如果将空域作为一个整体构建索引，在检索时会返回较多不相关数据，导致过滤率较低。因此我们转换另一个思路，对构成空域的各个线段进行索引，提升过滤率。  
![索引方式调整][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70 1]

## 四、射线法计算空间关系 ##

索引细化后，我们无法直接判断点与原空域的空间关系，这里我们改用几何光学里的[射线法][Link 1]进行空间计算。  
射线法的原理是：从目标点向无穷远处发射一条直线，如果直线与多边形的交点个数为奇数，则该点位于多边形的内部。  
![射线法原理][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70 2]  
结合R树检索与射线法计算过程后，整体算法调整如下：

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

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

## 五、特殊处理 ##

**射线法特殊场景**

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

1.  射线与线段交点刚好为两线段的交接处；  
    ![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70 3]  
    此时射线同时穿过两线段，但实际交点只有一个。因此在进行相交判断时，需要按左开右闭图形进行计算，同时构成图形的线段数据需要按照有向图排列。
2.  交点在两线段交接处，且两线段在射线同侧；  
    ![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70 4]  
    此时虽然射线与线段有交点，但实际并未发生多边形内外区域穿越，这类交点需要排除不纳入统计。
3.  射线起点在线段上；  
    ![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70 5]  
    相交的一种特殊场景，可以按正常相交进行计算。
4.  射线经过平行线；  
    ![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70 6]  
    同第二种场景，由于未发生内外区域穿越，这类交点同样需要排除。

**地理坐标特殊场景**

对于地理坐标计算时，当出现坐标跨越换日线时，此时再用平面几何的方式进行计算会出现异常。例如计算点坐标（170, 0）水平射线与线段(170, 10), (-170, -10)是否相交：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70 7]  
在地理坐标系里，坐标之间始终按照劣弧方式进行连接，因此(170, 10), (-170, -10)表示的是跨越换日线的线段，而平面坐标系里则不存在这一限制。此时可以对图形坐标进行水平扩展，例如(170, 10), (-170, -10)转换为(170, 10), (190, -10)。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70 8]  
**曲线的计算**

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

## 参考资料： ##

[地理坐标系与投影坐标系辨析][Link 2]  
[投影类型][Link 3]  
[Wiki: Point in polygon][Wiki_ Point in polygon]

阮伟聪

[MySQL]: https://blog.csdn.net/vipshop_fin_dev/article/details/109151777
[20210620155144271.gif_pic_center]: /images/20221005/4158285e4ba24839851de1e7c2f9dffe.png
[20210620155207895.gif_pic_center]: /images/20221005/7a1821e7bb02424e8a84b1592e6f840a.png
[20210620155231480.gif_pic_center]: /images/20221005/484e9fd8607345a0b475dacf5461fa22.png
[20210620155252927.gif_pic_center]: /images/20221005/98192f3494de4f8486a8d36614836466.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70]: /images/20221005/ff303ecf610541c1848bfdc9ade4119a.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70_pic_center]: /images/20221005/fcf4d04e7c184af383a7f253b28d47ae.png
[R]: /images/20221005/ab14fde1c6c1412caff29cc133854038.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70 1]: /images/20221005/f045c62da77148bda5e93fb3cf3e0a1c.png
[Link 1]: https://rosettacode.org/wiki/Ray-casting_algorithm
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70 2]: /images/20221005/4220196e3ffa41538db21b800b7f6bed.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70 3]: /images/20221005/3365ce5510e14834914fdd9a5d84d462.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70 4]: /images/20221005/f88c2dc7ac524f23ab27ccba9971fda7.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70 5]: /images/20221005/aacf4d88b8694c5d94f4522042898bd6.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70 6]: /images/20221005/029157b8e4c3450495520faa6ad262ed.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70 7]: /images/20221005/79f6537e6af840e09ef50e93c61263e3.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L20wXzQ3MDIyNzc2_size_16_color_FFFFFF_t_70 8]: /images/20221005/2041681e6c3146479374aab721fa31a7.png
[20210620210516572.png]: /images/20221005/49d7ae5a3084476e96d1fc959fa590ac.png
[Link 2]: https://blog.csdn.net/wxqqq666/article/details/91542325
[Link 3]: https://desktop.arcgis.com/zh-cn/arcmap/10.3/guide-books/map-projections/projection-types.htm
[Wiki_ Point in polygon]: https://en.wikipedia.org/wiki/Point_in_polygon#Ray_casting_algorithm