geopandas进行空间查询速度优化的几种方法（R树等）

背景（请跳过）：

带有地理坐标的大数据为我们研究人类移动模式提供了新的研究契机。但是在处理大量地理数据的过程中，传统的地理信息数据处理工具，如Arcgis，Qgis等，已经不能满足人们对于海量数据的处理需求。以arcpy为例，其封闭的环境、较低的处理效率以及复杂的函数组织结构让大量小白们望而却步。Geopandas为熟悉python pandas包的人们进行地理数据处理提供了较为低成本的学习方式。其官方文档为：geopandas

空间查询优化方法：

采用geopandas包进行简单处理非常方便。但是对大量数据进行处理时，处理效率较低。因而人们提供了各种方法对地理空间数据查询（如包含查询，相交查询等）进行了优化。下面我们实验几种优化方法：

1、R树：

R树的基础知识链接：百度百科、维基百科

基本代码：

spatial_index = gdf.sindex     
#建立R树索引. gdf,geopandas dataframe 数据
possible_matches_index = list(spatial_index.intersection(polygon.bounds)) 
possible_matches = gdf.iloc[possible_matches_index]   
#通过R树进行简单包含查询
precise_matches = possible_matches[possible_matches.intersects(polygon)]    
#将简单查询结果进行进行匹配，得到精确查询结果

该代码参考网址。同时文章中包含几何面较复杂的时候一些简单的处理方法。

测试：

我们测试为：统计北京五环之内各个网格内上车的数量，其中网格大小为4km*4km, 网格个数为64，出租车数目为204700，未加空间索引的代码如下：

gridgdf=gpd.read_file(grid_path)
taxi_data=pd.read_csv(taxi_path)
taxi_data=taxi_data.set_geometry('O_point')
#数据处理，注意taxi_data和gridgdf的坐标系要一致

for index,temp_grid in gridgdf.iterrows():
    taxi_data.loc[taxi_data.within(temp_grid['geometry']),'O_id']=temp_grid['id']
#查询

通过R树建立空间索引并查询的代码如下：

gridgdf=gpd.read_file(grid_path)
taxi_data=pd.read_csv(taxi_path)
taxi_data=taxi_data.set_geometry('O_point')
spatial_index = taxi_data.sindex

for index,temp_grid in gridgdf.iterrows():
    possible_matches_index=list(spatial_index.intersection(temp_grid.geometry.bounds))
    possible_matches=taxi_data.iloc[possible_matches_index]
    prices_mathces_index=possible_matches.intersects(temp_grid.geometry).index
    taxi_data.loc[prices_mathces_index,'O_id']=temp_grid['id']

未建立空间索引测试结果为62s，建立空间索引的测试查询时间为5s。由于采用的是正方形网格的查询，因而速度提升不明显。当采用复杂几何查询时，采用空间索引得到的优化速度将会更好。

2、矢量化：

shaply包的vectorized包含着一些对查询的优化，通过shaply和geopandas一起协作可以达到较好的优化效果。我们同样以空间点查询进行说明，代码来源

import numpy as np, shapely.vectorized as sv
from shapely.geometry import Point, Polygon
polygon = Polygon([(10,10), (10,100), (100,100), (100, 10)])
xy = np.array([(n, n) for n in range(10000)])
points = [Point(x, y) for x, y in xy]
#初始化数据

def contains_py(polygon, points):
    return np.array([polygon.contains(point) for point in points])
%timeit contains_py(polygon, points)
#传统查询方法——遍历所有的点，并比较

%timeit sv.contains(polygon, x=xy[:,0], y=xy[:,1])
#通过shapely.vectorized方法进行比较的结果

作者在原博客中说，传统方法的查询时间为40ms/次，而shapely.vectorized方法为3ms/次。相较于传统方法该方法提升明显。我们同样在上文我们所说的数据中进行了实验，实验代码如下：

def is_within(df,polygon,xlabel,ylabel):
    xlist=df[xlabel].values
    ylist=df[ylabel].values
    is_list=sv.within(polygon,xlist,ylist)
    index=df.index
    series=pd.Series(data=is_list,index=index)
    return series

gridgdf=gpd.read_file(grid_path)
taxi_data=pd.read_csv(taxi_path)
taxi_data=taxi_data.set_geometry('O_point')

for index,temp_grid in gridgdf.iterrows():
    is_within_series=is_within(df=taxi_data,polygon=temp_grid.geometry,xlabel='O_lon',ylabel='O_lat')
    taxi_data.loc[is_within_series, 'O_id'] = temp_grid['id']

实验结果查询时间为5s，与采用R树空间索引时间相近。该方法计算较快，但是shapely.vectorized中只能判断contain和touch两种空间关系，可以判断的类型较少。

进一步测试与建议

我们采用同样数据，但是将空间查询网格大小逐渐缩小，分别为2km*2km，1km*1km，500m*500m计算空间索引方法与矢量化方法的计算时间变化。得到的结果如下：

栅格大小	矢量方法	R树
2km	20s	5s
1km	28s	13s
500m	340s	17s

通过对比我们发现，随着空间网格的逐渐变小，R树的查询速率要快于矢量方法。因此对于两种方法的比较后，我们可以总结如下：
- 当数据量较小情况下:如果采用的数据不是dataframe格式，采用矢量方法运算速度较快，且较方便。
- 当数据量较大的情况下：采用geopandas+R树的方法查询较快。