3

I'm building a Qgis 2.18 plugin with python. I have to check if at least one geometry of a linestring vector layer (hundreds of features) intersect the geometries of other two linestring vector layers (hundreds of thousands features). I've written my algorithm, and it works but my problem is the time of execution of the algorithm.

I'm asking if there is a way to speed up a bit my code, maybe using other python data structures, modules or pyqgis functions that I haven't discover yet. I put all my features in a list for every layer to iterate over them.

A snippet of my code:

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)

Any suggestion? I've read abuot Shapely or GeoPandas but I'm not sure if they'll help me.

3

Probably is a very common task on GIS world, so I post my code solution.

With spatial indices as suggested by Joseph I bring down algorithm time from 24 minutes to 80 seconds.

First scratch of a solution:

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)

It's a draft, but works. Obviusly that code could be a little more improved and be more pythonic.

  • 1
    Thanks for posting your solution! You certainly improved the time taken :) – Joseph Mar 7 '18 at 10:02

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.