如何在Python中基于Shapely实现更快的点面相交检测？

我正在用Python复现GIS软件的功能——将51个质心点与9万+地块多边形匹配，检测点是否位于多边形内部。目前基于Shapely编写的脚本可以正常运行，但耗时约5分钟。考虑到后续需要处理更多的质心点和地块数据，希望找到优化速度的方法。

原代码 #

import json, csv, time
from shapely.geometry import Point, Polygon

start = time.time()

with open('M035TaxPar_CY20_FY20.csv') as csv_file:
    csv_reader = csv.reader(csv_file, delimiter=',')
    line_count = 0
    for row in csv_reader:
        if line_count == 0:
            #print(f'Column names are {", ".join(row)}')
            line_count += 1
        else:
            my_coord = tuple([float(row[10]), float(row[11])])
            my_coord = Point(my_coord)
            line_count += 1
            
            with open('boston_parcels_all.json', 'r') as f:
                array = json.load(f)
                
            dict_size = len(array.get("features"))
            for i in range(0, dict_size -1):
                sub_dict = array.get("features")[i] 
                geometry_dict = sub_dict.get("geometry")
                current_Pol = geometry_dict.get("coordinates")
                
                attribute_data = sub_dict.get("properties")
                parcel_id = attribute_data.get("MAP_PAR_ID")
                
                try:
                    final_Pol = Polygon(current_Pol[0])
                except AttributeError:
                    continue
                except ValueError:
                    continue
                if not my_coord.within(final_Pol):
                    continue
                else:
                    if parcel_id == None:
                        continue
                    else:
                        print(f"My pol: {parcel_id}, my point: {row[0]}")
                
    print(f'Processed {line_count} lines.')
end = time.time()
total_time = end - start
print("\n"+ str(total_time))

数据示例 #

CSV文件（M035TaxPar_CY20_FY20.csv） #

MAP_PAR_ID,LOC_ID,POLY_TYPE,MAP_NO,SOURCE,PLAN_ID,LAST_EDIT,BND_CHK,NO_MATCH,TOWN_ID,X_coord,Y_coord
0301290001,F_772282_2959076,FEE,3,ASSESS,,20191107,,N,35,-71.07030535235137,42.36703589444949
0301286000,F_772818_2959719,FEE,3,ASSESS,,20191107,,N,35,-71.0682938996932,42.36876187375208

JSON文件（boston_parcels_all.json） #

{
    "type" : "FeatureCollection",
    "name" : "M035TaxPar_CY20_FY20",
    "features" : [
        {
            "type" : "Feature",
            "geometry" : {
                "type" : "Polygon",
                "coordinates" : [
                    [
                        [ -71.074767681, 42.3489358757 ],
                        [ -71.0745071763, 42.3490067112 ],
                        [ -71.074852421, 42.3497214512 ],
                        [ -71.0757149372, 42.3494869825 ],
                        [ -71.0753495281, 42.3487373608 ],
                        [ -71.074767681, 42.3489358757 ]
                    ]
                ]
            },
            "properties" : {
                "SHAPE_Leng" : 324.191190544,
                "SHAPE_Area" : 6483.24124923,
                "MAP_PAR_ID" : "0401134000",
                "LOC_ID" : "F_771022_2952578",
                "POLY_TYPE" : "FEE",
                "MAP_NO" : "4",
                "SOURCE" : "ASSESS",
                "LAST_EDIT" : 20191107,
                "NO_MATCH" : "N",
                "TOWN_ID" : 35
            }
        },
        {
            "type" : "Feature",
            "geometry" : {
                "type" : "Polygon",
                "coordinates" : [
                    [
                        [ -71.0554563931, 42.3547109221 ],
                        [ -71.0550339868, 42.3548837812 ],
                        [ -71.0557334699, 42.3555313301 ],
                        [ -71.0564003565, 42.3554356917 ],
                        [ -71.0562325294, 42.3545663782 ],
                        [ -71.0554563931, 42.3547109221 ]
                    ]
                ]
            },
            "properties" : {
                "SHAPE_Leng" : 351.344209198,
                "SHAPE_Area" : 7627.60108948,
                "MAP_PAR_ID" : "0304410000",
                "LOC_ID" : "F_776224_2954721",
                "POLY_TYPE" : "FEE",
                "MAP_NO" : "3",
                "SOURCE" : "ASSESS",
                "LAST_EDIT" : 20191107,
                "NO_MATCH" : "N",
                "TOWN_ID" : 35
            }
        }
    ]
}

