-2

I have two spatial datasets: municipality

municipality geometry
1000         POLYGON ((4.04305 47.56271, 4.04345 47.56303, ...
1001         POLYGON ((-0.24733 45.88792, -0.24715 45.88820...
1002         POLYGON ((-0.30449 45.91087, -0.30446 45.91120...

points

points geometry
200    POINT (3617775.075 3205012.273)
201    POINT (3617529.179 3205007.754)
202    POINT (4375420.232 4093242.822)

I want to estimate how many points are inside each municipality. Some of these points will not be inside any municipality, some municipalities will have several points. Here is my code:

import geopandas as gdp
from shapely.geometry import Point
from shapely.geometry.polygon import Polygon
import pandas as pd

points= gdp.read_file("points.shp").to_crs("epsg:4326")
municipality= gdp.read_file("municipality.shp")

df= gdp.tools.sjoin(municipality, points, predicate="within", how='left')
df.index_right.unique()
array([nan])

I find it very weird that none of the points are inserted in a municipality. What I am doing wrong? Any help will be highly appreciated.

MG Fern
  • 75
  • 9

1 Answers1

1

If you want to have the number of points inside each municipality you can do:

municipality['num_points'] = municipality.geometry.apply(lambda x: x.contains(points.geometry).sum())

If you want to see each point to what municipality it belongs to you can execute:

points_by_municipality = gpd.sjoin(municipality, points).groupby('municipality', group_keys=True).apply(lambda x: x)

Hope it works for you.

kithuto
  • 455
  • 2
  • 11