Create normal vector and calculate distance between start point and intersection

Question

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

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:

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?

Answer 1

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.

Answer 2

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:

   

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

   

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.