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

Question

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)

score 17 · Accepted Answer · edited Mar 31 '19 at 21:51

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.

This helped very much in explaining what was going on. We're actually using a GeoDataFrame, and so I was casting the problem in a somewhat more simple way. Your example works very well, but it's more convenient for us to keep the polygons (and their related attributes) in a GeoDataFrame for other reasons, otherwise the shapely way works nicely. — Thomas Pingel, May 27 '15 at 02:56
Good! I added your example to my answer which I find to be more comprehensive (with the ref to the doc). You might consider marking the question as solved and/or upvoting. — Fabzi, May 27 '15 at 18:03

Magnus · Answer 2 · 2023-04-14T11:55:37.553

11

Since geopandas underwent many performance-enhancing changes recently, answers here are outdated. Geopandas 0.8 introduced many changes that makes handling large datasets a lot faster.

    import geopandas as gpd
    from shapely.geometry import Polygon, Point
    
    points = gpd.GeoSeries([Point(.5,.5), Point(.5,1), Point(1,1), Point(10,10)])
    polys = gpd.GeoSeries([Polygon([(0,0), (0,2), (2,2), (2,0)])])
    
    point_gdf = gpd.GeoDataFrame({'geometry': points})
    poly_gdf = gpd.GeoDataFrame({'geometry': polys})
    
    gpd.overlay(point_gdf, poly_gdf, how='intersection')

edited Apr 14 '23 at 11:55

answered Nov 23 '20 at 13:45

Magnus

305
3
12

This is the only solution that works cleanly for me. – clifgray Dec 30 '20 at 22:19
1

NotImplementedError: overlay currently only implemented for GeoDataFrames – Trenton Mar 21 '23 at 19:46
Thank's for the hint. Geopandas must have dropped support for GeoSeries in the overlay function between 0.8 and 0.12. I've updated the answer so the GeoSeries are turned into GeoDataframes and added a point that is not within the Polygon for testing. – Magnus Apr 14 '23 at 11:57

score 3 · Answer 3 · answered May 27 '15 at 02:54

One way around this seems to be to either get a particular entry (which doesn't work for my application, but may work for someone else's:

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

Another way (more useful for my application) is to intersect with a unary union of the features for the second layer:

points.intersects(poly.unary_union)

score 3 · Answer 4 · edited Aug 18 '23 at 15:55

3

I think the fastest way to do it is by using geopandas.sjoin.

import geopandas as gpd

gpd.sjoin(pts, poly, how='left', predicate='intersects')

Check example: link

edited Aug 18 '23 at 15:55

adrien.ludwig

3
4

answered May 18 '22 at 15:56

S.Perera

874
3
9
24

score 2 · Answer 5 · answered Jan 17 '19 at 16:44

You can easily check which points lies inside a polygon with this simple function below:

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

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

5 Answers5