Populating line attribute fields with attributes from point features snapped to their start_point or end_point

Question

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.

Answer 1

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

Answer 2

In Python you could:

1. starting on line1, get the first point of the line, there is, the point not connected with the line2, this is the point A1.
2. the point B2 is the point on line1 that has a connection with another line.
3. to decide which line is line2 you have to consider if the line is connected with another or no or, a little better, the line2 is the line with more lines connected after it.
4. once selected which line is line2, the point C1 is the point on line2 as opposed to the point connected on line1
5. to decide which line is line3 you choose the line that has a connection with the lowest label, B01 or C01

You have to implement functions to grasp these situations, that is a significant amount of code, maybe 200 lines, without considering the necessary tests. This will be a good problem to practice code.