通过循环将多边形内的所有点作为值添加到字典对象

Question

我有一个包含 13 个多边形的 shapefile，以及来自与多边形相同的一般区域的 ~33K 地址的地址列表。我通过 Google 地理编码 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对象。

我错过了什么？感谢您的帮助。

Answer 1

我建议您避免完全循环并为您的坐标和地址数据创建第二个 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 数据框现在将包含每个多边形内的所有地址然后 "update a dict where the key is the actual name of the polygon, and the values are all of the points that match within that polygon" 您可以使用 groupby 和 to_dict

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

我假设您的地址的列名称是 address，所以如果不是这样，请更改。

有关空间连接的详细信息，请参阅 geopandas docs on merging data。

通过循环将多边形内的所有点作为值添加到字典对象

Adding all points within polygon as value to dictionary object via loop

python

shapely

geopandas