Count Points in Polygon and write result to (Geo)Dataframe-CodePudding

I want to count how many points there are per Polygon

# Credits of this code go to: https://stackoverflow.com/questions/69642668/the-indices-of-the-two-geoseries-are-different-understanding-indices/69644010#69644010
import pandas as pd
import numpy as np
import geopandas as gpd
import shapely.geometry
import requests

# source some points and polygons
# fmt: off
dfp = pd.read_html("https://www.latlong.net/category/cities-235-15.html")[0]
dfp = gpd.GeoDataFrame(dfp, geometry=dfp.loc[:,["Longitude", "Latitude",]].apply(shapely.geometry.Point, axis=1))
res = requests.get("https://opendata.arcgis.com/datasets/69dc11c7386943b4ad8893c45648b1e1_0.geojson")
df_poly = gpd.GeoDataFrame.from_features(res.json())
# fmt: on

Now I sjoin the two. I use df_poly first, in order to add the points dfp to the GeoDataframe df_poly.

df_poly.sjoin(dfp)

Now I want to count how many points there are per polygon. I thought

df_poly.sjoin(dfp).groupby('OBJECTID').count()

But that does not add a column to the GeoDataframe df_poly with the count of each group.

CodePudding user response：

This is a follow on to this question The indices of the two GeoSeries are different - Understanding Indices

right_index of spatial join gives index of polygon as polygon was on right of spatial join
hence the series gpd.sjoin(dfp, df_poly).groupby("index_right").size().rename("points") can then be simply joined to the polygon GeoDataFrame to give how many points were found
note how="left" to ensure it's a left join, not an inner join. Any polygons with no points with have NaN you may want to fillna(0) in this case.

import pandas as pd
import numpy as np
import geopandas as gpd
import shapely.geometry
import requests

# source some points and polygons
# fmt: off
dfp = pd.read_html("https://www.latlong.net/category/cities-235-15.html")[0]
dfp = pd.concat([dfp,dfp]).reset_index(drop=True)
dfp = gpd.GeoDataFrame(dfp, geometry=dfp.loc[:,["Longitude", "Latitude",]].apply(shapely.geometry.Point, axis=1))
res = requests.get("https://opendata.arcgis.com/datasets/69dc11c7386943b4ad8893c45648b1e1_0.geojson")
df_poly = gpd.GeoDataFrame.from_features(res.json())
# fmt: on

df_poly.join(
    gpd.sjoin(dfp, df_poly).groupby("index_right").size().rename("points"),
    how="left",
)

CodePudding user response：

You need to add one of the columns from the output of count() back into the original DataFrame using merge. I have used the geometry column and renamed it to n_points:

df_poly.merge(
    df_poly.sjoin(
        dfp
    ).groupby(
        'OBJECTID'
    ).count().geometry.rename(
        'n_points'
    ).reset_index())

CodePudding user response：

Building on the answere Fergus McClean provided, this can even be done in less code:

df_poly.merge(df_poly.sjoin(dfp).groupby('OBJECTID').size().rename('n_points').reset_index())

However, the method (.join()) proposed by Rob Raymond to combine the two dataframes keeps the entries that have no count.