Python结合shapely将GeoJSON多面数据存入DataFrame并判断点面包含关系

问题原因 #

该Shapely弃用警告并非shape()方法本身导致，核心触发逻辑是：pandas构造DataFrame时会自动迭代传入的可迭代对象，而MultiPolygon等多部件几何属于可迭代类型，pandas尝试将其拆分为多行存储的迭代动作，触发了shapely的弃用提示。
直接存储坐标再构造几何报错的原因是：GeoJSON中Polygon和MultiPolygon的coordinates结构存在差异，手动构造需要先判断几何类型，处理复杂度高，远不如直接用shape()解析稳妥，只要避免pandas迭代几何对象即可解决警告问题。

无警告修改方案 #

只需要在构造每行DataFrame时，把name和area的值都套一层列表，让pandas识别为单个单元格值，不会触发迭代动作：

import pandas as pd
from shapely.geometry import shape

# 修改后的DataFrame构造代码
choro = pd.concat([
    pd.DataFrame({
        'name': [block['properties']['name']],
        'area': [shape(block['geometry'])]
    }) 
    for block in geojson['features']
], ignore_index=True)

点面判断代码 #

原有判断逻辑不需要修改，area列存储的是可直接调用shapely方法的Polygon/MultiPolygon对象：

for index, area in choro.iterrows():
    polygon = area['area']
    if polygon.contains(point):
        print('Point belongs to ', area['name'])

更高效的可选方案（推荐空间场景使用） #

如果经常处理空间数据，可以直接使用geopandas构造GeoDataFrame，代码更简洁，同时内置了批量空间运算能力，不需要循环判断：

import geopandas as gpd

# 直接从GeoJSON构造GeoDataFrame
choro = gpd.GeoDataFrame.from_features(geojson['features'])
# 批量判断点所属区域，直接返回匹配的行
matched_area = choro[choro.contains(point)]
if not matched_area.empty:
    print('Point belongs to ', matched_area.iloc[0]['name'])

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