3

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?

  • To run this method you also need Rtree python module and libspatialindex library. – xunilk Jan 14 '17 at 19:52
  • Your 2nd code snippet raises that exception because you overwrite the geopandas import wih pandas import pandas as gpd – user2856 Jan 14 '17 at 20:19
  • yes that mistake but not work again – Chris Papas Jan 14 '17 at 20:26
7

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).

enter image description here

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:

enter image description here

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.

enter image description here

  • the counts field is the points where is within in polygons? – Chris Papas Jan 14 '17 at 20:15
  • No. It was created to generate random points for my answer in another question in this site. – xunilk Jan 14 '17 at 21:05
  • ,but this spatial join dont count the points where is within in polygons?can i export pointInPoly to new shapefile ? – Chris Papas Jan 14 '17 at 21:12
  • Yes, you can (see my Editing Note) – xunilk Jan 14 '17 at 21:20
  • thnx you are so good,the new shapefile is the point geometry?on the first spatial join how to export the polygons where greater than 1 point or 0 points ?i am comfused now haha – Chris Papas Jan 14 '17 at 21:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.