【数据技术】利用Python判断坐标点是否在多边形内

本文讲述利用Python判断坐标点是否在多边形内的思路及代码，并以厦门AOI数据为例，测试代码的功能及性能。

本文完整代码的获取方式在文末，有需要的小伙伴自取。

利用Python中的GeoPandas功能包，能够实现判断坐标点是否在多边形，具体代码如下

下图为代码的运行结果

testPoint1被判定在多边形内，testPoint2被判定不在多边形内

下图为代码的地图可视化

红色线框为testPolygon，蓝色三角为testPoint1，绿色三角为testPoint2

利用厦门AOI数据，对代码的功能及性能进行测试，下图为厦门AOI数据

厦门AOI数据，共计5542个多边形

以厦门岛为边界每隔50米创建一个测试坐标点，具体代码如下

# 创建测试坐标点startLon = 118.062706# 起始经度endLon = 118.213566# 结束经度startLat = 24.420724# 起始纬度endLat = 24.565696# 结束纬度lonLength = 0.0005# 经度间隔latLength = 0.0005# 纬度间隔testPointList = []  # 存放测试坐标点的列表for i in range(0, 9999999999):    currentLon = startLon + lonLength * i    if currentLon > endLon:break    for j in range(0, 9999999999):        currentLat = startLat + latLength * j        if currentLat > endLat:break        testPointList.append(Point(currentLon, currentLat))print("测试坐标点数量：", len(testPointList))

下图为测试坐标点数据集

测试点数据集，共87580个坐标点

针对每个坐标点，遍历全部5542个AOI多边形，判断当前坐标点是否在多边形内，如果在则跳出循环，如果不在则继续

testPolygonList = []  # 存放多边形的列表# 利用geopandas读取geojson文件，也支持读取shp文件loadData_XM = gp.read_file("D:\\厦门AOI.geojson") for i inrange(0, len(loadData_XM)):    testPolygonList.append(loadData_XM.loc[i, "geometry"])print("多边形数量：", len(testPolygonList))time1 = datetime.now()  # 记录计算开始时间点for i inrange(0, len(testPointList)):    currentPoint = testPointList[i]    judgeResult = "false"# 记录坐标点是否在多边形内的变量，true为在，false为不在for j inrange(0, len(testPolygonList)):if testPolygonList[j].contains(currentPoint):  # 利用contains()可以判断多边形和点的关系            judgeResult = "true"breaktime2 = datetime.now()  # 记录计算结束时间点print("耗时：", (time2 - time1).seconds, "秒")

下图为代码运行结果，总共耗时943秒

下图为代码运行结果的地图可视化

红点为被判定在多边形内的点，黑点为被判定不在多边形内的点

上述代码已实现利用Python判断坐标点是否在多边形内的功能，但如果遇到大规模点数据或者大规模的多边形则需要等待较长时间，针对这种场景对代码进行优化。

优化的思路是在执行判断坐标点是否在多边形内的空间计算之前，先判断坐标点是否在多边形的外接矩形内，如果坐标点不在多边形的外接矩形内则不可能在多边形内，这样能减少空间计算的次数。

红框为多边形，紫框为多边形的外接矩形，绿色三角不在外接矩形内则不需要进行空间计算

下面为执行上述操作的代码

# 读取厦门AOI数据并建立外接矩形testPolygonList = []  # 存放多边形的列表testRectangleList = []  # 存放外接矩形的列表loadData_XM = gp.read_file("D:\\厦门AOI.geojson")for i inrange(0, len(loadData_XM)):    testPolygonList.append(loadData_XM.loc[i, "geometry"])    minLon = loadData_XM.loc[i, "geometry"].envelope.boundary.coords[0][0]  # 当前多边形的最小经度    minLat = loadData_XM.loc[i, "geometry"].envelope.boundary.coords[0][1]  # 当前多边形的最大经度    maxLon = loadData_XM.loc[i, "geometry"].envelope.boundary.coords[2][0]  # 当前多边形的最小纬度    maxLat = loadData_XM.loc[i, "geometry"].envelope.boundary.coords[2][1]  # 当前多边形的最大纬度    testRectangleList.append([minLon, minLat, maxLon, maxLat])print("多边形数量：", len(testPolygonList))time1 = datetime.now()  # 记录计算开始时间点for i inrange(0, len(testPointList)):    currentPoint = testPointList[i]    currentLon = currentPoint.x    currentLat = currentPoint.y    judgeResult = "false"for j inrange(0, len(testPolygonList)):        # 先判断坐标点是否在外接矩形内，如果在则进行空间计算if currentLon >= testRectangleList[j][0] and currentLon <= testRectangleList[j][2] and currentLat >= \                testRectangleList[j][1] and currentLat <= testRectangleList[j][3]:if testPolygonList[j].contains(currentPoint):                judgeResult = "true"breaktime2 = datetime.now()  # 记录计算结束时间点print("耗时：", (time2 - time1).seconds, "秒")

基于厦门AOI数据测试优化后代码的性能，结果耗时70秒，较之前的方法提升了10倍性能

为方便初学者使用，将上述方法及代码进行整理（获取方式见文末），目前仅需替换坐标点文件及多边形文件就可以实现判断坐标点是否在多边形内的功能。

替换坐标点文件及多边形文件

输出判断结果，TRUE为在多边形内，FALSE为不在多边形内

同时，也支持计算每个多边形内的坐标点数量

多边形上的数字为在其范围内的坐标点数量