3

I have one line-layer which pictures a street net and a point-layer which defines certain locations. Please see in the picture below: enter image description here

What I am trying to do is to create a normal vector through a point and its nearest line for each point. Then I would like to calculate the distance between the line's starting point and the intersection of the normalvector. This should look like this: enter image description here

As I have lots of points, I need to find a way to solve this generically.

I know that the plugin NNJoin gives me the nearest edge to a point, but how do I create normal vectors and how can I then calculate the lengths?

Do I have to write a python script to realize this or is there an other way?

  • Ouch, it is closed now... let me suggest v.net.connect, it is found in the QGIS Processing Toolbox | GRASS 7 | Vector. (1) You will obtain new Network line layer. (2) Then you can select the first segment which represents Start - Intersection interval. (3) Save them as new layer. (4) Calculate $length. – Kazuhito Feb 19 '18 at 20:49
  • I tried it with the v.net.connect function, I keep getting the error: Vector map not found although I updated my QGIS version to QGIS 2.18.16. – applebrown Feb 20 '18 at 11:17
  • Does your shapefile have any blank space in its name, by any chance? – Kazuhito Feb 20 '18 at 11:21
  • 1
    This sounds like a job for pgRouting, have you explored this? docs.pgrouting.org/2.5/en/… – cm1 Feb 28 '18 at 17:10
  • 2
    It is a matter of efficiency, PostGIS is one of the many tools avaivable in QGIS and defenetively a powerful. Spatial databases are very suited for your task and we should consider them as part of the QGIS features – Marco Mar 5 '18 at 12:34
3
+50

Here is Python solution (for QGIS Python console) which is mostly created with the help of the following links:

  1. Create Intersecting lines and remove dangles
  2. How do I find vector line bearing in QGIS or GRASS?
  3. How to create points in a specified distance along the line in QGIS?
  4. Redirecting stdout to “nothing” in python
  5. QGIS 2.18 Network analysis library

Following steps are involved:

  1. Find point on nearest road to given point and create a linestring.
  2. Extend the linestring by a small distance to make sure the linestring is intersecting the road. (I recognized that there are less intersections without using the extended linestring.)
  3. Finding the shortest path from the intersection point to the linestring start point using dijkstra() method from the QGIS network analysis library.

Test data:

  1. Road layer with 7220 features:

    enter image description here

  2. Random points layer with 7220 features (using the Toolbox: Bounding boxes from road layer, Random points in layer bounds (with minimum distance = 10))

    enter image description here

Python code:

import math
import os
import sys
from contextlib import contextmanager
from PyQt4.QtCore import *
from PyQt4.QtGui import *

from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *

# define input layer points and roads
p_lyr = QgsMapLayerRegistry.instance().mapLayersByName('points')[0]
l_lyr = QgsMapLayerRegistry.instance().mapLayersByName('roads')[0]

# function to create the azimuth of line segments   
def azimuth(point1, point2):
   return point1.azimuth(point2)

# project a point in azimuth direction using direction cosines
def cosdir_azim(azim):
   azim = math.radians(azim)
   cosa = math.sin(azim)
   cosb = math.cos(azim)
   return cosa,cosb

# used to get line segments
def pair(list):
    for i in range(1, len(list)):
        yield list[i-1], list[i]

# create function to hide QGIS output ("True") while creating QgsSpatialIndex for every feature (takes a while for x*1000 features)
@contextmanager
def silence_stdout():
    new_target = open(os.devnull, "w")
    old_target, sys.stdout = sys.stdout, new_target
    try:
        yield new_target
    finally:
        sys.stdout = old_target

# creating index silently
with silence_stdout():
    lines = [feature for feature in l_lyr.getFeatures()]
    lines_spIndex = QgsSpatialIndex()
    for elem in lines:
        lines_spIndex.insertFeature(elem)

