大型GeoJSON物体中求多边形坐标交点的最有效方法

Question

通过在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)

Answer 1

我会使用空间连接。

鉴于这一虚假数据：

我会用“inner”谓词来连接它：

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