1

I an deleting one of the intersecting polygons from a vector layer using PyQGIS code mentioned below. It is too slow as I have many features in the layer. Is there a way to speed it up?

index= QgsSpatialIndex(layer.getFeatures())
request = QgsFeatureRequest()
clause = QgsFeatureRequest.OrderByClause('id', ascending=True)
orderby = QgsFeatureRequest.OrderBy([clause])
request.setOrderBy(orderby)
for feat in layer.getFeatures(request):
    feats=[f for f in layer.getFeatures(request)]
    for f in feats:
        if index.intersects(f.geometry().boundingBox()):
            layer.select(feat.id())
 selIDs=layer.selectedFeatureIds()
 layer.startEditing()
 layer.deleteFeatures(selIDs1)
 layer.commitChanges()
2
  • 1
    Have you already read this answer? gis.stackexchange.com/questions/411045/… It gives you several recommendations. Commented Dec 20, 2021 at 16:37
  • I tried to implement QgsRequest and QgsFeatureIndex but I think it's not helping much as I need find intersection between features of the same layer or I am not getting it correct it seems.- @Germán Carrillo Commented Dec 20, 2021 at 17:54

1 Answer 1

5

@GermánCarrillo already provided a pretty useful link to a great answer. Here is something more specific about your script, where there are a few issues:

  1. Since you are looking for intersecting features, there is no need to sort the layer by ids (or is there something you did not mention?)
  2. You are iterating over the whole layer for each feature and checking the index (incorrectly) afterwards, so there is pretty no performance gain in the index usage. Instead, for each feature, iterate over the index only.
  3. You can store the features geometries directly inside the index by using flags=QgsSpatialIndex.FlagStoreFeatureGeometries so you can avoid an inperformant feature request and instead get the geometry from the index much faster
  4. Dont forget to exclude self-intersections; otherwise your entire layer will be deleted
  5. If the index matches this does not mean that the features actually intersect, that only means the bounding box intersects, so you need to do an actual intersects-test on the matches.
  6. There is no need to do a selection if you only want to delete the features

Try this:

layer = iface.activeLayer()
index = QgsSpatialIndex(layer.getFeatures(), flags=QgsSpatialIndex.FlagStoreFeatureGeometries)
featids = []
for feat in layer.getFeatures():
    for fid in index.intersects(feat.geometry().boundingBox()): # iterate over the index-matches. The index returns the IDs of the features where the boundingbox intersects
        if fid == feat.id(): # ignore self-intersections
            continue
        fgeom = index.geometry(fid) # get the feature geometry by the id from the index
        if fgeom.intersects(feat.geometry()): # now check if not only the bounding box intersects, but if the actual features geometries intersect
            featids.append(feat.id()) # if so append the id to a list for selection afterwards
layer.startEditing()
layer.deleteFeatures(featids)
layer.commitChanges()

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.