Skip to content

Commit

Permalink
[processing] Fixes for Service Area algorithms
Browse files Browse the repository at this point in the history 
- 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.

(cherry-picked from ccb72eb)
  • Loading branch information
@nyalldawson
nyalldawson committed Apr 6, 2018
1 parent d93862c commit 017f6a5
Show file tree
Hide file tree
Showing 23 changed files with 480 additions and 76 deletions.
3 changes: 2 additions & 1 deletion python/analysis/network/qgsgraphanalyzer.sip.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
148 changes: 111 additions & 37 deletions python/plugins/processing/algs/qgis/ServiceAreaFromLayer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,10 +37,12 @@
QgsFeatureSink,
QgsFeatureRequest,
QgsGeometry,
QgsGeometryUtils,
QgsFields,
QgsPointXY,
QgsField,
QgsProcessing,
QgsProcessingParameterBoolean,
QgsProcessingParameterEnum,
QgsProcessingParameterPoint,
QgsProcessingParameterField,
Expand Down Expand Up @@ -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'))
Expand Down Expand Up @@ -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'
Expand All @@ -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))
Expand Down Expand Up @@ -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():
Expand All @@ -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
Loading

0 comments on commit 017f6a5

Please sign in to comment.