2

I have the following problem: I have a list of shapely points and a list of shapely polygons. Now I want to check in which polygon a given point is.

At the moment I am using the following code, which seems not very clever:


# polygons_df is a pandas dataframe that contains the geometry of the polygons and the usage of the polygons (landuses in this case, e.g. residential)

# point_df is a pandas dataframe that contains the geometry of the points and the usage of the point (landuses in this case, e.g. residential)

# polylist is my list of shapely polygons

# pointlist is my list of shapely points 

from shapely.geometry import Point, Polygon
import pandas as pd
import geopandas as gpd

i = 0
while i < len(polygons_df.index):
    j = 0
    while j < len(point_df.index):
        if polylist[i].contains(point):
            point.at[j, 'tags.landuse'] = polygons_df.iloc[i]['tags.landuse']
        else:
            pass
        j += 1
    i += 1

Can I somehow speed this up? I have more than 100.000 points and more than 10.000 polygons and these loops take a while. Thanks!

Kevin
  • 161
  • 9
  • Are `polylist` the geometries of `polygons_df`, and `pointlist` the geometries of `point_df`? And there is also a `point` variable you didn't explain (general note: it is always recommended to provide a reproducible code example, one that can run on itself) – joris Apr 18 '19 at 15:48
  • 1
    But in general, I think you are looking for a "point-in-polgyon" spatial join, which you can do with geopandas.sjoin function: `geopandas.sjoin(point_df, polygons_df, op='within')` (see https://geopandas.readthedocs.io/en/latest/mergingdata.html#spatial-joins) – joris Apr 18 '19 at 15:50
  • To speed-up the operation, consult this: https://geoffboeing.com/2016/10/r-tree-spatial-index-python/. – swatchai Apr 19 '19 at 01:29
  • @swatchai the `geopandas.sjoin` spatial join uses an rtree spatial index under the hood for you – joris Apr 19 '19 at 07:41
  • Thank you guys very much, the geopandas.sjoin function is doing exactly what I want! – Kevin Apr 23 '19 at 11:41
  • Possible duplicate of [Looking for a fast way to find the polygon a point belongs to using Shapely](https://stackoverflow.com/questions/20297977/looking-for-a-fast-way-to-find-the-polygon-a-point-belongs-to-using-shapely) – Georgy May 31 '19 at 12:30

1 Answers1

1

I know a solution was found in the comments for the particular problem, but to answer a related question of how to check if an array of points is inside a shapely Polygon, I found the following solution:

>>> poly = Polygon([(0,0), (1,0), (0,1)])
>>> contains = np.vectorize(lambda p: poly.contains(Point(p)), signature='(n)->()')
>>> contains(np.array([[0.5,0.49],[0.5,0.51],[0.5,0.52]]))
array([ True, False, False])

I don't know that this neccesarily speeds up the calculations, but at least you can avoid the for-loop.

Dith
  • 151
  • 8