计算多边形中的点并将结果写入 (Geo)Dataframe

Posted 2023-03-23

【中文标题】计算多边形中的点并将结果写入 (Geo)Dataframe

【英文标题】：Count Points in Polygon and write result to (Geo)Dataframe

【发布时间】：2021-12-07 04:31:15

【问题描述】：

我想计算每个多边形有多少个点

# Credits of this code go to: https://***.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

现在我sjoin这两个。我首先使用df_poly，以便将点dfp 添加到GeoDataframe df_poly。

df_poly.sjoin(dfp)

现在我想计算每个polygon 有多少points。 我以为

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

但这不会将column 添加到GeoDataframe df_poly 以及每个count 的count。

【参考方案1】：

您需要使用合并将count() 的输出中的一列添加回原始DataFrame。我使用了几何列并将其重命名为n_points：

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

【讨论】：

此答案有效，但您能否向希望获得理解的人解释一下？

【参考方案2】：

这是这个问题的后续The indices of the two GeoSeries are different - Understanding Indices

空间连接的

right_index

给出了多边形的索引，因为多边形位于空间连接的右侧 因此可以将系列gpd.sjoin(dfp, df_poly).groupby("index_right").size().rename("points") 简单地连接到多边形

GeoDataFrame

以给出找到的点数 注意how="left" 以确保它是左连接，而不是内连接。在这种情况下，任何没有点的多边形都有NaN，您可能需要fillna(0)。

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",
)

【参考方案3】：

基于 Fergus McClean 提供的答案，这甚至可以用更少的代码完成：

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

但是，Rob Raymond 提出的将两个dataframes 结合起来的方法（.join()）保留了没有计数的条目。