# set up memory layer for the normal vector
d_lyr = QgsVectorLayer('LineString', 'normalVector', 'memory')
QgsMapLayerRegistry.instance().addMapLayer(d_lyr)
prov = d_lyr.dataProvider()
# adding three attributes (holding point_id, road_id and the distance)
prov.addAttributes( [ QgsField("point_id", QVariant.Int), QgsField("road_id", QVariant.Int), QgsField("distance",QVariant.Int)])

# function to create normal vector (geometry) and distance
def normalVector(distance):
    vect = []
    feat = []
    for points in p_lyr.getFeatures():
        # find closest point to line
        minDistPoint = min([l.geometry().closestSegmentWithContext(QgsPoint(points.geometry().asPoint())) for l in lines])[1]
        feat = QgsFeature()
        # create line from point to minDistPoint
        feat.setGeometry(QgsGeometry.fromPolyline([QgsPoint(points.geometry().asPoint()), QgsPoint(minDistPoint[0], minDistPoint[1])]))
        line = feat.geometry().asPolyline()
        for seg_start, seg_end in pair(line):
            cosa, cosb = cosdir_azim(azimuth(seg_start, seg_end))
            lenght = distance
            feat = QgsFeature()
            # extend line for intersection with road feature
            feat.setGeometry(QgsGeometry.fromPolyline([seg_start,QgsPoint(seg_end.x()+(lenght*cosa), seg_end.y()+(lenght*cosb))]))
            for id in lines_spIndex.intersects(feat.geometry().boundingBox()):
                sp_geom = lines[id].geometry().asPolyline()
                start_point = QgsPoint(sp_geom[0])
                sp_x, sp_y = start_point[0], start_point[1]
                if feat.geometry().intersects(lines[id].geometry()):
                    intersection_point = feat.geometry().intersection(lines[id].geometry()).asPoint()
                    end_point = QgsPoint(intersection_point)
                    ep_x, ep_y = end_point[0], end_point[1]
                    # building graph and calculating shortest path
                    director = QgsLineVectorLayerDirector(l_lyr, -1, '', '', '', 3)
                    properter = QgsDistanceArcProperter()
                    director.addProperter(properter)
                    crs = qgis.utils.iface.mapCanvas().mapRenderer().destinationCrs()
                    builder = QgsGraphBuilder(crs)
                    pStart = QgsPoint(ep_x, ep_y )
                    pStop = QgsPoint(sp_x, sp_y)
                    tiedPoints = director.makeGraph(builder, [pStart, pStop])
                    graph = builder.graph()
                    tStart = tiedPoints[0]
                    tStop = tiedPoints[1]
                    idStart = graph.findVertex(tStart)
                    idStop = graph.findVertex(tStop)
                    (tree, cost) = QgsGraphAnalyzer.dijkstra(graph, idStart, 0)
                    if tree[idStop] == -1:
                        print "Path not found"
                    else:
                        p = []
                        curPos = idStop
                        while curPos != idStart:
                            p.append(graph.vertex(graph.arc(tree[curPos]).inVertex()).point())
                            curPos = graph.arc(tree[curPos]).outVertex()
                        p.append(tStart)
                    vect = QgsFeature()
                    vect.setGeometry(QgsGeometry.fromPolyline(p))
                    vect.setAttributes([int(points["id"]), int(lines[id]["id"]), vect.geometry().length()])
        prov.addFeatures([vect])
    d_lyr.updateExtents()
    d_lyr.triggerRepaint()
    d_lyr.updateFields()

# run function
normalVector(0.1)

Results:

For the amount of data as described above it takes about 16 minutes of processing time. For 11 times I got Path not found.

