4

I have a aggregate process that takes about 9h to run collecting data from a point dataset to put into a line.

The process is

aggregate(layer:='Roads_10m_Chainage.tab', aggregate:='max', expression:="distance" * (attribute(@parent, 'RoadMntnc') = attribute($currentfeature, 'RoadMntnc')),filter:= within(start_point(geometry(@parent)),buffer($geometry,5,10)))

This ensures the value from the point data has the same RoadMntnc value as the line it's looking at the start point of. The buffer ensures it collects the closest value -this works as the chainage is 10m so there will always be a point within 5m.

I was hoping to speed up the process by using a PyQGIS selection of the road line and chainage point datasets having the same RoadMntnc value rather than it crunching through the whole dataset.

So far I have

from operator import itemgetter

def Update_Field(fieldName, ntext): # Function to simplify updating fields
    print (ntext)
    fieldIndex = layer.fields().indexFromName(fieldName)
    dpr.changeAttributeValues({feature.id(): {fieldIndex: ntext}})
    layer.commitChanges()

roadlist=[]
layer = iface.activeLayer()
dpr = layer.dataProvider()

chainage_layer=QgsVectorLayer(r'Z:\My Drive\Mangoesmapping\Spatial Projects\2019\DSC\126_Road_Updates\Working\Output\DSC_Roads_10m_Chainage.tab', "Chainage", "ogr")
#for feature in layer.getFeatures(): # Create list of unique roads to work through
    #if feature['RoadMntnc'] not in roadlist:
    #    roadlist.append(feature['RoadMntnc'])
roadlist=['Kingston_Rd'] #Manual road list for testing
for road in roadlist: # for each road get key attributes and sort
    print (road)
    layer.selectByExpression('"RoadMntnc"= \'{0}\''.format(road), QgsVectorLayer.SetSelection)
    selection=layer.selectedFeatures()
    l_selection=layer.selectedFeatures()
    i=1
    road_parts=[]
    chainage_layer.selectByExpression('"RoadMntnc"= \'{0}\''.format(road), QgsVectorLayer.SetSelection)
    chainage_selection=chainage_layer.selectedFeatures()
    for f in selection:
        for l in chainage_selection:
            if f.geometry().intersects(l.geometry()): #how do I check if a 1m buffer of the start point has a chainge point within it?
                print (l.id(),l['RoadMntnc'], l['distance'])
#if the above works then update record
                #ntext=str(item['Distance'])
                #print (ntext)
                #fieldName='Start_Chainage'
                #Update_Field(fieldName,ntext)

All values are printed in the last line - how to I change the if statement to just show the distance value from the chainage value within 1m of the start of the line?

How do I write the expression to utilize the two selections or modify the aggregate function to use the two selections or use an envelope to collect the required value?

UPDATE: Using César Arquero answer

The buffer to check needs to be around the start point of each road segment. All I need is the chainage feature that intersects this buffer - it should just be one point, but if more than one I need the lowest value.

from operator import itemgetter

def Update_Field(fieldName, ntext): # Function to simplify updating fields
    print (ntext)
    fieldIndex = layer.fields().indexFromName(fieldName)
    dpr.changeAttributeValues({feature.id(): {fieldIndex: ntext}})
    layer.commitChanges()

roadlist=[]
layer = QgsProject.instance().mapLayersByName('DSC_Roads')[0]
chainage_layer= QgsProject.instance().mapLayersByName('DSC_Roads_10m_Chainage')[0]

for feature in layer.getFeatures(): # Create list of unique roads to work through
    if feature['RoadMntnc'] not in roadlist:
        roadlist.append(feature['RoadMntnc'])
#
roadlist=['Dean_Rd'] #for testing
#-----------------------------------------------------------------------
#OPTION 1
for road in roadlist: # for each road get key attributes and sort
    i=1
    road_parts=[]
    request = QgsFeatureRequest(QgsExpression('"RoadMntnc"= \'{0}\''.format(road)))
    l_selection = layer.getFeatures(request)
    chainage_selection = chainage_layer.getFeatures(request)

    for l in l_selection:
        geom = l.geometry().asPolyline()
        print ("Start: " + str(geom[0]))
        print ("End: " + str(geom[-1]))
        for c in chainage_selection:
            buffer = geom[0].buffer(5,5) #5m buffer with 5 segments
            if c.geometry().intersects(buffer):
                print(l.id(),c['RoadMntnc'],c['distance'])
