Pyqgis iterar a través de cientos de miles de características, el tema de la velocidad

Question

Estoy construyendo una Qgis 2.18 plugin con python. Tengo que comprobar si al menos una geometría de un linestring capa vectorial (cientos de características) se cruzan las geometrías de los otros dos linestring capas vectoriales (cientos de miles). He escrito mi algoritmo, y funciona, pero mi problema es el tiempo de ejecución del algoritmo.

Me estoy preguntando si existe una manera de acelerar un poco el código, tal vez usando python estructuras de datos, módulos o pyqgis funciones que no he descubrir todavía. Puse todas mis funciones en una lista para cada capa para iterar sobre ellos.

Un fragmento de mi código:

lay_to_check=  QgsVectorLayer( 'to_check.shp' "Layer to check" , "ogr" )
lay1 = QgsVectorLayer( 'lay1.shp', "LAY1" , "ogr" )
lay2 = QgsVectorLayer( 'lay2.shp', "LAY2" , "ogr" )

#create two lists of features of lay1 and lay2 that intersects the extent of the layer to check
ext = lay_to_check.extent()
g = lambda x:x.geometry().intersects(ext)


request_lay1= QgsFeatureRequest().setSubsetOfAttributes( [''], lay1.fields() )
list_lay1= filter( g, lay1.getFeatures(request_lay1) )

request_lay2 = QgsFeatureRequest().setSubsetOfAttributes( ['field2'], lay2.fields() )
list_lay2 = filter( g, lay2.getFeatures(request_lay2) )

#for every feature I control if intersect every feature of other two list
request = QgsFeatureRequest().setSubsetOfAttributes( [''], lay_to_check.fields() )
for feature in lay_to_check.getFeatures(request):

    control = []
    geom = feature.geometry()
    p = lambda x:x.geometry().intersects(geom)


    if list_lay1:
        list_int1 = filter( p, list_lay1 )
        if list_int1:
            control.append( True )
        else:
            control.append( False )

    if list_lay2:
         list_int2 = filter( p, list_lay2 )
         if list_int2:
             list_attribute = [ feature['field2'] for feat in list_int2 ]
             control.append( max(list_attribute) )
         else:
             control.append(0)

Alguna sugerencia? He leído ubicado aproximadamente a bien formada o GeoPandas pero no estoy seguro de si nos van a ayudar a mí.

Answer 1

Probablemente es una tarea muy habitual en el mundo GIS, así que he puesto mi código de la solución.

Con índices espaciales como sugiere José I bajar el algoritmo de tiempo de 24 minutos a 80 segundos.

Primer scratch de una solución:

lay_to_check=  QgsVectorLayer( 'to_check.shp' "Layer to check" , "ogr" )
lay1 = QgsVectorLayer( 'lay1.shp', "LAY1" , "ogr" )
lay2 = QgsVectorLayer( 'lay2.shp', "LAY2" , "ogr" )

ext = lay_to_check.extent()
g = lambda x:x.geometry().intersects(ext)


request_lay1= QgsFeatureRequest().setSubsetOfAttributes( [''], lay1.fields() )
list_lay1= filter( g, lay1.getFeatures(request_lay1) )

request_lay2 = QgsFeatureRequest().setSubsetOfAttributes( ['field2'], lay2.fields() )
list_lay2 = filter( g, lay2.getFeatures(request_lay2) )

allfeatures1 = {feature.id(): feature for feature in lay1.getFeatures(request_lay1)}
index1 = QgsSpatialIndex()
map(index1.insertFeature, list_lay1)

allfeatures2 = {feature.id(): feature for feature in lay2.getFeatures(request_lay2)}
index12= QgsSpatialIndex()
map(index2.insertFeature, list_lay2)


request = QgsFeatureRequest().setSubsetOfAttributes( [''], lay_to_check.fields() )
for feature in lay_to_check.getFeatures(request):

    control = []

    if list_lay1:
        ids = index1.intersects(feature.geometry().boundingBox())
        for id in ids:
            f = allfeaturesI[id]
            if feature.geometry().intersects(f.geometry()):
                list_int1.append(True)
        if list_int1:
            control.append( True )
        else:
            control.append( False )

    if list_lay2:
        ids = index2.intersects(feature.geometry().boundingBox())
        for id in ids:
            f = allfeaturesI[id]
            if feature.geometry().intersects(f.geometry()):
                list_int2.append(f)
         if list_int2:
             list_attribute = [ feat['field2'] for feat in list_int2 ]
             control.append( max(list_attribute) )
         else:
             control.append(0)

Es un proyecto, pero funciona. Claro que el código podría ser un poco más mejorado y ser más python.