Creating spatial join between points and polygons in GeoPandas?

Question

I want to create a simple spatial join between points and polygons using GeoPandas but I think GeoPandas has bug?

GeoPandas code:

from geopandas import gpd 
import geopandas 
points = geopandas.GeoDataFrame.from_file('points.shp') # or geojson etc 
polys = geopandas.GeoDataFrame.from_file('polygons.shp') 
pointInPoly = gpd.sjoin(points, polys, how='left',op='within') 

Error :

Traceback (most recent call last): 
  File "/home/username/testshapely/sumpointsinsidepolygon/testgeo.py", line 7, in <module> 
    pointInPoly = gpd.sjoin(points, polys, how='left',op='within') 
  File "/usr/local/lib/python2.7/dist-packages/geopandas/tools/sjoin.py", line 57, in sjoin 
    r_idx = np.concatenate(idxmatch.values) 
ValueError: need at least one array to concatenate 

And if I change the imports with the some code:

import geopandas 
import pandas as pd 
import geopandas as gpd 
from geopandas import GeoDataFrame, read_file 
from geopandas.tools import sjoin 
from shapely.geometry import Point, mapping,shape 
import pandas as gpd 

I get the following error:

    pointInPoly = gpd.sjoin(points, polys, how='left',op='within') 
AttributeError: 'module' object has no attribute 'sjoin' 

Any ideas why?

Answer 1

There is not any bug in 'sjoin' method. To realize this, you also need Rtree python module (in my case installed with easy_install) and libspatialindex library (from my Debian Linux repository). After installation of these libraries, I ran my adapted version of your code:

from geopandas import gpd 
points = gpd.GeoDataFrame.from_file('/home/zeito/pyqgis_data/random_points.shp') # or geojson etc 
polys = gpd.GeoDataFrame.from_file('/home/zeito/pyqgis_data/polygon8.shp') 
pointInPoly = gpd.sjoin(points, polys, op='within') 

and it worked.

This is the head of pointInPoly object:

pointInPoly.head()
   Crime_Type  Year                                     geometry  id  value  \
8     Robbery  2012   POINT (416395.3002364286 4459234.39989487)   8   1552   
12    Robbery  2015  POINT (424032.4184008946 4450754.956628661)  12   1411   
18      Theft  2011   POINT (412338.806060733 4456573.265242064)  18   1477   
14    Robbery  2011  POINT (471673.6468123283 4456495.133554206)  14   2474   
11      Theft  2015  POINT (364073.5052346897 4438710.980806596)  11   1919   

    index_right  FID          area  cno  counts  hedgerow        l   max  \
8             1    1  2.805022e+08    1     500         1  65657.8  2330   
12            1    1  2.805022e+08    1     500         1  65657.8  2330   
18            1    1  2.805022e+08    1     500         1  65657.8  2330   
14            2    2  3.737159e+08    2     100        68  76584.0  2962   
11            5    5  1.448294e+08    5      20         2  47325.7  2571   

     min         n  perimeter  values  
8   1351  131315.6  65657.871      15  
12  1351  131315.6  65657.871      15  
18  1351  131315.6  65657.871      15  
14  1438  153168.0  76584.033      15  
11  1534   94651.4  47325.772      15

By using the QuickWKT plugin of QGIS, it can be observed that the point with id = 8 (in red) is one of the all six possibilities (see next image).

Editing Note:

Saving pointInPoly as shapefile:

pointInPoly.to_file('/home/zeito/pyqgis_data/pointInPoly.shp')

and loaded it (in red) at the Map Canvas of QGIS:

Now, it's easier to observe that there are only six points that match op='whithin' sjoin parameter (there is one missed point at head and I do not why).

On the other hand, at next image, you can observe one indirectly way to get the count (into 'index_righ' field created by default) for each polygon feature.