1

I have polygons and points.

What I do is a very simple operation - count points in the polygon. Primary I am working in Python, however I use QGIS to check and visualize the results.

I have this code where I using GeoPandas Spatial Join and a loop to find which point belongs to which small polygon. To check I did the same thing in QGIS using Count Points in Polygon tool, but the number of points are a little different from the result I obtained with Python.

import geopandas as gpd

polygons = gpd.read_csv(r'E:\...\Polygons.shp')
points = gpd.read_file(r'E:\...\points.shp')

points_join = gpd.sjoin(points, polygons)[['polygon_ID', 'point_Ht', 'point_D', 'geometry']]

for idx, val in polygons.index:
    df_p = points_join[points_join['polygon_ID'] == val]
    
    number_of_points = df_p['polygon_ID'].count()

    polygons.loc[idx, 'NumberPoints'] = number_of_points
    

Here is the point count using both methods, and as we can see that there are no significant differences between both methods, however the number of points do not match.

enter image description here

enter image description here

Has anyone has experienced something similar working with counting points in polygon in Python and QGIS? At first I assumed that the problem could be the coordinate system, however I use the same files with the same CRS each time.

1
  • 1
    "In" and "On" are different relationships. Have you checked to see if the "missing" points aren't really within any polygons?
    – Vince
    Dec 27, 2021 at 12:55

2 Answers 2

2

I get the same count in geopandas and QGIS. Your code is a bit complicated, try using groupby:

import geopandas as gpd

dfpoint = gpd.read_file(filename=r'/home/bera/Desktop/GIStest/1000points.shp')
dfpoly = gpd.read_file(filename=r'/home/bera/Desktop/GIStest/polygons.shp')
dfsjoin = gpd.sjoin(left_df=dfpoly, right_df=dfpoint, how='left') #Left or polygons with zero points are discarded
dfcount = dfsjoin.groupby('polyid')['pointid'].count().rename('pointcount').reset_index()

#   polyid  pointcount
#0       1          97
#1       2          70
#2       3           0

enter image description here

0
0

@BERA answer is correct and works smoothly, however as I need the loop for further calculations I found another way simply adding op = 'within' when joining points and polygons, like this:

import geopandas as gpd

polygons = gpd.read_csv(r'E:\...\Polygons.shp')
points = gpd.read_file(r'E:\...\points.shp')

points_join = gpd.sjoin(points, polygons, op='within')

for idx in polygons.index:
    df_p = points_join[points_join['polygon_ID'] == idx]
    
    number_of_points = df_p['polygon_ID'].count()

    polygons.loc[idx, 'NumberPoints'] = number_of_points

Point number matches QGIS!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.