[processing] Fixes for Service Area algorithms

- Output interpolated points when travel cost falls mid-way along
an edge
- Output all intermediate reachable points also
- Make outputting upper/lower bound points optional, and non-default.
Now by default we just output all definitely reachable points and
the interpolated points along edges which correspond to the travel cost.
This allows the output to be used to correctly generate service areas
e.g. by concave/convex polygons and all reachable nodes will be
included in the area.
- Allow algorithm to optionally output a line layer (and make the
point layer optional too, and default to just the line layer output)
containing all reachable line segments (including interpolated
segments of lines when the travel cost sits midway along that
edge). This output is more easily understandably for users.

Author: nyalldawson

python/analysis/network/qgsgraphanalyzer.sip.in

@@ -30,7 +30,8 @@ Solve shortest path problem using Dijkstra algorithm
 :param source: source graph
 :param startVertexIdx: index of the start vertex
 :param criterionNum: index of the optimization strategy
-:param resultTree: array that represents shortest path tree. resultTree[ vertexIndex ] == inboundingArcIndex if vertex reachable, otherwise resultTree[ vertexIndex ] == -1
+:param resultTree: array that represents shortest path tree. resultTree[ vertexIndex ] == inboundingArcIndex if vertex reachable, otherwise resultTree[ vertexIndex ] == -1.
+Note that the startVertexIdx will also have a value of -1 and may need special handling by callers.
 :param resultCost: array of the paths costs
 %End
 

python/plugins/processing/algs/qgis/ServiceAreaFromLayer.py

@@ -37,10 +37,12 @@
                        QgsFeatureSink,
                        QgsFeatureRequest,
                        QgsGeometry,
+                       QgsGeometryUtils,
                        QgsFields,
                        QgsPointXY,
                        QgsField,
                        QgsProcessing,
+                       QgsProcessingParameterBoolean,
                        QgsProcessingParameterEnum,
                        QgsProcessingParameterPoint,
                        QgsProcessingParameterField,
@@ -75,7 +77,9 @@ class ServiceAreaFromLayer(QgisAlgorithm):
     SPEED_FIELD = 'SPEED_FIELD'
     DEFAULT_SPEED = 'DEFAULT_SPEED'
     TOLERANCE = 'TOLERANCE'
+    INCLUDE_BOUNDS = 'INCLUDE_BOUNDS'
     OUTPUT = 'OUTPUT'
+    OUTPUT_LINES = 'OUTPUT_LINES'
 
     def icon(self):
         return QIcon(os.path.join(pluginPath, 'images', 'networkanalysis.svg'))
@@ -146,14 +150,24 @@ def initAlgorithm(self, config=None):
                                                    self.tr('Topology tolerance'),
                                                    QgsProcessingParameterNumber.Double,
                                                    0.0, False, 0, 99999999.99))
