12

I am using the code below to find a country (and sometimes state) for millions of GPS points. The code currently takes about one second per point, which is incredibly slow. The shapefile is 6 MB.

I read that geopandas uses rtrees for spatial joins, making them incredibly efficient, but this does not seem to work here. What am I doing wrong? I was hoping for a thousand points per second or so.

The shapefile and csv can be downloaded here (5MB): https://www.dropbox.com/s/gdkxtpqupj0sidm/SpatialJoin.zip?dl=0

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 time


#parameters
shapefile="K:/.../Shapefiles/Used/World.shp"
df=pd.read_csv("K:/.../output2.csv",index_col=None,nrows=20)# Limit to 20 rows for testing    

if __name__=="__main__":
    start=time.time()
    df['geometry'] = df.apply(lambda z: Point(z.Longitude, z.Latitude), axis=1)
    PointsGeodataframe = gpd.GeoDataFrame(df)
    PolygonsGeodataframe = gpd.GeoDataFrame.from_file(shapefile)
    PointsGeodataframe.crs = PolygonsGeodataframe.crs
    print time.time()-start
    merged=sjoin(PointsGeodataframe, PolygonsGeodataframe, how='left')
    print time.time()-start
    merged.to_csv("K:/01. Personal/04. Models/10. Location/output.csv",index=None)
    print time.time()-start
16

adding the argument op='within' in the sjoin function speeds up the point-in-polygon operation dramatically.

Default value is op='intersects', which I guess would also lead to correct result, but is 100 to 1000 times slower.

7

The question asks how to take advantage of r-tree in geopandas spatial joins, and another responder correctly points out that you should use 'within' instead of 'intersects'. However, you can also take advantage of an r-tree spatial index in geopandas while using intersects/intersection, as demonstrated in this geopandas r-tree tutorial:

spatial_index = gdf.sindex
possible_matches_index = list(spatial_index.intersection(polygon.bounds))
possible_matches = gdf.iloc[possible_matches_index]
precise_matches = possible_matches[possible_matches.intersects(polygon)]
5

What's likely going on here is that only the dataframe on the right is fed into the rtree index: https://github.com/geopandas/geopandas/blob/master/geopandas/tools/sjoin.py#L48-L55 Which for an op="intersects" run would mean the Polygon was fed into the index, so for every point, the corresponding polygon is found through the rtree index.

But for op="within", the geodataframes are flipped since the operation is actually the inverse of contains: https://github.com/geopandas/geopandas/blob/master/geopandas/tools/sjoin.py#L41-L43

So what happened when you switched the op from op="intersects" to op="within" is that for every polygon, the corresponding points are found through the rtree index, which in your case sped up the query.

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.