我正在做一个坐标映射相关的项目，核心需求是判断坐标点是否属于指定的多边形集合。目前项目的数据规模相当大：要处理约1000万个坐标点，匹配超过1亿个多边形。我已经查过相关问题，这次不是重复提问——我的场景是动态坐标点与静态多边形的匹配，和之前看到的批量静态点匹配场景不太一样。

为了聚焦问题，我先缩小了测试范围：用单个坐标点匹配200万个建筑轮廓多边形，分别试了两种方法：

方法1：基于Shapely的实现 #

from shapely.geometry import shape, Point
import json

f = open('path/to/file.geojson', 'r')
data = json.loads(f.read())
point = Point(42.3847, -71.127411)
for feature in data['features']:
    polygon = shape(feature['geometry'])
    if polygon.contains(point):
        print(polygon)

遍历200万个多边形花了大概30秒，这个效率完全达不到项目要求。

方法2：基于mplPath的实现 #

import matplotlib.path as mplPath
import numpy as np
from tqdm import tqdm
import json

f = open('path/to/file.geojson', 'r')
data = json.loads(f.read())

building_arrays = [np.array(data['features'][i]['geometry']['coordinates'][0]) for i, v in enumerate(tqdm(data['features']))]
bbPath_list = [mplPath.Path(building) for building in tqdm(building_arrays)]
for b in tqdm(bbPath_list):
    if b.contains_point((-71.1273842, 42.3847423)):
        print(b)

这个方法耗时约6秒，比Shapely快了不少，但放到项目整体的千万级坐标点+亿级多边形的规模里看，还是慢得没法接受。

现在想请教各位大佬：

• 有没有更高效的实现方式？目前我暂时不想用PySpark这类分布式计算方案（打算把它当作最后的备选），但实在不行的话也可以考虑。
• 能不能用向量化计算替代逐一遍历多边形的方式来提升效率？后续我会更新用numba优化后的测试结果。

内容的提问来源于stack exchange，提问作者mmz