时间戳

首页 > 解决方案 > 在大型 GeoJSON 对象中查找多边形坐标交点的最有效方法

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

问题描述

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

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

我通过在 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有任何改进。

标签： pythonpython-3.xcomputational-geometrygeopandasshapely

解决方案

我会使用空间连接。

鉴于此假数据：

在此处输入图像描述

我将它与“内部”谓词一起加入：

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

推荐阅读