#if the above works then update record
#                ntext=str(l['distance'])
#                print(ntext)
#                fieldName='Start_Chainage'
#                Update_Field(fieldName,ntext)

UPDATE 2 For some reason it's only updating roads that have one section and not multiple and also where there more than one point within the buffer it's not putting in any value. I have cleaned up the code to try and see why it's not functioning as expected...

def Check_Intersection (c_sel, l, start_end): #Find points within start/end of line
    geom = l.geometry().asPolyline()
    g = QgsGeometry.fromPointXY(geom[start_end])
    buffer = g.buffer(10,4) #10m buffer with 4 segments
    #buffer_to_iface(buffer) #add buffer to map for testing
    for c in c_sel: # for each chainage point -ideally should be envelope in buffer area only
        #print (c['RoadMntnc'],c['distance']) #show value for testing
        if c.geometry().within(buffer):
            print(l.id(),c['RoadMntnc'],int(c['distance']))
            #Update record
#                layer.startEditing()
#                ntext=int(c['distance'])
#                fieldName='Start_Chainage2'
#                FieldIndex = l.fields().indexFromName(fieldName)
#                layer.changeAttributeValue(l.id(), FieldIndex, ntext)
#                layer.commitChanges()
            break
def  buffer_to_iface(buffer):       #buffer to iface for testing
    buff_layer = QgsVectorLayer('Polygon', 'buffer' , 'memory')
    prov = buff_layer.dataProvider()
    feat = QgsFeature()
    feat.setGeometry(buffer)
    prov.addFeatures([feat])
    buff_layer.updateExtents()
    QgsProject.instance().addMapLayers([buff_layer])[0]

roadlist=[]
layer = QgsProject.instance().mapLayersByName('DSC_Roads')[0]
chainage_layer= QgsProject.instance().mapLayersByName('DSC_Roads_10m_Chainage')[0]


for feature in layer.getFeatures(): # Create list of unique roads to work through
    if feature['RoadMntnc'] not in roadlist:
        roadlist.append(feature['RoadMntnc'])
#
#roadlist=['Kingston_Rd'] #for testing
#-----------------------------------------------------------------------

for road in roadlist: # for each road select line and point features
    request = QgsFeatureRequest(QgsExpression('"RoadMntnc"= \'{0}\''.format(road)))
    l_sel = layer.getFeatures(request)
    c_sel = chainage_layer.getFeatures(request)

    for l in l_sel: # repeat for each road segment
        print (str(l.id())+l['RoadMntnc'])
        start_end=0 #get start points
        Check_Intersection (c_sel,l,start_end)

print ('Completed')

layer.updateFields()
2
  • I don't see clearly where you want to go exactly. Could you share your layers to try the exact aggregate function and see how it goes ? Commented Sep 13, 2019 at 11:44
  • @CésarArquero -sample data is in drive.google.com/open?id=1RL7VZL_giZuZ9xVfDFgk409_RIErIB0m - Basically I am just trying to find a quicker process to get point data attributes based on the start/end of a line dataset where both have a common attribute value (a road name in this case).
    – GeorgeC
    Commented Sep 13, 2019 at 12:11

1 Answer 1

3
+50

