géopandas jointure spatiale extrêmement lente

J'utilise le code ci-dessous pour trouver un pays (et parfois un état) pour des millions de points GPS. Le code prend actuellement environ une seconde par point, ce qui est incroyablement lent. Le fichier de formes fait 6 Mo.

J'ai lu que les géopandas utilisent rtrees pour les jointures spatiales, ce qui les rend incroyablement efficaces, mais cela ne semble pas fonctionner ici. Qu'est-ce que je fais mal? J'espérais mille points par seconde environ.

Le shapefile et le csv peuvent être téléchargés ici (5 Mo): 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

l'ajout de l'argument op = 'within' dans la fonction sjoin accélère considérablement l'opération de point dans le polygone.

La valeur par défaut est op = 'intersects', ce qui, je suppose, conduirait également à un résultat correct, mais est 100 à 1000 fois plus lent.

La question demande comment tirer parti de r-tree dans les jointures spatiales des géopandas, et un autre intervenant souligne correctement que vous devez utiliser «dedans» au lieu de «intersecte». Cependant, vous pouvez également profiter d'un index spatial r-tree dans les géopandas tout en utilisant intersects/ intersection, comme démontré dans ce tutoriel r-tree 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)]

Ce qui se passe probablement ici, c'est que seule la trame de données de droite est introduite dans l'index rtree: https://github.com/geopandas/geopandas/blob/master/geopandas/tools/sjoin.py#L48-L55 Lequel pour un op="intersects"run signifierait que le polygone a été introduit dans l'index, donc pour chaque point, le polygone correspondant est trouvé via l'index rtree.

Mais pour op="within", les géodonnées sont inversées car l'opération est en fait l'inverse de contains: https://github.com/geopandas/geopandas/blob/master/geopandas/tools/sjoin.py#L41-L43

Donc, ce qui s'est passé lorsque vous avez basculé opde op="intersects"à op="within"est que pour chaque polygone, les points correspondants sont trouvés via l'index rtree, ce qui dans votre cas a accéléré la requête.