如何加速 geopandas 空间连接？

Question

我有两个 geopandas 数据框。对于左框架中的每一行，我想找到右框架中的哪些行在空间上与该行重叠以及重叠程度。一旦获得这些信息，我就能够根据重叠程度进行空间连接。

不幸的是，我在大量多边形上执行此操作：一个州的所有美国人口普查区（德克萨斯州有 5,265 个）和大量大小相似（但不重合）的美国人口普查区块（德克萨斯州有其中约 914,231 个）。

我正在寻找一种方法来更快地做到这一点。到目前为止我的代码如下。

使用的数据集可以从美国人口普查中获取：街区数据、区域数据。

#!/usr/bin/env python3

import geopandas as gpd
import geopandas_fast_sjoin as gpfsj
import time
import os
import pickle
import sys

os.environ["GDAL_DATA"] = "/usr/share/gdal"

TRACT_FILE    = "./data/tracts/tl_2010_{fips}_tract10.shp"
BLOCK_FILE    = "./data/blocks/tabblock2010_{fips}_pophu.shp"
PROJECTION    = '+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs'

print("Reading data...")
start_time = time.time()
tracts = gpd.read_file(TRACT_FILE.format(fips=48))
blocks = gpd.read_file(BLOCK_FILE.format(fips=48))
print('Time: ', time.time()-start_time )

print("Converting coordinate reference systems...")
start_time = time.time()
tracts = tracts.to_crs(PROJECTION)
blocks = blocks.to_crs(PROJECTION)
print('Time: ', time.time()-start_time )

print("Performing spatial join...")
start_time = time.time()
joined = gpd.sjoin(tracts, blocks, how='left')
print('Time: ', time.time()-start_time )

print("Calculating areas of intersection...")
start_time = time.time()
joined['area_of_intersect'] = [row['geometry'].intersection(blocks.loc[row['index_right']]['geometry']).area for i,row in joined.iterrows()]
print('Time: ', time.time()-start_time )

Answer 1

有几种优化可以使此操作更快：在 C++ 中完成所有工作而不涉及 Python，使用空间索引快速识别交叉点的候选对象，使用准备好的几何图形快速检查候选对象，并在可用的情况下并行化整个操作核心。

所有这些都可以在 Python 中完成，但会带来一些开销，但在我的测试中，在数 GB 数据集上使用多处理模块会导致 Python 耗尽可用内存。在 Windows 上，这可能是不可避免的，在 Linux 上，写时复制应该可以防止这种情况发生。也许可以通过仔细的编程来完成，但使用 Python 的全部意义在于不必担心这些细节。因此，我选择将计算转移到C++。

为了实现这一目标，我使用pybind11构建了一个新的 Python 模块，该模块接受来自 geopandas 的几何图形列表并生成三个列表的输出：(1) 左侧几何图形的行索引；(2) 右手几何的行索引；(3) 两者之间重叠的面积（仅当 >0 时）。

例如，对于几何图形 left=[A,B,C,D] 和 right=[E,F,G,H] 的输入，令：

• E 完全在 A 中
• F 与 A 和 B 重叠
• G 和 H 没有重叠

然后返回看起来像：

List1 List2 List3
A     E     Area(E)
A     F     AreaIntersection(A,F)
B     F     AreaIntersection(B,F)

在我的机器上，该sjoin操作花费了 73 秒，计算交叉点花费了 1,066 秒，总共 1139 秒（19 分钟）

在我的 12 核机器上，下面的代码需要 50 秒才能完成所有这些工作。

因此，对于需要交叉区域的空间连接，这仅节省了一点时间。但是，对于需要交叉区域的空间连接，这可以节省大量时间。换句话说，计算所有这些交叉点需要大量工作！

在进一步的测试中，在不使用准备好的加速几何图形的情况下计算相交面积在 12 个核心上花费了 287 秒。因此，并行化交叉点可实现 4 倍的加速，而并行化准备好的几何图形可实现 23 倍的加速。

生成文件

all:
    $(CXX) -O3 -g -shared -std=c++11 -I include `python3-config --cflags --ldflags --libs` quick_join.cpp -o geopandas_fast_sjoin.so -fPIC -Wall -lgeos -fopenmp 

快速加入.cpp

List1 List2 List3
A     E     Area(E)
A     F     AreaIntersection(A,F)
B     F     AreaIntersection(B,F)

测试.py

all:
    $(CXX) -O3 -g -shared -std=c++11 -I include `python3-config --cflags --ldflags --libs` quick_join.cpp -o geopandas_fast_sjoin.so -fPIC -Wall -lgeos -fopenmp