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)Best Answer
According to the documentation:
Binary operations can be applied between two GeoSeries, in which casethe operation is carried out elementwise. The two series will bealigned 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 True1 True2 Truedtype: boolOr 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 atHow to deal with rounding errors in Shapely for issues that may arise with points on boundaries.
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 gpdfrom shapely.geometry import Polygon, Pointpoints = 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')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)I think the fastest way to do it is by using geopandas.sjoin.
import geopandas as gpdgpd.sjoin(pts, poly, how='left', predicate='intersects')Check example: link
You can easily check which points lies inside a polygon with this simple function below:
import geopandasfrom 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 Polygonpnt = geopandas.GeoDataFrame(geometry=[point], index=['A'])return(pnt.within(shape).iloc[0])for p in [p1, p2, p3]:print(point_inside_shape(p, g))