-
+        params.append(QgsProcessingParameterBoolean(self.INCLUDE_BOUNDS,
+                                                    self.tr('Include upper/lower bound points'),
+                                                    defaultValue=False))
         for p in params:
             p.setFlags(p.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
             self.addParameter(p)
 
-        self.addParameter(QgsProcessingParameterFeatureSink(self.OUTPUT,
-                                                            self.tr('Service area (boundary nodes)'),
-                                                            QgsProcessing.TypeVectorPoint))
+        lines_output = QgsProcessingParameterFeatureSink(self.OUTPUT_LINES,
+                                                         self.tr('Service area (lines)'),
+                                                         QgsProcessing.TypeVectorLine, optional=True)
+        lines_output.setCreateByDefault(True)
+        self.addParameter(lines_output)
+
+        nodes_output = QgsProcessingParameterFeatureSink(self.OUTPUT,
+                                                         self.tr('Service area (boundary nodes)'),
+                                                         QgsProcessing.TypeVectorPoint, optional=True)
+        nodes_output.setCreateByDefault(False)
+        self.addParameter(nodes_output)
 
     def name(self):
         return 'serviceareafromlayer'
@@ -176,6 +190,10 @@ def processAlgorithm(self, parameters, context, feedback):
         defaultSpeed = self.parameterAsDouble(parameters, self.DEFAULT_SPEED, context)
         tolerance = self.parameterAsDouble(parameters, self.TOLERANCE, context)
 
+        include_bounds = True # default to true to maintain 3.0 API
+        if self.INCLUDE_BOUNDS in parameters:
+            include_bounds = self.parameterAsBool(parameters, self.INCLUDE_BOUNDS, context)
+
         fields = startPoints.fields()
         fields.append(QgsField('type', QVariant.String, '', 254, 0))
         fields.append(QgsField('start', QVariant.String, '', 254, 0))
@@ -240,12 +258,11 @@ def processAlgorithm(self, parameters, context, feedback):
         feedback.pushInfo(QCoreApplication.translate('ServiceAreaFromLayer', 'Calculating service areas…'))
         graph = builder.graph()
 
-        (sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
-                                               fields, QgsWkbTypes.MultiPoint, network.sourceCrs())
+        (point_sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
+                                                     fields, QgsWkbTypes.MultiPoint, network.sourceCrs())
+        (line_sink, line_dest_id) = self.parameterAsSink(parameters, self.OUTPUT_LINES, context,
+                                                         fields, QgsWkbTypes.MultiLineString, network.sourceCrs())
 
-        vertices = []
-        upperBoundary = []
-        lowerBoundary = []
         total = 100.0 / len(snappedPoints) if snappedPoints else 1
         for i, p in enumerate(snappedPoints):
             if feedback.isCanceled():
@@ -255,35 +272,92 @@ def processAlgorithm(self, parameters, context, feedback):
             origPoint = points[i].toString()
 
             tree, cost = QgsGraphAnalyzer.dijkstra(graph, idxStart, 0)
-            for j, v in enumerate(cost):
-                if v > travelCost and tree[j] != -1:
-                    vertexId = graph.edge(tree[j]).fromVertex()
-                    if cost[vertexId] <= travelCost:
-                        vertices.append(j)
-
-            for j in vertices:
-                upperBoundary.append(graph.vertex(graph.edge(tree[j]).toVertex()).point())
-                lowerBoundary.append(graph.vertex(graph.edge(tree[j]).fromVertex()).point())
-
-            geomUpper = QgsGeometry.fromMultiPointXY(upperBoundary)
-            geomLower = QgsGeometry.fromMultiPointXY(lowerBoundary)
-
-            feat.setGeometry(geomUpper)
-
-            attrs = source_attributes[i]
-            attrs.extend(['upper', origPoint])
-            feat.setAttributes(attrs)
-            sink.addFeature(feat, QgsFeatureSink.FastInsert)
-
-            feat.setGeometry(geomLower)
-            attrs[-2] = 'lower'
-            feat.setAttributes(attrs)
-            sink.addFeature(feat, QgsFeatureSink.FastInsert)
 
-            vertices[:] = []
-            upperBoundary[:] = []
-            lowerBoundary[:] = []
+            vertices = set()
+            area_points = []
+            lines = []
+            for vertex, start_vertex_cost in enumerate(cost):
+                inbound_edge_index = tree[vertex]
+                if inbound_edge_index == -1 and vertex != idxStart:
+                    # unreachable vertex
+                    continue
+
+                if start_vertex_cost > travelCost:
+                    # vertex is too expensive, discard
+                    continue
+
+                vertices.add(vertex)
+                start_point = graph.vertex(vertex).point()
+
+                # find all edges coming from this vertex
+                for edge_id in graph.vertex(vertex).outgoingEdges():
+                    edge = graph.edge(edge_id)
+                    end_vertex_cost = start_vertex_cost + edge.cost(0)
+                    end_point = graph.vertex(edge.toVertex()).point()
+                    if end_vertex_cost <= travelCost:
+                        # end vertex is cheap enough to include
+                        vertices.add(edge.toVertex())
+                        lines.append([start_point, end_point])
+                    else:
+                        # travelCost sits somewhere on this edge, interpolate position
+                        interpolated_end_point = QgsGeometryUtils.interpolatePointOnLineByValue(start_point.x(), start_point.y(), start_vertex_cost,
+                                                                                                end_point.x(), end_point.y(), end_vertex_cost, travelCost)
+                        area_points.append(interpolated_end_point)
+                        lines.append([start_point, interpolated_end_point])
+
+            for v in vertices:
+                area_points.append(graph.vertex(v).point())
+
+            feat = QgsFeature()
+            if point_sink is not None:
+                geomPoints = QgsGeometry.fromMultiPointXY(area_points)
+                feat.setGeometry(geomPoints)
+                attrs = source_attributes[i]
+                attrs.extend(['within', origPoint])
+                feat.setAttributes(attrs)
+                point_sink.addFeature(feat, QgsFeatureSink.FastInsert)
+
+                if include_bounds:
+                    upperBoundary = []
+                    lowerBoundary = []
+
+                    vertices = []
+                    for vertex, c in enumerate(cost):
+                        if c > travelCost and tree[vertex] != -1:
+                            vertexId = graph.edge(tree[vertex]).fromVertex()
+                            if cost[vertexId] <= travelCost:
+                                vertices.append(vertex)
+
+                    for v in vertices:
+                        upperBoundary.append(graph.vertex(graph.edge(tree[v]).toVertex()).point())
+                        lowerBoundary.append(graph.vertex(graph.edge(tree[v]).fromVertex()).point())
+
+                    geomUpper = QgsGeometry.fromMultiPointXY(upperBoundary)
+                    geomLower = QgsGeometry.fromMultiPointXY(lowerBoundary)
+
+                    feat.setGeometry(geomUpper)
+                    attrs[-2] = 'upper'
+                    feat.setAttributes(attrs)
+                    point_sink.addFeature(feat, QgsFeatureSink.FastInsert)
+
+                    feat.setGeometry(geomLower)
+                    attrs[-2] = 'lower'
+                    feat.setAttributes(attrs)
+                    point_sink.addFeature(feat, QgsFeatureSink.FastInsert)
+
+            if line_sink is not None:
+                geom_lines = QgsGeometry.fromMultiPolylineXY(lines)
+                feat.setGeometry(geom_lines)
+                attrs = source_attributes[i]
+                attrs.extend(['lines', origPoint])
+                feat.setAttributes(attrs)
+                line_sink.addFeature(feat, QgsFeatureSink.FastInsert)
 
             feedback.setProgress(int(i * total))
 
-        return {self.OUTPUT: dest_id}
+        results = {}
+        if point_sink is not None:
+            results[self.OUTPUT] = dest_id
+        if line_sink is not None:
+            results[self.OUTPUT_LINES] = line_dest_id
+        return results