在大型 GeoJSON 对象中查找多边形坐标交点的最有效方法

Question

我正在做一个需要坐标映射的项目 - 确定坐标点是否存在于一系列多边形中。映射的数量非常大——跨越 100 多个多边形的约 1000 万个坐标。

在继续之前，我已经看过问题 here and 。这个问题不是多余的，因为它涉及动态点和静态多边形。

我通过在 200 万个多边形的子集上映射单个坐标来缩小该问题的项目范围。这是我使用的代码：

from shapely.geometry import shape, Point

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 秒（太多时间）。

我也试过使用 mplPath 如下：

import matplotlib.path as mplPath

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 秒。一个改进，但考虑到我需要的映射量仍然有点慢。

有没有更快的方法来实现这样的映射？我不喜欢使用 PySpark 和分布式计算，因为我认为这是一个核选项，但如果需要，我愿意使用它。是否可以矢量化计算而不是遍历多边形？我将生成一个更新，显示使用 numba 是否有任何改进。

Answer 1

我会使用空间连接。

给定这个假数据：

我会用“within”谓词加入它：

from shapely.geometry import Point, Polygon
import geopandas

polys = geopandas.GeoDataFrame({
    "name": ["foo", "bar"],
    "geometry": [
        Polygon([(5, 5), (5, 13), (13, 13), (13, 5)]),
        Polygon([(10, 10), (10, 15), (15, 15), (15, 10)]),
    ]
})

pnts = geopandas.GeoDataFrame({
    "pnt": ["A", "B", "C"],
    "geometry": [
        Point(3, 3), Point(8, 8), Point(11, 11)
    ]
})

result = geopandas.sjoin(pnts, polys, how='left', op='within')

我得到：

pnt                  geometry  index_right name
  A   POINT (3.00000 3.00000)          NaN  NaN
  B   POINT (8.00000 8.00000)          0.0  foo
  C POINT (11.00000 11.00000)          0.0  foo
  C POINT (11.00000 11.00000)          1.0  bar

在大型 GeoJSON 对象中查找多边形坐标交点的最有效方法

Most efficient way to find polygon-coordinate intersections in large GeoJSON object

python

computational-geometry

python-3.x

shapely

geopandas