如何在Python中查找多边形或多多边形中包含哪些点？

Question

我有 2 个表，第一个称为“filtered_commerce”，第二个称为“filtered_secteur”。 “filtered_commerce”表包含一个“zone_population”列，其中包含表示地理区域的多边形和多边形。 Filtered_sector 表包含一个包含点的“centroid_degres”列和一个包含点的标识符的“INS9”列。我想找到每个区域中有哪些点，并将每个区域中包含的点的 INS9 保存在“filtered_commerce”的新列中。区域是多边形或多多边形，保存如下：POLYGON ((4.589822 51.324047,....))，点保存为：POINT (4.440742, 51.179215)

我试过这个：

# Fonction pour trouver les INS9 des points dans chaque zone
def find_points_in_zones(df_commerce, df_secteur):
    results = {}
    for commerce_index, row_commerce in df_commerce.iterrows():
        zone = row_commerce['zone_population']
        print(f"Zone {commerce_index}: {zone}")
        points_in_zone = []
        for secteur_index, row_secteur in df_secteur.iterrows():
            point = row_secteur['centroid_degres']
            print(f"    Point {secteur_index}: {point}")
            if zone.contains(point):
                print('il y a un point!')
                points_in_zone.append(row_secteur['INS9'])
        results[commerce_index] = points_in_zone
    return results

points_in_zones = find_points_in_zones(filtered_commerce, filtered_secteur)

# Afficher les résultats
for zone_index, points in points_in_zones.items():
    print(f"Zone {zone_index}: {points}")

但是好像没有区域包含点。我猜问题是“zone.contains(point)”，因为我们从来没有输入过这个if。我也尝试过“.intercepts”，但没有成功。

在项目的其他部分，我使用了类似的功能，但是当我在这里尝试类似的功能时，它不起作用：

polygone_concurrence = isochrone_commerce_double.iloc[0]['geometry']
filtered_commerce.at[index, 'zone_concurrence'] = polygone_concurrence
concurrents = filtered_commerce[filtered_commerce.geometry.within(polygone_concurrence)]    
concurrents_ids = concurrents['UNITID'].tolist()
filtered_commerce.at[index, 'concurrence'] = concurrents_ids

谢谢您的帮助。

Answer 1

创建地理数据框和多边形的空间连接点数据：

import geopandas as gpd
from shapely import wkt
import matplotlib.pyplot as plt

polygon_file = r"/home/bera/Desktop/GIStest/filtered_commerce.csv"
poly = gpd.read_file(polygon_file)
# poly.head(1)
#   zone_id                                               area geometry
# 0       1  MultiPolygon (((20.20966752 63.85528855, 20.26...     None

point_file = r"/home/bera/Desktop/GIStest/filtered_sector.csv"
point = gpd.read_file(point_file)
# point.head(1)
#   INS9                 centroid_degres geometry
# 0    1  Point (20.25524523 63.8059068)     None

#Create geometries, and set coordinate system
poly["geometry"] = poly.apply(lambda x: wkt.loads(x["area"]), axis=1).set_crs(4326)
point["geometry"] = point.apply(lambda x: wkt.loads(x["centroid_degres"]), axis=1).set_crs(4326)

fig, ax = plt.subplots(figsize=(15, 15))
poly.plot(ax=ax, color="blue", edgecolor="white")
point.plot(ax=ax, color="black")

#Spatial Join points to polygons
sjoin = gpd.sjoin(left_df=poly, right_df=point[["INS9","geometry"]], how="left")
sjoin.head()
#   zone_id                                               area  ... index_right  INS9
# 0       1  MultiPolygon (((20.20966752 63.85528855, 20.26...  ...         2.0     3
# 0       1  MultiPolygon (((20.20966752 63.85528855, 20.26...  ...        24.0    25
# 0       1  MultiPolygon (((20.20966752 63.85528855, 20.26...  ...        10.0    11
# 1       2  MultiPolygon (((20.29507989 63.85004948, 20.33...  ...         6.0     7
# 1       2  MultiPolygon (((20.29507989 63.85004948, 20.33...  ...        12.0    13