geopandas junção espacial extremamente lenta

Estou usando o código abaixo para encontrar um país (e algumas vezes estado) para milhões de pontos de GPS. Atualmente, o código leva cerca de um segundo por ponto, o que é incrivelmente lento. O shapefile é de 6 MB.

Eu li que as geopandas usam rtrees para junções espaciais, tornando-as incrivelmente eficientes, mas isso não parece funcionar aqui. O que estou fazendo de errado? Eu esperava mil pontos por segundo, mais ou menos.

O shapefile e o csv podem ser baixados aqui (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

adicionar o argumento op = 'inside' na função sjoin acelera drasticamente a operação point-in-polygon.

O valor padrão é op = 'intersects', o que eu acho que também levaria ao resultado correto, mas é 100 a 1000 vezes mais lento.

A pergunta pergunta como tirar proveito do r-tree nas junções espaciais das geopandas, e outro respondedor indica corretamente que você deve usar 'dentro' em vez de 'cruza'. No entanto, você também pode tirar proveito de um índice espacial de árvore r em geopandas enquanto usa intersects/ intersection, conforme demonstrado neste tutorial sobre árvore r de geopandas :

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

O que é provável acontecendo aqui é que apenas a trama de dados à direita é alimentado no índice RTREE: https://github.com/geopandas/geopandas/blob/master/geopandas/tools/sjoin.py#L48-L55 que para um op="intersects"executar significaria que o polígono foi alimentado no índice; portanto, para cada ponto, o polígono correspondente é encontrado através do índice rtree.

Mas para op="within", os quadros de geodata são invertidos, pois a operação é realmente o inverso de contains: https://github.com/geopandas/geopandas/blob/master/geopandas/tools/sjoin.py#L41-L43

Então, o que aconteceu quando você trocou opde op="intersects"para op="within"é que, para cada polígono, os pontos correspondentes são encontrados no índice rtree, que no seu caso acelerou a consulta.