I added two options to your code:

  1. Just added a buffer to the line in order to get the attributes of the intersected point (I still have doubts about the path you are trying).
  2. This second option is based on the aggregated process from your 1st line of code. In this case, the point layer is filtered by road name and then, aggregation is performed.

    from operator import itemgetter
    
    def Update_Field(fieldName, ntext): # Function to simplify updating fields
        print (ntext)
        fieldIndex = layer.fields().indexFromName(fieldName)
        dpr.changeAttributeValues({feature.id(): {fieldIndex: ntext}})
        layer.commitChanges()
    
    roadlist=[]
    layer = QgsProject.instance().mapLayersByName('DSC_Roads_10m_Chainage')[0]
    chainage_layer= QgsProject.instance().mapLayersByName('DSC_Roads')[0]
    
    for feature in layer.getFeatures(): # Create list of unique roads to work through
        if feature['RoadMntnc'] not in roadlist:
            roadlist.append(feature['RoadMntnc'])
    
    #-----------------------------------------------------------------------
    #OPTION 1
    for road in roadlist: # for each road get key attributes and sort
        i=1
    road_parts=[]
    request = QgsFeatureRequest(QgsExpression('"RoadMntnc"= \'{0}\''.format(road)))
    l_selection = layer.getFeatures(request)
    chainage_selection = chainage_layer.getFeatures(request)
    for l in l_selection:
        for c in chainage_selection:
            buffer = l.geometry().buffer(1,20) #1m buffer with 20 segments
            if l.geometry().intersects(buffer):
                print(l.id(),l['RoadMntnc'],l['distance'])
    #if the above works then update record
    #                ntext=str(l['distance'])
    #                print(ntext)
    #                fieldName='Start_Chainage'
    #                Update_Field(fieldName,ntext)
    
    #-----------------------------------------------------------------------
    #OPTION 2
    Alist=[]
    for road in roadlist: # for each road get key attributes and sort
        print(road)
        layer.setSubsetString('') #remove any filter on the layer
        layer.setSubsetString('"RoadMntnc"= \'{0}\''.format(road))
        A = processing.run("qgis:aggregate", {
            'INPUT':layer,
            'GROUP_BY':'\"distance\" * (attribute(@parent, \'RoadMntnc\') = attribute($currentfeature, \'RoadMntnc\'))',
            'AGGREGATES':[
                {'aggregate': 'concatenate', 'delimiter': ',', 'input': '"RoadMntnc"', 'length': 100, 'name': 'RoadMntnc', 'precision': 0, 'type': 10},
                {'aggregate': 'concatenate', 'delimiter': ',', 'input': '"RoadName"', 'length': 50, 'name': 'RoadName', 'precision': 0, 'type': 10}, 
                {'aggregate': 'concatenate', 'delimiter': ',', 'input': '"Locality"', 'length': 30, 'name': 'Locality', 'precision': 0, 'type': 10},
                {'aggregate': 'concatenate', 'delimiter': ',', 'input': '"Rd_Control"', 'length': 30, 'name': 'Rd_Control', 'precision': 0, 'type': 10},
                {'aggregate': 'sum', 'delimiter': ',', 'input': '"Chainage"', 'length': 10, 'name': 'Chainage', 'precision': 0, 'type': 6}, 
                {'aggregate': 'sum', 'delimiter': ',', 'input': '"GIS_Length"', 'length': 10, 'name': 'GIS_Length', 'precision': 0, 'type': 6},
                {'aggregate': 'sum', 'delimiter': ',', 'input': '"Segments"', 'length': 10, 'name': 'Segments', 'precision': 0, 'type': 6},
                {'aggregate': 'concatenate', 'delimiter': ',', 'input': '"Publish"', 'length': 10, 'name': 'Publish', 'precision': 0, 'type': 10},
                {'aggregate': 'concatenate', 'delimiter': ',', 'input': '"Comments"', 'length': 254, 'name': 'Comments', 'precision': 0, 'type': 10},
                {'aggregate': 'concatenate', 'delimiter': ',', 'input': '"Org_Ref"', 'length': 20, 'name': 'Org_Ref', 'precision': 0, 'type': 10}, 
                {'aggregate': 'concatenate', 'delimiter': ',', 'input': '"Rd_Id"', 'length': 10, 'name': 'Rd_Id', 'precision': 0, 'type': 10},
                {'aggregate': 'sum', 'delimiter': ',', 'input': '"QRA_Ref"', 'length': 10, 'name': 'QRA_Ref', 'precision': 0, 'type': 6},
                {'aggregate': 'maximum', 'delimiter': ',', 'input': '"distance"', 'length': 0, 'name': 'distance', 'precision': 0, 'type': 6},
                {'aggregate': 'sum', 'delimiter': ',', 'input': '"angle"', 'length': 0, 'name': 'angle', 'precision': 0, 'type': 6}],
             'OUTPUT':'TEMPORARY_OUTPUT'})
        Alist.append(A['OUTPUT']) #add the output to the list
    
    #Remove filters and merge all layers
    layer.setSubsetString('') #remove the fiter
    Amerged = processing.run("native:mergevectorlayers", {
        'LAYERS':Alist,
        'CRS':QgsCoordinateReferenceSystem('EPSG:25833'),
        'OUTPUT':'TEMPORARY_OUTPUT'})
    
    #OUTPUT to iface
    QgsProject.instance().addMapLayers([Amerged['OUTPUT']])
    

