unión espacial geopandas extremadamente lenta

Question

Estoy utilizando el siguiente código para encontrar un país (y a veces un estado) para millones de puntos GPS. Actualmente, el código tarda aproximadamente un segundo por punto, lo que es increíblemente lento. El shapefile es de 6 MB.

He leído que geopandas utiliza rtrees para las uniones espaciales, haciéndolas increíblemente eficientes, pero esto no parece funcionar aquí. ¿Qué estoy haciendo mal? Esperaba mil puntos por segundo más o menos.

El shapefile y el csv pueden descargarse aquí (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

Answer 1

Añadir el argumento predicate='within' en la función sjoin acelera drásticamente la operación point-in-polygon.

El valor por defecto es predicate='intersects', que supongo que también llevaría a un resultado correcto, pero es de 100 a 1000 veces más lento.

Answer 2

Lo que probablemente ocurre aquí es que sólo el marco de datos de la derecha se introduce en el índice rtree: https://github.com/geopandas/geopandas/blob/master/geopandas/tools/sjoin.py#L48-L55 Que para un predicate="intersects" ejecutar significaría que el polígono se introdujo en el índice, por lo que para cada punto, el polígono correspondiente se encuentra a través del índice rtree.

Pero para predicate="within" los marcos de datos geográficos se invierten, ya que la operación es en realidad la inversa de contains : https://github.com/geopandas/geopandas/blob/master/geopandas/tools/sjoin.py#L41-L43

Entonces, ¿qué pasó cuando cambió el predicate de predicate="intersects" a predicate="within" es que para cada polígono, los puntos correspondientes se encuentran a través del índice rtree, lo que en tu caso aceleró la consulta.

Answer 3

La pregunta se refiere a cómo aprovechar las ventajas de r-tree en las uniones espaciales de geopandas, y en otra respuesta se señala correctamente que se debe utilizar 'within' en lugar de 'intersects'. Sin embargo, también se puede aprovechar un índice espacial r-tree en geopandas utilizando intersects / intersection como se demuestra en este 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)]

Answer 4

En relación con el tema de las geopandas debatido en esta respuesta y visto ici en la última versión:

Si tiene dos geodataframes y desea realizar un intersects y sabe qué índice espacial debe consultarse, puede unir de pequeño a grande y utilizar how="right" para forzar el uso de este índice:

small_gdf.sjoin(large_gdf, predicate="intersects", how="right")