10

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?

Taras
  • 32,823
  • 4
  • 66
  • 137
Stephen
  • 385
  • 1
  • 4
  • 13

4 Answers4

24

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)]
Matt Payne
  • 41
  • 1
  • 8
joris
  • 3,903
  • 22
  • 28
5

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.

4

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')
sash_wash
  • 51
  • 2
3

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

Jane Kathambi
  • 131
  • 2