I have two input layers. The first one is a line-layer named "graph" and the second one is a point-layer named "nodes".
The graph layer holds 177350 line features which are representing a streetmap. Loaded into QGIS it looks like this:
Each street (or line) holds some information: The id of the street (edgeid), the node id of the startnode, the node id of the endnode and some other attributes as shown below:
Every node in the point-layer has a unique nodeid, x-, y- coordinates and a type (the type attribute is not relevant):
What I want to do is to calculate the azimuth for every street (line) in the graph layer. The outputlayer of my python script should look exactly like the graph-layer, just with an additional azimuth attribute. Like this:
To achieve this, I have written an python script which loops through all the features of the graph layer. For every street, I look up the x-,y-coordinates of the startnode, and of the endnode. I will need the x-,y-coordinates for calculate the azimuth later on in the script. I am doing this using QgsFeatureRequest. But as there are 177350 lines to go through, I will have to request 177350 * 2 times (177350 times for startnodes and 177350 times for endnodes). And it just takes too long.
I have already tried to set the "NoGeometry"-Flag and the setSubsetOfAttributes-method as it was suggested here. But it seems to be useless. What am I doing wrong here? How can I speed up the process of searching for the nodeids in the nodes-layer?
i = 0 n = len(feats) # loop through graph layer and find node coordinates for from/to nodes for i, feat1 in enumerate(feats): progress.setPercentage(int(100 * i / n)) fromNodeId = feat1.attributes()[idx_from] toNodeId = feat1.attributes()[idx_to] # search for fromNode by fromNodeId exp = QgsExpression('nodeid = \'' + fromNodeId +'\' ') request = QgsFeatureRequest(exp) request.setFlags(QgsFeatureRequest.NoGeometry) request.setSubsetOfAttributes(['nodeid', 'x', 'y'], inputNodes.fields()) request.setLimit(1) print "time 1" for feat2 in inputNodes.getFeatures(request): from_x = feat2.attributes()[idx_x] from_y = feat2.attributes()[idx_y] print "time 2" # search for specific ToNode by ToNodeId exp = QgsExpression('nodeid = \'' + toNodeId +'\' ') request = QgsFeatureRequest(exp) request.setFlags(QgsFeatureRequest.NoGeometry) request.setSubsetOfAttributes(['nodeid', 'x', 'y'], inputNodes.fields()) request.setLimit(1) print "time 3" for feat2 in inputNodes.getFeatures(request): to_x = feat2.attributes()[idx_x] to_y = feat2.attributes()[idx_y] print "time 4" fromPoint = QgsPoint(from_x, from_y) toPoint = QgsPoint(to_x, to_y) newLine = QgsFeature(outLayer.pendingFields()) (...) azimuth = fromPoint.azimuth(toPoint) newLine.setAttribute(AZIMUTH, azimuth) (...)
EDIT: I have added some timing print outputs. From print "time 1" to print "time 2" it takes approximately 4 seconds. The same with print "time 3" to print "time 4". This means one runthrough the outher for loop takes approximately 8 seconds (I know measuring the time with prints isn't the best way. But I just want to show the bottlenecks of the python code).
If you need some more information on my question, please feel free to leave a comment down below!