Looking for GeoPandas QGIS equivalent of ArcGIS "Select by location" tool

Question

I have two GeoDataFrames:

1. A single multi-part polygon
2. Many points

I have read them in with the following:

import geopandas as gpd

poly = gpd.read_file('C:\\Users\\srcha\\Desktop\\folder\\poly.shp')
points = gpd.read_file('c:\\Users\\srcha\\Desktop\\folder\\points.shp')

I want to run the equivalent of "Select by location" in ArcGIS with GeoPandas and shapely tools to select points that are within the polygon.

How should I go about this?

Answer 1

If poly is a GeoDataFrame with a single geometry, extract this:

polygon = poly.geometry[0]

Then, you can use the within method to check which points are within the polygon:

points.within(polygon)

returning a boolean True/False values which can be used to filter the original dataframe:

subset = points[points.within(polygon)]

Answer 2

In case the polygon layer contains many polygon features, the solution can be extended, like so:

poly = gpd.read_file('C:/Users/srcha/Desktop/folder/poly.shp')
points = gpd.read_file('c:/Users/srcha/Desktop/folder/points.shp')

poly['dummy'] = 'dummy'  # add dummy column to dissolve all geometries into one
geom = poly.dissolve(by='dummy').geometry[0]  # take the single union geometry
subset = points[points.within(geom)]

I'm not fond of the dummy-trick, but I have no other idea how to do this using geopandas.

Answer 3

It is much faster if you use spatial join:

subset = gpd.sjoin(points, polygon, how='inner', op='within')

GeoPandas Version 0.10.0 or later:

subset = gpd.sjoin(points, polygon, how='inner', predicate='within')

Answer 4

If you convert the poly boundary GeoDataFrame to a single polygon, eg with:

poly_union = poly.geometry.unary_union

then you can use boolean filtering with the within operator to select points within and outside of poly_union:

within_poly_union = points[points.geometry.within(poly_union)]
outside_poly_union = points[~points.geometry.within(poly_union)]

Alternatively, you could use the disjoint spatial predicate:

outside_poly_union = points[points.geometry.disjoint(poly_union)]
inside_poly_union = points[~points.geometry.disjoint(poly_union)]

Answer derived from here