PyQGIS – Populating Line Attribute Fields with Attributes from Point Features

field-calculatorpyqgis

I have a line layer (MultiLineString) and 4 point layers (MultiPoint).

My line layer has 2 attribute fields, start_pt and end_pt, I need them populated with an identifier attribute (string data populated with values like A01, B01, B02, C01, etc. so it's not the regular $id value) from any of the potential 4 point features that can be snapped to the line feature.

Here's a sample project:

This is what line layer looks like:

Here's what the populated attribute fields should be like: (filled manually)

I'm looking for a mean to do this through either field calculator expressions OR PyQGIS as I need it as part of a plugin to automate it. So I can't use DB/SQL or anything else, basically.

What I have used so far was a calculator expression (tried changing start_point/end_point and array_first and array_last):

if (
within( array_first (overlay_touches( 'point1',$geometry)), buffer (start_point ($geometry), 0.00001)),
    array_first (overlay_touches( 'point1',"id")),
    array_last (overlay_touches( 'point2',"id")))

Which I would then assign to a variable with QgsExpression() and populate my tables with context.SetFeature(feature)

Back when I was limited to 2 point layers it used to work with relative success but as of right now it doesn't do the trick mainly because within() in and of itself is limited to 2 values for either True or False so it doesn't work for more than 2 layers.

I'm still a beginner with Python and PyQGIS so I'm not sure how to make it work with it either. I know I could use a spatial index as a mean to detect closes features from a given point layer that are close to my line feature's start_point or end_point (as it is shown at pyqgis: add attributes of points to attributtable of lines) but as I understand, only one spatial index is used at a time, so I'm not sure about that either.

How can I, for instance, retrieve the "id" attribute from the retrieved nearestNeighbor from a spatial index?

In the thread linked about, there's this bit of code which theoretically explains how to retrieve the attribute id of the given feature. But it obviously only works with a single layer in mind.

fs =  point_layer.getFeatures(QgsFeatureRequest().setFilterFid(nearest_to_start))
point_feature = QgsFeature()
fs.nextFeature(point_feature)
p_id = point_feature['id']

Is there a way to use multiple spatial indices to somehow achieve my goal?

I'm updating this question as the solution provided in the edit by @pigreco doesn't really solve the problem. Using his help I got on the right track. This is more of an elabored comment than an answer as I do need to provide screenshots and code.

Here's an example following his method:

Overlay_nearest returns wrong values whenever two points are from the same layer. I have tried setting the limit to 1 and max_distance to 0 but the results remain sensibly the same. Replacing the first overlay_nearest by overlay_intersects (which ensures only the interescting geometry gets returned) does limit the amount of values (thus getting rid of cases like A03 replacing D04) instead I get NULL values instead of wrong values.

Here's the code I used:

-- select id
with_variable('feature',
-- search for the nearest points
overlay_intersects(layer:=
-- search for the closest layer
with_variable('in_layer',array('point1','point2','point3','point4'), -- point layer list
expression:=with_variable('in_dist',
array_foreach(@in_layer,
distance(overlay_nearest(@element,$geometry)[0],
start_point($geometry))), -- change start or end_point
array_get(@in_layer, array_find(@in_dist, array_sort(@in_dist)[0]))))
-- search for the closest layer
,expression:= id )
-- search for the nearest points
, if(array_length(
@feature)>1,
@feature[0], -- 0: start_point; 1: end_point
@feature[0]) --
)

Here are the results I get:

I have also tried with an aggregate expression:

aggregate(layer:=with_variable('in_layer',array('point1','point2','point3','point4'),
          expression:=with_variable('in_dist',
          array_foreach(@in_layer, distance(overlay_nearest(@element,$geometry)[0], end_point($geometry))), 
                 array_get(@in_layer, array_find(@in_dist, array_sort(@in_dist)[0])))),
          aggregate:='max', expression:=id, 
          filter:=intersects( $geometry, end_point(geometry(@parent) )))

Which gives me these results:

Although and to be quite honest, I'm fairly confused at what is returned by the array_forreach(array_get()) and the relation with my intersects filter which takes the current feature's geometry.

Best Answer

Using the QGIS (>= 3.18) field calculator:

for the start_pt field:

overlay_nearest(layer:=
with_variable('in_layer',array('cat_A','cat_B','cat_C','cat_D'),
with_variable('in_dist',
        array_foreach(@in_layer,distance(overlay_nearest(@element,$geometry)[0], start_point($geometry))),
array_get(@in_layer, array_find(@in_dist, array_min(@in_dist)))))
,expression:= id )[0]

for the end_pt field:

overlay_nearest(layer:=
with_variable('in_layer',array('cat_A','cat_B','cat_C','cat_D'),
with_variable('in_dist',
        array_foreach(@in_layer,distance(overlay_nearest(@element,$geometry)[0], end_point($geometry))),
array_get(@in_layer, array_find(@in_dist, array_min(@in_dist)))))
,expression:= id )[0]

The previous expressions work for the case described, but if a linear element had points belonging to the same vector point, the expressions would no longer work, as it would not distinguish the start_point from the end_point;

I add the expression that solves the problem:

-- select id
with_variable('feature',
-- search for the nearest points
    overlay_nearest(layer:=
    -- search for the closest layer
        with_variable('in_layer',array('cat_A','cat_B','cat_C','cat_D'), -- point layer list
            with_variable('in_dist',
                array_foreach(@in_layer,
                       distance(overlay_nearest(@element,$geometry)[0],
                end_point($geometry))), -- change start or end_point
        array_get(@in_layer, array_find(@in_dist, array_sort(@in_dist)[0]))))
    -- search for the closest layer
        ,expression:= id, limit:=2, max_distance:=0.1 )
-- search for the nearest points
, if(array_length(
        @feature)>1,
        @feature[1], -- 0: start_point; 1: end_point
        @feature[0]) -- 
)
-- select id

Related Solutions

[GIS] Model Builder Error populating attributes with Calculate Field

Ok, the solution was very simple. Adding double quotes around the %Value% and/or %Name% variable lets model builder recognize the iterator output as an inline variable. So the Calculate field expression that worked was: "%Value%"

QGIS PyQGIS: Snapping Lines’ Start/End Points to Nearest Features

Here is a solution for MultiLineStrings and MultiPoints. It moves the start- and endpoints of every part of the MultiLineString to the closest point.

In case you do not want to move start- and enpoints of every part, but only of every feature, you may want to change these two lines accordingly (simply comment/uncomment in the code below):

if part.startPoint() == geom.vertexAt(vert): to if feat.geometry().vertexAt(0) == geom.vertexAt(vert):
elif part.endPoint() == geom.vertexAt(vert): to elif feat.geometry().vertexAt([i for i, f in enumerate(feat.geometry().vertices())][-1]) == geom.vertexAt(vert): (not very convenient, but simply doing feat.geometry().vertexAt(-1) directly doesnt work...)

Here is the full code with further explanations as comments:

points = QgsProject.instance().mapLayersByName('points')[0]
lines = QgsProject.instance().mapLayersByName('lines')[0]
tolerance = 10000 # snapping tolerance in CRS units

if QgsWkbTypes.isMultiType(points.wkbType()): # if point-input is of type multigeometry
    res = processing.run("native:multiparttosingleparts",{'INPUT':points,'OUTPUT':'TEMPORARY_OUTPUT'}) # convert to singlepoint
    singlepoints = res['OUTPUT']
else: # if point-input is of type singlegeometry 
    singlepoints = points # leave the inputlayer as it is

# create the spatial index of the pointlayer
points_idx = QgsSpatialIndex(singlepoints.getFeatures())

with edit(lines):
    for feat in lines.getFeatures():
        geom = feat.geometry()
        vert = 0 # reset vertices-count on every new feature
        vert_dict = {} # clear vertices dictionary on every new feature
        for part in geom.parts(): # iterate through all parts of each feature
            for vertex in part.vertices(): # iterate through all vertices of each part
                # if feat.geometry().vertexAt(0) == geom.vertexAt(vert): # if the current vertex is a startpoint of the feature
                if part.startPoint() == geom.vertexAt(vert): # if the current vertex is a startpoint of the features part, then:
                    nearestneighbor_start = points_idx.nearestNeighbor(QgsPointXY(part.startPoint()), neighbors=1, maxDistance=tolerance) # find the closest point to the current startpoint
                    if len(nearestneighbor_start) > 0:# only adjust the startpoint if there is a nearest point within the maxdistance
                        nearpoint_startgeom = singlepoints.getFeature(nearestneighbor_start[0]).geometry() # get the geometry of the nearest point in pointlayer
                        vert_dict[vert] = nearpoint_startgeom.asPoint() # add the index of the current vertex as key and the geometry of the closest point as value to the dictionary
                # elif feat.geometry().vertexAt([i for i, f in enumerate(feat.geometry().vertices())][-1]) == geom.vertexAt(vert): # if the current vertex is an endpoint of the current feature
                elif part.endPoint() == geom.vertexAt(vert): # if the current vertex is an endpoint of the current features part
                    nearestneighbor_end = points_idx.nearestNeighbor(QgsPointXY(part.endPoint()), neighbors=1, maxDistance=tolerance)
                    if len(nearestneighbor_end) > 0:
                        nearpoint_endgeom = singlepoints.getFeature(nearestneighbor_end[0]).geometry()
                        vert_dict[vert] = nearpoint_endgeom.asPoint()
                else:
                    pass # if current vertex is not a start or end point skip it...
                vert += 1 # increase vertices-counter
        for vertindex, newpoint in vert_dict.items(): # for every feature iterate over the just created dictionary (vertindex (=dict key) is the start or endpoint we want to move and newpoint (=dict value) the position we want to move it to)
            QgsVectorLayerEditUtils(lines).moveVertex(newpoint.x(),newpoint.y(),feat.id(),vertindex) # use QgsVectorLayerEditUtils to edit the vertex, for more details see https://qgis.org/pyqgis/3.4/core/QgsVectorLayerEditUtils.html#qgis.core.QgsVectorLayerEditUtils
# https://gis.stackexchange.com/questions/410768/snapping-lines-start-point-end-points-to-the-closest-point-features-to-them/410778#410778

For further background you may want to read Modifying specific vertices of MultiLineString using PyQGIS

Example for snapping start- and endpoints of all parts of all feature:

Example for snapping start- and enpoints only for all features:

Best Answer

Related Solutions

[GIS] Model Builder Error populating attributes with Calculate Field

QGIS PyQGIS: Snapping Lines’ Start/End Points to Nearest Features

Related Question