判断点是否在多边形内

判断某个点是否在某多边形内是常见的问题，比如某点是否在某个城市内、某城市位于哪个国家等。但地球是一个球面，对于比较小的多边形，采用平面的算法通常不会有问题。但是在某些情况下（尤其是多边形非常大的时候），直接应用平面上的算法会有问题。其中最常见的就是多边形横跨 180° 经线的场景。

本文将以 Python 作为示例代码，探索解决这种问题的方案。

数据来源

使用模拟的数据可能不如使用真实数据更能显示其应用性。这里我选取的之前的文章曾经使用过的 Nuvel 板块模型数据，作为测试的多边形。测试的点选择了位置上有代表性的 9 个城市。这些城市的坐标为：

>>> cities = {
...     'Beijing': (39.9075, 116.3880), 'Dubai': (25.2631, 55.2972),
...     'Honolulu': (22.3089, -157.8261), 'Johannesburg': (-26.1833, 27.9989),
...     'Mumbai': (19.0728, 72.8825), 'New York': (40.7144, -74.0060),
...     'Paris': (48.8583, 2.2944), 'Rio': (-22.9083, -43.1964),
...     'Sydney': (-33.9167, 150.8833)
... }

而板块数据，作者做了适当调整。原因在于其原数据和城市坐标表示西经的方式不同：板块数据使用大于 180 的值表示西经，而城市坐标采用负数表示。因此，我将板块数据中大于 180 的经度减去 360，使其与城市坐标一致。板块和城市的分布如下图。

工具包

对于点与多边形的问题，目前主要使用的有两个 Python 包：matplotlib.path 和 Shapely。使用这两个包都能解决问题，接下来我将分别演示每一种的使用方法。我个人偏好于 Shapely，因为其处理图形关系方面功能更强大一些。

下面分别使用两个包测试阿拉伯板块（Arabian）与各城市的关系。使用 matplotlib.path 的代码如下：

>>> import numpy as np
>>> from matplotlib.path import Path
>>> arabian = np.genfromtxt('standard-plates/Arabian.txt', names=['lon', 'lat'],
...                         dtype=float, comments=':')
...
>>> polygon1 = Path([(lon, lat) for lon, lat in arabian])
>>> for city, coord in cities.items():
...     if polygon1.contains_point([coord[1], coord[0]]):
...         print(f'{city} is in Arabian plate!')

使用 Shapely 的代码如下：

>>> import numpy as np
>>> from shapely.geometry import Point
>>> from shapely.geometry.polygon import Polygon
>>> arabian = np.genfromtxt('standard-plates/Arabian.txt', names=['lon', 'lat'],
...                         dtype=float, comments=':')
...
>>> polygon2 = Polygon([(lon, lat) for lon, lat in arabian])
>>> for city, coord in cities.items():
...     point = Point(coord[1], coord[0])
...     if polygon2.contains(point):
...         print(f'{city} is in Arabian plate!')

两份代码的执行结果都是打印出：

Dubai is in Arabian plate!

这是正确的。但别高兴得太早，如果把上述的阿拉伯板块（Arabian）替换为太平洋板块（Pacific），你就会发现问题了。两份代码都打印出：

Beijing is in Pacific plate!
Dubai is in Pacific plate!
Mumbai is in Pacific plate!
New York is in Pacific plate!
Paris is in Pacific plate!

这是因为，太平洋板块这个多边形横跨了 180 度经线。我们都知道，在大地坐标系中，东西经 180 度位于同一条经线，但这两个包应用的是平面算法。所以程序不认为这个板块的坐标是连续的。解决这个问题的方法有两个：一是应用传统制图学的方法，将多边形和点位投影到平面上，然后使用投影平面上的坐标进行位置关系计算；二是预处理跨 180 度经线的多边形，将其沿 180 度经线分割成小的区块，分别与点位计算位置关系。

有一个名为 pyproj 的包可以处理地图投影问题。但我采用了第二种方案。因为我担心地图投影有可能会将事情变得复杂：选择适合的地图投影和投影参数将是个麻烦。幸好本示例中需要分割的板块还不太多（只有 4 个），它们是：太平洋板块（Pacific）、澳大利亚板块（Australian）、北美洲板块（North American）和南极洲板块（Antarctic）。将它们沿 180 度经线分割后的图形如下。

我已将本文使用的代码、分割前后的板块数据整理为一个压缩文件，你可以下载验证使用。最后使用分割后的边界重新测试：

>>> import json
>>> from shapely.geometry import Point
>>> from shapely.geometry.polygon import Polygon
>>> plates = {}
>>> for name in ['African', 'Antarctic', 'Arabian', 'Australian', 'Caribbean',
...              'Cocos', 'Eurasian', 'Indian', 'Juan', 'Nazca', 'North_Am',
...              'Pacific', 'Philippine', 'Scotia', 'South_Am']:
...     with open(f'new/{name}.json') as j_f:
...         plates[name] = json.load(j_f)
...
>>> boundaries = {}
>>> for name, info in plates.items():
...     boundaries[name] = []
...     for frag in info['fragments']:
...         poly = Polygon(frag)
...         boundaries[name].append(poly)
...
>>> for city, coord in cities.items():
...     point = Point(coord[1], coord[0])
...     for name, bound in boundaries.items():
...         if any(frag.contains(point) for frag in bound):
...             print(f'{city} is in {name}!')

以上代码的输出为：

Beijing is in Eurasian!
Dubai is in Arabian!
Honolulu is in Pacific!
Johannesburg is in African!
Mumbai is in Indian!
New York is in North_Am!
Paris is in Eurasian!
Rio is in South_Am!
Sydney is in Australian!

完美解决问题！