1

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.

  • The pointLayer layer is always the same? – mgri Jan 13 '17 at 14:07
  • Yes, the pointLayer is always the same, well it will always have the same name and same attribute name. – Felza Jan 13 '17 at 14:15
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!).

  • Well, it took almost 3 times the time, and also didn't work. No features were added to the polygon layer. If it is any help knowing, the data I'm testing with has 8300 features in the larger polygon layer and the point layer has 4200 features. – Felza Jan 13 '17 at 14:59
  • Try writing %f instead of %s. However, I'll do some test on my own because you didn't provide any sample file – mgri Jan 13 '17 at 15:12
  • The ID is a string. I would provide some test data, however the data I have is sensitive and as it is my first time working with GIS I don't know how to generate my own test data large enough to test with. Sorry for the inconvenience, but thank you for helping. – Felza Jan 13 '17 at 15:25
  • @Felza please have a look to my new answer and let me know if it works for you. – mgri Jan 14 '17 at 0:18

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.