enter image description here

  • 1
    definetely better speed! one thing: this will return the distance based on a straight line between intersection point and start point...could be important to know (imagine a serpentine road: the straight distance between intersection and start point could be some 20m, while the road would actually wind up that hill for 5km) – ThingumaBob Mar 7 '18 at 12:40
  • 1
    @applebrown - I've edited my answer. In addition to the length the resulting normalVector layer holds also the point id. In the part vect.setAttributes([int(p["ogc_fid"]), vect.geometry().length()]) you have to change "ogc_fid" to the your point id field name (hopefully there is one). – Stefan Mar 7 '18 at 12:42
  • 1
    To check the duplicate points (there should be the problem with majority of distances) create a virtual layer: SELECT id, count(*) AS count FROM normalVector GROUP BY 1 ORDER BY count. Click No Geometry. The reason could be the extended line which intersects more than one linestring (when the linestring are really close together). You could decrease the distance (e.g. to 0.001). – Stefan Mar 7 '18 at 12:45
  • 1
    @ThingumaBob Ok. applebrown looks for the length of the piece of road which is between the intersection point and the start point. I misinterpret "[] ... distance between the line's starting point and the intersection of the normalvector." I will change this. – Stefan Mar 7 '18 at 12:50
  • 1
    @ThingumaBob, applebrown - I've edited my answer. Now, it works not only for straight lines. – Stefan Mar 12 '18 at 9:15
5

I absolutely recommend to check out PostgreSQL/PostGIS, especially if your datasets are large and dynamic!

However, if you want to stick to QGIS for now, try the Virtual Layers. They provide full integration of SpataLite functions to be used on any existing vector layer opened in the project.

The following query will identify the closest line to each point via MIN(ST_Distance(...)) in conjunction with the GROUP BY statement and calculates the distance to the line's start point (in direction of digitizing) with ST_Line_Locate_Point (this function does those two things you are looking for: finding the closest point on the line, i.e. the normal vector intersection, and returning the distance to the start point of the line as a fraction of it's total lenght; thus the part where the ST_Length(...) is multiplied by that fraction):

SELECT sub.pt_id,
       sub.ln_id,
       ST_Length(sub.ln_geom) * ST_Line_Locate_Point(sub.ln_geom, sub.pt_geom) AS length
FROM (
    SELECT p.<id> AS pt_id,
           l.<id> AS ln_id,
           p.geometry AS pt_geom,
           l.geometry AS ln_geom,
           MIN(ST_Distance(p.geometry, l.geometry))
    FROM <points_layer> AS p
      CROSS JOIN <lines_layer> AS l
    GROUP BY p.<id>
) AS sub

This assumes your geometries (i.e. the respective layers) use a projection with meter as unit! Don´t forget to fill in the correct layer and column names.

This returns the id of the point, the id of it's closest line and the distance as the layers attribute table.

Note: the subquery is a rather inefficient construct, yet the only approach I know of, without using the VirtualKNN tables (SpatiaLite starts support of VirtualKNNI with version 4.4.0 and QGIS 2.X and 3.0 come only with SpatiaLite 4.3.0); The whole query might get slow on larger tables.

  • Thank you for your answer, but I would really like to solve this with QGIS. Isn't there any other way with QGIS than with the v.net.connect function? – applebrown Mar 2 '18 at 9:17
  • @applebrown read about QGIS' Virtual Layers (link above)...this is actually a feature in QGIS, and a powerful one. No external DB or tools needed. – ThingumaBob Mar 2 '18 at 10:52
  • What do I have to substitute for sub.pt_id, sub.ln_id, sub.ln_geom and sub.pt_geom? My line layer is called graph and my point layer is called evis. Also, I do not have an unique id in the attributes for both layers. What do I write here? Sorry, I am not that experienced with SQL – applebrown Mar 2 '18 at 13:35
  • 1
    @applebrown you only need to substitute names in <> (i.e. <points_layer> and <line_layer> with evis and graph), all other names are aliases that you can alter as you like. for the p.<id> and l.<id>: QGIS adds a primary key column internally, you can adress that with _uid_ I think (like p._uid_), that should work. – ThingumaBob Mar 2 '18 at 13:43
  • It is calculating pretty slow as you said. I will let you know when I have the output :). – applebrown Mar 2 '18 at 13:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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