EDIT of option 1: buffer from start point of a line:

from operator import itemgetter
def Update_Field(fieldName, ntext): # Function to simplify updating fields
    print (ntext)
    fieldIndex = layer.fields().indexFromName(fieldName)
    dpr.changeAttributeValues({feature.id(): {fieldIndex: ntext}})
    layer.commitChanges()

roadlist=[]
layer = QgsProject.instance().mapLayersByName('DSC_Roads')[0]
chainage_layer= QgsProject.instance().mapLayersByName('DSC_Roads_10m_Chainage')[0]

for feature in layer.getFeatures(): # Create list of unique roads to work through
    if feature['RoadMntnc'] not in roadlist:
        roadlist.append(feature['RoadMntnc'])
#
roadlist=['Schilds_Rd'] #for testing
#-----------------------------------------------------------------------
#OPTION 1
for road in roadlist: # for each road get key attributes and sort
    i=1
    road_parts=[]
    request = QgsFeatureRequest(QgsExpression('"RoadMntnc"= \'{0}\''.format(road)))
    l_selection = layer.getFeatures(request)
    chainage_selection = chainage_layer.getFeatures(request)

    for l in l_selection:
        geom = l.geometry().asPolyline()
        print ("Start: " + str(geom[0]))
        print ("End: " + str(geom[-1]))
        for c in chainage_selection:
            g = QgsGeometry.fromPointXY(geom[0])
            buffer = g.buffer(5,5) #5m buffer with 5 segments

            #buffer to iface
            layer = QgsVectorLayer('Polygon', 'buffer' , 'memory')
            prov = layer.dataProvider()
            feat = QgsFeature()
            feat.setGeometry(buffer)
            prov.addFeatures([feat])
            layer.updateExtents()
            QgsProject.instance().addMapLayers([layer])

            if c.geometry().intersects(buffer):
                print(l.id(),c['RoadMntnc'],c['distance'])
#if the above works then update record
#                ntext=str(l['distance'])
#                print(ntext)
#                fieldName='Start_Chainage'
#                Update_Field(fieldName,ntext)
6
  • Thanks -with regards to Option 1 - the buffer just needs to be at the start of the each segment of the line to find the closed chainage point with the same road name. I have been trying to change it based on gis.stackexchange.com/a/304805/2891 but haven't been able to get it to work.
    – GeorgeC
    Commented Sep 16, 2019 at 4:39
  • do you know how to modify your suggestion to get the value from the point closest to the start point? - I have updated my code based on yours.
    – GeorgeC
    Commented Sep 17, 2019 at 4:29
  • 1
    You need to get a geometry from a QgsPointXY (you got it here: docs.qgis.org/3.4/en/docs/pyqgis_developer_cookbook/…). See the code above to check the update. Commented Sep 17, 2019 at 7:37
  • Thanks -that works well in most cases and seems a lot quicker. I want to add the buffer that's created to to the map for testing why it doesn't work sometimes - I have tried buffer_obj=buffer.geometry.asPolygon()' and QgsProject.instance().addMapLayer(buffer)` but both give errors.
    – GeorgeC
    Commented Sep 17, 2019 at 11:35
  • I've added that part to the code. May be you find a better solution in order to reduce the number of lines. Regards ¡ Commented Sep 17, 2019 at 21:46

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.