首页 > 解决方案 > 通过循环将多边形内的所有点作为值添加到字典对象

问题描述

我有一个包含 13 个多边形的 shapefile,以及来自与多边形相同的一般区域的约 33K 地址的地址列表。我通过谷歌地理编码 API 对地址进行了地理编码,现在我试图查看哪些地址在哪些多边形中,并按多边形的名称对它们进行分组。

我可以一次处理一个多边形,但我在循环中遗漏了一些东西。

到目前为止,这是我所在的位置:

# Import shapefile and convert coordinates to match address file
sf = gpd.read_file(MY_SHAPEFILE)
sf_geo = sf.to_crs(epsg=4326)

# Import geocoded addresses
address_data = pd.read_csv(ADDRESS_FILE)
# Create points from lat/lon coordinate columns in address file
geometry_points = [Point(xy) for xy in zip(address_data['longitude'], 
                   address_data['latitude'])]

# Create object from one of the Polygons
p = sf_geo.iloc[0].geometry

i = 0
for point in geometry_points:
    if point.within(p):
        i += 1
        print(i)
    else:
        continue   

以上在所有多边形上都可以正常工作。但是,我真正希望的是能够更新一个字典,其中键是多边形的实际名称,值是该多边形内匹配的所有点。然后,我可以将多边形名称添加到地址列表中。

# Extract the names of each polygon
area_names = list(sf_geo['Name'])
# Create dict of polygon name : geometry
for r in sf_geo:
    shape_dict = dict(zip(area_names, sf['geometry']))

# Initialize empty dict to hold list of all addresses within each polygon
polygon_contains_dict = {k: [] for k in area_names}

以上在打印时创建了这种格式的字典:

{'10 ppm': <shapely.geometry.polygon.Polygon object at 0x7fea194225d0>, '20 ppm': <shapely.geometry.polygon.Polygon object at 0x7fe9f2e23590>, ETC}

以及一个字典,其中键与 shape_dict 相同,但值是空列表。

我正在使用以下内容尝试遍历 shape_dict 中的所有键以及从地址创建的所有点,并更新一个列表,然后该列表将成为 polygon_contains_dict 中每个键的值:

for key, value in shape_dict.items():
    contains_list = []
    not_contained = []
    for point in geometry_points:
        if point.within(value):
            contains_list.append(point)
        else:
            not_contained.append(point)
    polygon_contains_dict[key] = contains_list

但是,这对 contains_list 和(显然)polygon_contains_dict 中的值都没有任何添加。所有点都被倾倒在 not_contained 中。

因为我知道点实际上在某些多边形内,所以我知道我遗漏了一些东西。geometry_points 中的所有点都是 Point 对象,而 shape_dict.values 中的所有多边形都是 Polygon 对象。

我错过了什么?谢谢你的帮助。

标签: pythongeopandasshapely

解决方案


我建议您避免完全循环并为您的坐标和地址数据创建第二个 geopandas 数据框,然后进行空间连接:

# Import geocoded addresses
address_data = pd.read_csv(ADDRESS_FILE)
# Create points from lat/lon coordinate columns in address file
geometry_points = [Point(xy) for xy in zip(address_data['longitude'], 
               address_data['latitude'])]
address_gpd=gpd.GeoDataFrame(address_data,crs={'init': 'epsg:4326'},geometry=geometry_points) # second geopandas frame

# inner spatial join with shapefile 
df=gpd.sjoin(sf_geo,address_gpd,how='inner',op='intersects') 

数据框现在将df具有每个多边形内的所有地址,然后“更新一个字典,其中键是多边形的实际名称,值是该多边形内匹配的所有点”,您可以使用 groupby 和 to_dict

df=df.groupby(['Name'])['address'].apply(list)
polygon_contains_dict=df.to_dict('index') 

我假设你的地址的列名是address这样的,如果不是这样的话。

有关空间连接的更多信息,请参阅有关合并数据的 geopandas 文档。


推荐阅读