CODEWITHSUNDEEP
pythongeographygeopandas
Thomas Pingel
Thomas Pingel

Reputation: 231

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

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)

Upvotes: 19

Views: 43671

Answers (5)

S.Perera
S.Perera

Reputation: 894

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

Upvotes: 3

Magnus
Magnus

Reputation: 377

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

Upvotes: 13

Fabzi
Fabzi

Reputation: 666

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.

Upvotes: 20

Ioannis Nasios
Ioannis Nasios

Reputation: 8527

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

Upvotes: 2

Thomas Pingel
Thomas Pingel

Reputation: 231

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)

Upvotes: 3

Related Questions