How to find which points are contained in a polygon or multipolygon in python?

Question

I have 2 tables, the first is called 'filtered_commerce' and the 2nd is called 'filtered_secteur'. The 'filtered_commerce' table contains a 'zone_population' column which contains polygons and multipolygons representing geographic areas. The filtered_sector table contains a 'centroid_degres' column which contains points and an 'INS9' column which contains the identifiers of the points. I would like to find which points are in each zone and save for each zone the INS9 of the points included in it in a new column of the 'filtered_commerce'. Areas are polygons or multipolygons saved as follows: POLYGON ((4.589822 51.324047,....)) and points are saved as: POINT (4.440742, 51.179215)

I tried this :

# 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}")

but it is as if no area contained points. I guess the problem is 'zone.contains(point)' because we never enter in this if. I also tried the '.intercepts' and it didn't worked.

In an other part of the project, I used a similar function but when I tried something similar here, it didn't work :

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

Thank you for your help.

Answer 1

It looks like you have WKT strings.

You can convert them to real geometries with shapely, create geodataframes and spatial join point data to the polygons:

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