如何找到与 geopandas 中的多边形相交的点？答案

【问题标题】：How to find which points intersect with a polygon in geopandas?如何找到与 geopandas 中的多边形相交的点？
【发布时间】：2015-08-05 00:07:13
【问题描述】：

我一直在尝试在地理数据框中使用“相交”功能，以查看哪些点位于多边形内。但是，只有帧中的第一个特征会返回为真。我做错了什么？

from geopandas.geoseries import *

p1 = Point(.5,.5)
p2 = Point(.5,1)
p3 = Point(1,1)

g1 = GeoSeries([p1,p2,p3])
g2 = GeoSeries([p2,p3])

g = GeoSeries([Polygon([(0,0), (0,2), (2,2), (2,0)])])

g1.intersects(g) # Flags the first point as inside, even though all are.
g2.intersects(g) # The second point gets picked up as inside (but not 3rd)

【问题讨论】：

标签： python geography geopandas

【解决方案1】：

根据documentation：

二元运算可以在两个 GeoSeries 之间应用，在这种情况下该操作是按元素进行的。这两个系列将通过匹配索引对齐。

你的例子不应该工作。因此，如果您想测试每个点是否位于单个多边形中，您必须这样做：

poly = GeoSeries(Polygon([(0,0), (0,2), (2,2), (2,0)]))
g1.intersects(poly.ix[0])

输出：

    0    True
    1    True
    2    True
    dtype: bool

或者，如果您想测试特定 GeoSeries 中的所有几何图形：

points.intersects(poly.unary_union)

Geopandas 依靠 Shapely 进行几何工作。有时直接使用它是有用的（并且更容易阅读）。以下代码也可以像宣传的那样工作：

from shapely.geometry import *

p1 = Point(.5,.5)
p2 = Point(.5,1)
p3 = Point(1,1)

poly = Polygon([(0,0), (0,2), (2,2), (2,0)])

for p in [p1, p2, p3]:
    print(poly.intersects(p))

你也可以看看 How to deal with rounding errors in Shapely 对于边界上的点可能出现的问题。

【讨论】：

这对解释发生了什么非常有帮助。我们实际上使用的是 GeoDataFrame，所以我用一种更简单的方式来解决这个问题。您的示例效果很好，但是出于其他原因，我们将多边形（及其相关属性）保存在 GeoDataFrame 中更方便，否则这种匀称的方式效果很好。
好！我将您的示例添加到我的答案中，我发现它更全面（参考文档）。您可以考虑将问题标记为已解决和/或赞成。

【解决方案2】：

解决此问题的一种方法似乎是获取特定条目（这不适用于我的应用程序，但可能适用于其他人的应用程序：

from geopandas.geoseries import *

p1 = Point(.5,.5)
p2 = Point(.5,1)
p3 = Point(1,1)

points = GeoSeries([p1,p2,p3])

poly = GeoSeries([Polygon([(0,0), (0,2), (2,2), (2,0)])])

points.intersects(poly.ix[0])

另一种方法（对我的应用程序更有用）是与第二层特征的一元联合相交：

points.intersects(poly.unary_union)

【讨论】：

【解决方案3】：

您可以使用下面这个简单的函数轻松检查多边形内的点：

import geopandas
from shapely.geometry import *

p1 = Point(.5,.5)
p2 = Point(.5,1)
p3 = Point(1,1)

g = Polygon([(0,0), (0,2), (2,2), (2,0)])

def point_inside_shape(point, shape):
    #point of type Point
    #shape of type Polygon
    pnt = geopandas.GeoDataFrame(geometry=[point], index=['A'])
    return(pnt.within(shape).iloc[0])

for p in [p1, p2, p3]:
    print(point_inside_shape(p, g))

【讨论】：

【解决方案4】：

由于 geopandas 最近经历了许多提高性能的变化，这里的答案已经过时了。 Geopandas 0.8 引入了许多更改，使处理大型数据集的速度大大加快。

import geopandas
from shapely import Polygon

p1 = Point(.5,.5)
p2 = Point(.5,1)
p3 = Point(1,1)

points = GeoSeries([p1,p2,p3])

poly = GeoSeries([Polygon([(0,0), (0,2), (2,2), (2,0)])])

geopandas.overlay(points, poly, how='intersection')

【讨论】：

这是唯一适合我的解决方案。