How to create a small polygon vector layer from a larger layer efficiently in QGIS?

Question

I am creating a plugin for QGIS in python which needs to load in parts of a larger polygon vector layer. This layer is taken from a shapefile and then a new smaller vector layer is created by adding all the polygons(features) with the same ID as an already loaded point layer. The way I am doing this right now is like this:

for file in [f for f in listdir(folderpath) if isfile(join(folderpath, f))]:
#Ignore the fact that it will load all shapefiles of a folder
if file[-4:] == ".shp":
    large = QgsVectorLayer(folderpath + "\\" + file, file[:-4], "ogr")
    #Modify layer to only contain relevant polygons
    pointLayer = QgsMapLayerRegistry.instance().mapLayersByName("Points")[0]
    features = pointLayer.getFeatures()
    newFeatures = []
    for f in features:
        expr = QgsExpression( "\"PolygonID\"='{}'".format(f['PointID']))
        for polyF in large.getFeatures( QgsFeatureRequest(expr)):
            newFeatures.append(polyF)
    small = QgsVectorLayer("Polygon?crs=epsg:4326&field=PolygonID:string(40)", large.name(), "memory")
    pr = small.dataProvider()
    pr.addFeatures(newFeatures)
    small.updateExtents()
    #Load in the "small polygon layer" in the bottom of the layer table
    QgsMapLayerRegistry.instance().addMapLayer(small, False)
    QgsProject.instance().layerTreeRoot().addLayer(small)

and it works, however it is extremely slow (because of large the larger layers are). There might be a much easier way to do this, however this my first ever project working with both GIS and Python.

Answer 1

Since my previous answer didn't help you, I completely edited it because I propose a new approach.

If I have understood your question, you simply want to select all the features from the polygon which store the same ID as the point layer and save them in a new layer. Well, you may solve your problem with pure Python: instead of parsing thousands of features with thousands of values (which need huge time consuming), you may simple recur to dictionaries and sets from their values (see in the following).

I created two sample shapefiles:

• one point layer (containing 4200 features) which stores random values from 0 to 4200 in the "PointID" field;
• one polygon layer (containing 7000 features) which stores random values from 0 to 9000 in the "PolygonID" field.

Using your code (I edited it a little in the first line and I added timing):

import os
from datetime import datetime
import time

start_time = datetime.now()
folderpath = 'C:/Users/xxx/Desktop/'
for file in [f for f in os.listdir(folderpath)]:
    if file[-4:] == ".shp":
        large = QgsVectorLayer(folderpath + "\\" + file, file[:-4], "ogr")
        #Modify layer to only contain relevant polygons
        pointLayer = QgsMapLayerRegistry.instance().mapLayersByName("Points")[0]
        features = pointLayer.getFeatures()
        newFeatures = []
        for f in features:
            expr = QgsExpression( "\"PolygonID\"='{}'".format(f['PointID']))
            for polyF in large.getFeatures( QgsFeatureRequest(expr)):
                newFeatures.append(polyF)
        small = QgsVectorLayer("Polygon?crs=epsg:4326&field=PolygonID:string(40)", large.name(), "memory")
        pr = small.dataProvider()
        pr.addFeatures(newFeatures)
        small.updateExtents()
        #Load in the "small polygon layer" in the bottom of the layer table
        QgsMapLayerRegistry.instance().addMapLayer(small, False)
        QgsProject.instance().layerTreeRoot().addLayer(small)

end_time = datetime.now()
print (end_time - start_time)

the duration time for the analysis was 0:02:55.270000 (more or less 3 minutes).

Using my solution, instead:

import os
from datetime import datetime
import time

start_time = datetime.now()
folderpath = 'C:/Users/xxx/Desktop/'
for file in [f for f in os.listdir(folderpath)]:
    if file[-4:] == ".shp":
        large = QgsVectorLayer(folderpath + "\\" + file, file[:-4], "ogr")
        dict_a = {}
        dict_feat_a = {}
        for feat in large.getFeatures():
            dict_a[feat.id()] = feat['PolygonID']
            dict_feat_a[feat.id()] = feat

        pointLayer = QgsMapLayerRegistry.instance().mapLayersByName("Points")[0]
        dict_b = {}
        for ft in pointLayer.getFeatures():
            dict_b[ft.id()] = ft['PointID']

        int_list = list(set(dict_a.values()) & set(dict_b.values()))
        newFeatures = []
        for k, v in dict_a.iteritems():
            if v  in int_list:
                newFeatures.append(dict_feat_a[k])

        small = QgsVectorLayer("Polygon?crs=epsg:4326&field=PolygonID:string(40)", large.name(), "memory")
        pr = small.dataProvider()
        pr.addFeatures(newFeatures)
        small.updateExtents()
        #Load in the "small polygon layer" in the bottom of the layer table
        QgsMapLayerRegistry.instance().addMapLayer(small, False)
        QgsProject.instance().layerTreeRoot().addLayer(small)

end_time = datetime.now()
print (end_time - start_time)

the duration time for the analysis was 0:00:00.562000 (less than 1 second!).