6

I have a point layer of settlements and a LineString layer of a road network. The settlements have been snapped to the network so they are directly connected by it. Settlements are not always at the start or end of the line - they may be halfway along it etc. I would like to output a list of (or otherwise identify) which settlements are directly connected to each other by road.

network of four points

So for example, in the above image we would identify that point A was connected to B and C, but not D. The desired result in full would be the table below.

Point Connected to
A B, C
B A, C, D
C A, B, D
D B, C
4
  • But connected points are alway on the same line? So in your example e.g. between A and B there is always just one single line, not several separate lines?
    – Babel
    Jun 15, 2023 at 9:14
  • "Connected" means what exactly? If between A and B there is another point X on the same line as A and B, then A and B are still considered to be connected or both of them will only be connected to X?
    – Babel
    Jun 15, 2023 at 9:32
  • I guess a different description would be neighbours on the network. So if your point X was added, they would be 'connected' to x and not each other. Just like C interrupts the A and D in the image.
    – Will
    Jun 15, 2023 at 9:58
  • Please, do not forget about "What should I do when someone answers my question?"
    – Taras
    Jul 4, 2023 at 10:46

2 Answers 2

3

You can use QGIS expressions and check for each point which other points are on the same line. For this to work, you have first to make sure that each network line connects only two points: see the screenshot where each separate line has another color and below for how to achieve that.

The label (2nd line) is based on the expression below and automatically generates the name of the neighboring points: enter image description here

To prepare this, proceed as follows:

  1. Create a new (empty) multipart line layer, copy and paste all lines from the network there.

  2. Select and merge all features so that the whole network consists of one single multipart feature.

  3. Split the lines with the points. See here how to do so.

Now you're ready to use the expression on the point layer (replace network in line 13 with the name of the layer created in step 3; and replace name at the end of line 19 and the 4th last line with the name of the attribute field that contains the name of the point):

array_to_string (
    array_distinct(
        array_filter (
            string_to_array(
                array_to_string (
                    array_foreach (
                        generate_series (1, num_geometries(collect($geometry))),
                        with_variable(
                            'iteration',
                            @element,
                            array_to_string (
                                array_foreach(
                                    overlay_intersects ('Split',$geometry),  -- change layer name here
                                    if (
                                        intersects (
                                            buffer (@element,0.1),
                                            geometry(get_feature_by_id(@layer,@iteration))
                                        ),
                                        attribute (get_feature_by_id (@layer,@iteration),'name'), -- change fieldname here
                                        'delete'
                                    )
                                )
                            )
                        )
                    )
                )
            ),
            @element <>'delete' and @element <> name  -- change fieldname here
        )
    )
)
2

If you are open to a Python script:

from itertools import product
from collections import defaultdict
from qgis.core import QgsProject, QgsFeatureRequest

point_layer = QgsProject.instance().mapLayersByName("point_layer")[0]
line_layer = QgsProject.instance().mapLayersByName("line_layer")[0]

dict_result = defaultdict(list)

for line in line_layer.getFeatures():
    request = QgsFeatureRequest()
    request.setDistanceWithin(line.geometry(), 0.001)
    points = list(point_layer.getFeatures(request))
    pt_points = [point.geometry().asPoint() for point in points]
    start = None
    end = None

    for pt in line.geometry().asPolyline():
        try:
            index = pt_points.index(pt)
        except ValueError:
            continue
        if start is None:
            start = points[index]
            continue
        if end is None:
            end = points[index]
            dict_result[start["id_field"]].append(end["id_field"])
            dict_result[end["id_field"]].append(start["id_field"])
            start = end
            end = None

for point_key, point_values in dict_result.items():
    print(point_key, f"[{', '.join(point_values)}]")

enter image description here

result:

A [B]
B [A, C]
C [B, E, D]
E [C, F]
D [C]
F [E, G]
G [F]
2
  • Thanks. This is producing accurate but incomplete results, both in not including all neighbouring points on the network and not completing for all points. The former suggests an issue with the underlying network, though I don't see a systematic difference that explains it at this point. I initially got a "MultiLineString geometry cannot be converted to a polyline." error but resolved that with a simple application of the Multiparts to singleparts function.
    – Will
    Jun 15, 2023 at 10:10
  • This does what you stated, I would double check that you have snapped your points to line vertices and that all lines are not multilinestring with the processing tool
    – Kalak
    Jun 15, 2023 at 10:22

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.