How to find which points intersect with a polygon in geopandas?

Thomas Pingel picture Thomas Pingel · May 22, 2015 · Viewed 22.2k times · Source

I've been trying to use the "intersects" feature on a geodataframe, looking to see which points lie inside a polygon. However, only the first feature in the frame will return as true. What am I doing wrong?

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)

Answer

Fabzi picture Fabzi · May 25, 2015

According to the documentation:

Binary operations can be applied between two GeoSeries, in which case the operation is carried out elementwise. The two series will be aligned by matching indices.

Your examples are not supposed to work. So if you want to test for each point to be in a single polygon you will have to do:

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

Outputs:

    0    True
    1    True
    2    True
    dtype: bool

Or if you want to test for all geometries in a specific GeoSeries:

points.intersects(poly.unary_union)

Geopandas relies on Shapely for the geometrical work. It is sometimes useful (and easier to read) to use it directly. The following code also works as advertised:

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

You might also have a look at How to deal with rounding errors in Shapely for issues that may arise with points on boundaries.