3

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:

The graph layer holds 177350 line features which are representing a streetmap. Loaded into QGIS it looks like this:

enter image description here

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:

enter image description here

The "nodes"-layer:

The nodes layer holds 87667 point features which are representing the start- or endnodes of the streets in the graph-layer. Loaded into QGIS and overlayed with the graph-layer it looks like this: enter image description here

Every node in the point-layer has a unique nodeid, x-, y- coordinates and a type (the type attribute is not relevant):

enter image description here

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:

enter image description here

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!

  • How long does it take – BERA Feb 8 at 11:50
  • 1
    16 days approximately – applebrown Feb 8 at 11:51
  • With my coordinates from the start- and endnode, I want to calculate the azimuth. I will add the next steps to the question body above. please have a look! – applebrown Feb 8 at 11:59
  • what exactly do you mean? should I describe the python code snippet more in detail? – applebrown Feb 8 at 12:16
  • @BERA I have now described the issue more in detail. Maybe you can look through it again? – applebrown Feb 8 at 12:41
2

A few different ways to do this but if your lines are drawn in the same order as the points you can do this which avoids most of the API and just uses expressions:

features = iface.activeLayer().getFeatures()
exp = QgsExpression("azimuth(start_point($geometry), end_point($geometry))")

for f in features:
    result = exp.evaluate(f)
    print(result)
    ...
    newLine['azimuth'] = result
  • the values I am getting from this azimuth calculation are very small (between 0 - 5?). Could you please explain to me what this azimuth is? What I calculated in my python script was the north-bounded azimuth from which I got values between -180 to +180. – applebrown Feb 11 at 14:10

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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