优化方案及代码 #

核心优化点 #

• 避免重复IO操作：原代码每个点都重新加载JSON文件，改为仅加载一次。
• 预处理多边形数据：提前创建所有有效Polygon对象并关联对应ID，避免循环中重复创建。
• 空间索引过滤：用R-tree索引快速筛选候选多边形，仅对可能包含点的多边形做精确判断，大幅降低计算量。

优化后代码（Shapely+Rtree） #

先安装依赖：pip install rtree shapely

import json, csv, time
from shapely.geometry import Point, Polygon
from rtree import index

start = time.time()

# 1. 预处理多边形与空间索引
polygons = []
idx = index.Index()

with open('boston_parcels_all.json', 'r') as f:
    feature_collection = json.load(f)
    for i, feature in enumerate(feature_collection.get("features")):
        geometry = feature.get("geometry")
        if geometry.get("type") != "Polygon":
            continue
        coords = geometry.get("coordinates")
        try:
            poly = Polygon(coords[0])
        except (AttributeError, ValueError):
            continue
        parcel_id = feature.get("properties").get("MAP_PAR_ID")
        if not parcel_id:
            continue
        polygons.append((poly, parcel_id))
        idx.insert(i, poly.bounds)

# 2. 处理CSV点数据
with open('M035TaxPar_CY20_FY20.csv') as csv_file:
    csv_reader = csv.reader(csv_file, delimiter=',')
    next(csv_reader)  # 跳过表头
    line_count = 0
    for row in csv_reader:
        line_count += 1
        x = float(row[10])
        y = float(row[11])
        point = Point(x, y)
        
        # 用索引筛选候选多边形
        candidate_indices = list(idx.intersection(point.bounds))
        for i in candidate_indices:
            poly, parcel_id = polygons[i]
            if point.within(poly):
                print(f"My pol: {parcel_id}, my point: {row[0]}")
                break  # 单匹配场景下找到即停止，多匹配则移除break

print(f'Processed {line_count} lines.')
end = time.time()
print(f"\nTotal time: {total_time:.2f} seconds")

更简洁高效方案：GeoPandas #

GeoPandas内置空间索引，适合大规模空间数据处理，代码更简洁。
安装依赖：pip install geopandas

import geopandas as gpd
import time

start = time.time()

# 加载多边形与点数据
gdf_parcels = gpd.read_file('boston_parcels_all.json')
gdf_points = gpd.read_file(
    'M035TaxPar_CY20_FY20.csv',
    GEOM_POSSIBLE_NAMES=['X_coord', 'Y_coord'],
    crs=gdf_parcels.crs
)
gdf_points['geometry'] = gpd.points_from_xy(gdf_points['X_coord'], gdf_points['Y_coord'])

# 空间连接匹配点与多边形
result = gpd.sjoin(gdf_points, gdf_parcels, how='left', predicate='within')

# 输出匹配结果
for _, row in result.dropna(subset=['MAP_PAR_ID_right']).iterrows():
    print(f"My pol: {row['MAP_PAR_ID_right']}, my point: {row['MAP_PAR_ID_left']}")

print(f'Processed {len(gdf_points)} lines.')
end = time.time()
print(f"\nTotal time: {end - start:.2f} seconds")

优化效果说明 #

原代码的核心性能瓶颈是重复加载JSON和无索引遍历所有多边形，优化后：

• 空间索引将单一点的匹配复杂度从O(N)降至O(logN)，9万+多边形场景下速度提升几十倍甚至上百倍。
• GeoPandas方案封装了底层优化，代码更简洁，后续扩展大规模数据更方便。

内容的提问来源于stack exchange，提问作者Robert