QGIS API Documentation 3.43.0-Master (e01d6d7c4c0)
qgsalgorithmserviceareafromlayer.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsalgorithmserviceareafromlayer.cpp
3 ---------------------
4 begin : July 2018
5 copyright : (C) 2018 by Alexander Bruy
6 email : alexander dot bruy at gmail dot com
7 ***************************************************************************/
8
9/***************************************************************************
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 ***************************************************************************/
17
19
20#include "qgsgeometryutils.h"
21#include "qgsgraphanalyzer.h"
22
24
25QString QgsServiceAreaFromLayerAlgorithm::name() const
26{
27 return QStringLiteral( "serviceareafromlayer" );
28}
29
30QString QgsServiceAreaFromLayerAlgorithm::displayName() const
31{
32 return QObject::tr( "Service area (from layer)" );
33}
34
35QStringList QgsServiceAreaFromLayerAlgorithm::tags() const
36{
37 return QObject::tr( "network,service,area,shortest,fastest" ).split( ',' );
38}
39
40QString QgsServiceAreaFromLayerAlgorithm::shortHelpString() const
41{
42 return QObject::tr( "This algorithm creates a new vector layer with all the edges or parts of "
43 "edges of a network line layer that can be reached within a distance "
44 "or a time, starting from features of a point layer. The distance and "
45 "the time (both referred to as \"travel cost\") must be specified "
46 "respectively in the network layer units or in hours." );
47}
48
49QString QgsServiceAreaFromLayerAlgorithm::shortDescription() const
50{
51 return QObject::tr( "Creates a vector layer with all the edges or parts of "
52 "edges of a network line layer that can be reached within a distance "
53 "or a time, starting from features of a point layer." );
54}
55
56QgsServiceAreaFromLayerAlgorithm *QgsServiceAreaFromLayerAlgorithm::createInstance() const
57{
58 return new QgsServiceAreaFromLayerAlgorithm();
59}
60
61void QgsServiceAreaFromLayerAlgorithm::initAlgorithm( const QVariantMap & )
62{
63 addCommonParams();
64 addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "START_POINTS" ), QObject::tr( "Vector layer with start points" ), QList<int>() << static_cast<int>( Qgis::ProcessingSourceType::VectorPoint ) ) );
65
66 auto travelCost = std::make_unique<QgsProcessingParameterNumber>( QStringLiteral( "TRAVEL_COST" ), QObject::tr( "Travel cost (distance for 'Shortest', time for 'Fastest')" ), Qgis::ProcessingNumberParameterType::Double, 0, true, 0 );
67 travelCost->setFlags( travelCost->flags() | Qgis::ProcessingParameterFlag::Hidden );
68 addParameter( travelCost.release() );
69
70 addParameter( new QgsProcessingParameterNumber( QStringLiteral( "TRAVEL_COST2" ), QObject::tr( "Travel cost (distance for 'Shortest', time for 'Fastest')" ), Qgis::ProcessingNumberParameterType::Double, 0, false, 0 ) );
71
72 auto includeBounds = std::make_unique<QgsProcessingParameterBoolean>( QStringLiteral( "INCLUDE_BOUNDS" ), QObject::tr( "Include upper/lower bound points" ), false, true );
73 includeBounds->setFlags( includeBounds->flags() | Qgis::ProcessingParameterFlag::Advanced );
74 addParameter( includeBounds.release() );
75
76 std::unique_ptr<QgsProcessingParameterNumber> maxPointDistanceFromNetwork = std::make_unique<QgsProcessingParameterDistance>( QStringLiteral( "POINT_TOLERANCE" ), QObject::tr( "Maximum point distance from network" ), QVariant(), QStringLiteral( "INPUT" ), true, 0 );
77 maxPointDistanceFromNetwork->setFlags( maxPointDistanceFromNetwork->flags() | Qgis::ProcessingParameterFlag::Advanced );
78 maxPointDistanceFromNetwork->setHelp( QObject::tr( "Specifies an optional limit on the distance from the points to the network layer. If a point is further from the network than this distance it will be treated as non-routable." ) );
79 addParameter( maxPointDistanceFromNetwork.release() );
80
81 auto outputLines = std::make_unique<QgsProcessingParameterFeatureSink>( QStringLiteral( "OUTPUT_LINES" ), QObject::tr( "Service area (lines)" ), Qgis::ProcessingSourceType::VectorLine, QVariant(), true );
82 outputLines->setCreateByDefault( true );
83 addParameter( outputLines.release() );
84
85 auto outputPoints = std::make_unique<QgsProcessingParameterFeatureSink>( QStringLiteral( "OUTPUT" ), QObject::tr( "Service area (boundary nodes)" ), Qgis::ProcessingSourceType::VectorPoint, QVariant(), true );
86 outputPoints->setCreateByDefault( false );
87 addParameter( outputPoints.release() );
88
89 auto outputNonRoutable = std::make_unique<QgsProcessingParameterFeatureSink>( QStringLiteral( "OUTPUT_NON_ROUTABLE" ), QObject::tr( "Non-routable features" ), Qgis::ProcessingSourceType::VectorPoint, QVariant(), true );
90 outputNonRoutable->setHelp( QObject::tr( "An optional output which will be used to store any input features which could not be routed (e.g. those which are too far from the network layer)." ) );
91 outputNonRoutable->setCreateByDefault( false );
92 addParameter( outputNonRoutable.release() );
93}
94
95QVariantMap QgsServiceAreaFromLayerAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
96{
97 loadCommonParams( parameters, context, feedback );
98
99 std::unique_ptr<QgsFeatureSource> startPoints( parameterAsSource( parameters, QStringLiteral( "START_POINTS" ), context ) );
100 if ( !startPoints )
101 throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "START_POINTS" ) ) );
102
103 // use older deprecated travel cost style if specified, to maintain old api
104 const bool useOldTravelCost = parameters.value( QStringLiteral( "TRAVEL_COST" ) ).isValid();
105 double travelCost = parameterAsDouble( parameters, useOldTravelCost ? QStringLiteral( "TRAVEL_COST" ) : QStringLiteral( "TRAVEL_COST2" ), context );
106
107 int strategy = parameterAsInt( parameters, QStringLiteral( "STRATEGY" ), context );
108 if ( strategy && !useOldTravelCost )
109 travelCost *= mMultiplier;
110
111 bool includeBounds = true; // default to true to maintain 3.0 API
112 if ( parameters.contains( QStringLiteral( "INCLUDE_BOUNDS" ) ) )
113 {
114 includeBounds = parameterAsBool( parameters, QStringLiteral( "INCLUDE_BOUNDS" ), context );
115 }
116
117 QVector<QgsPointXY> points;
118 QHash<int, QgsAttributes> sourceAttributes;
119 loadPoints( startPoints.get(), points, sourceAttributes, context, feedback );
120
121 feedback->pushInfo( QObject::tr( "Building graph…" ) );
122 QVector<QgsPointXY> snappedPoints;
123 mDirector->makeGraph( mBuilder.get(), points, snappedPoints, feedback );
124
125 feedback->pushInfo( QObject::tr( "Calculating service areas…" ) );
126 std::unique_ptr<QgsGraph> graph( mBuilder->takeGraph() );
127
128 QgsFields newFields;
129 newFields.append( QgsField( QStringLiteral( "type" ), QMetaType::Type::QString ) );
130 newFields.append( QgsField( QStringLiteral( "start" ), QMetaType::Type::QString ) );
131 QgsFields fields = QgsProcessingUtils::combineFields( startPoints->fields(), newFields );
132
133 QString pointsSinkId;
134 std::unique_ptr<QgsFeatureSink> pointsSink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, pointsSinkId, fields, Qgis::WkbType::MultiPoint, mNetwork->sourceCrs() ) );
135
136 QString linesSinkId;
137 std::unique_ptr<QgsFeatureSink> linesSink( parameterAsSink( parameters, QStringLiteral( "OUTPUT_LINES" ), context, linesSinkId, fields, Qgis::WkbType::MultiLineString, mNetwork->sourceCrs() ) );
138
139 QString nonRoutableSinkId;
140 std::unique_ptr<QgsFeatureSink> nonRoutableSink( parameterAsSink( parameters, QStringLiteral( "OUTPUT_NON_ROUTABLE" ), context, nonRoutableSinkId, startPoints->fields(), Qgis::WkbType::Point, mNetwork->sourceCrs() ) );
141
142 const double pointDistanceThreshold = parameters.value( QStringLiteral( "POINT_TOLERANCE" ) ).isValid() ? parameterAsDouble( parameters, QStringLiteral( "POINT_TOLERANCE" ), context ) : -1;
143
144 int idxStart;
145 QVector<int> tree;
146 QVector<double> costs;
147
148 int inboundEdgeIndex;
149 double startVertexCost, endVertexCost;
150 QgsPointXY startPoint, endPoint;
151 QgsGraphEdge edge;
152
153 QgsFeature feat;
154 QgsAttributes attributes;
155
156 const double step = snappedPoints.size() > 0 ? 100.0 / snappedPoints.size() : 1;
157 for ( int i = 0; i < snappedPoints.size(); i++ )
158 {
159 if ( feedback->isCanceled() )
160 {
161 break;
162 }
163
164 const QgsPointXY snappedPoint = snappedPoints.at( i );
165 const QgsPointXY originalPoint = points.at( i );
166
167 if ( pointDistanceThreshold >= 0 )
168 {
169 double distancePointToNetwork = 0;
170 try
171 {
172 distancePointToNetwork = mBuilder->distanceArea()->measureLine( originalPoint, snappedPoint );
173 }
174 catch ( QgsCsException & )
175 {
176 throw QgsProcessingException( QObject::tr( "An error occurred while calculating length" ) );
177 }
178
179 if ( distancePointToNetwork > pointDistanceThreshold )
180 {
181 feedback->pushWarning( QObject::tr( "Point is too far from the network layer (%1, maximum permitted is %2)" ).arg( distancePointToNetwork ).arg( pointDistanceThreshold ) );
182 if ( nonRoutableSink )
183 {
184 feat.setGeometry( QgsGeometry::fromPointXY( originalPoint ) );
185 attributes = sourceAttributes.value( i + 1 );
186 feat.setAttributes( attributes );
187 if ( !nonRoutableSink->addFeature( feat, QgsFeatureSink::FastInsert ) )
188 throw QgsProcessingException( writeFeatureError( nonRoutableSink.get(), parameters, QStringLiteral( "OUTPUT_NON_ROUTABLE" ) ) );
189 }
190
191 feedback->setProgress( i * step );
192 continue;
193 }
194 }
195
196 const QString originalPointString = originalPoint.toString();
197
198 idxStart = graph->findVertex( snappedPoint );
199
200 QgsGraphAnalyzer::dijkstra( graph.get(), idxStart, 0, &tree, &costs );
201
202 QgsMultiPointXY areaPoints;
203 QgsMultiPolylineXY lines;
204 QSet<int> vertices;
205
206 for ( int j = 0; j < costs.size(); j++ )
207 {
208 inboundEdgeIndex = tree.at( j );
209
210 if ( inboundEdgeIndex == -1 && j != idxStart )
211 {
212 // unreachable vertex
213 continue;
214 }
215
216 startVertexCost = costs.at( j );
217 if ( startVertexCost > travelCost )
218 {
219 // vertex is too expensive, discard
220 continue;
221 }
222
223 vertices.insert( j );
224 startPoint = graph->vertex( j ).point();
225
226 // find all edges coming from this vertex
227 const QList<int> outgoingEdges = graph->vertex( j ).outgoingEdges();
228 for ( int edgeId : outgoingEdges )
229 {
230 edge = graph->edge( edgeId );
231 endVertexCost = startVertexCost + edge.cost( 0 ).toDouble();
232 endPoint = graph->vertex( edge.toVertex() ).point();
233 if ( endVertexCost <= travelCost )
234 {
235 // end vertex is cheap enough to include
236 vertices.insert( edge.toVertex() );
237 lines.push_back( QgsPolylineXY() << startPoint << endPoint );
238 }
239 else
240 {
241 // travelCost sits somewhere on this edge, interpolate position
242 QgsPointXY interpolatedEndPoint = QgsGeometryUtils::interpolatePointOnLineByValue( startPoint.x(), startPoint.y(), startVertexCost, endPoint.x(), endPoint.y(), endVertexCost, travelCost );
243
244 areaPoints.push_back( interpolatedEndPoint );
245 lines.push_back( QgsPolylineXY() << startPoint << interpolatedEndPoint );
246 }
247 } // edges
248 } // costs
249
250 // convert to list and sort to maintain same order of points between algorithm runs
251 QList<int> verticesList = qgis::setToList( vertices );
252 areaPoints.reserve( verticesList.size() );
253 std::sort( verticesList.begin(), verticesList.end() );
254 for ( int v : verticesList )
255 {
256 areaPoints.push_back( graph->vertex( v ).point() );
257 }
258
259 if ( pointsSink )
260 {
261 QgsGeometry geomPoints = QgsGeometry::fromMultiPointXY( areaPoints );
262 feat.setGeometry( geomPoints );
263 attributes = sourceAttributes.value( i + 1 );
264 attributes << QStringLiteral( "within" ) << originalPointString;
265 feat.setAttributes( attributes );
266 if ( !pointsSink->addFeature( feat, QgsFeatureSink::FastInsert ) )
267 throw QgsProcessingException( writeFeatureError( pointsSink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
268
269 if ( includeBounds )
270 {
271 QgsMultiPointXY upperBoundary, lowerBoundary;
272 QVector<int> nodes;
273 nodes.reserve( costs.size() );
274
275 int vertexId;
276 for ( int v = 0; v < costs.size(); v++ )
277 {
278 if ( costs.at( v ) > travelCost && tree.at( v ) != -1 )
279 {
280 vertexId = graph->edge( tree.at( v ) ).fromVertex();
281 if ( costs.at( vertexId ) <= travelCost )
282 {
283 nodes.push_back( v );
284 }
285 }
286 } // costs
287
288 upperBoundary.reserve( nodes.size() );
289 lowerBoundary.reserve( nodes.size() );
290 for ( int n : std::as_const( nodes ) )
291 {
292 upperBoundary.push_back( graph->vertex( graph->edge( tree.at( n ) ).toVertex() ).point() );
293 lowerBoundary.push_back( graph->vertex( graph->edge( tree.at( n ) ).fromVertex() ).point() );
294 } // nodes
295
296 QgsGeometry geomUpper = QgsGeometry::fromMultiPointXY( upperBoundary );
297 QgsGeometry geomLower = QgsGeometry::fromMultiPointXY( lowerBoundary );
298
299 feat.setGeometry( geomUpper );
300 attributes = sourceAttributes.value( i + 1 );
301 attributes << QStringLiteral( "upper" ) << originalPointString;
302 feat.setAttributes( attributes );
303 if ( !pointsSink->addFeature( feat, QgsFeatureSink::FastInsert ) )
304 throw QgsProcessingException( writeFeatureError( pointsSink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
305
306 feat.setGeometry( geomLower );
307 attributes = sourceAttributes.value( i + 1 );
308 attributes << QStringLiteral( "lower" ) << originalPointString;
309 feat.setAttributes( attributes );
310 if ( !pointsSink->addFeature( feat, QgsFeatureSink::FastInsert ) )
311 throw QgsProcessingException( writeFeatureError( pointsSink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
312 } // includeBounds
313 }
314
315 if ( linesSink )
316 {
318 feat.setGeometry( geomLines );
319 attributes = sourceAttributes.value( i + 1 );
320 attributes << QStringLiteral( "lines" ) << originalPointString;
321 feat.setAttributes( attributes );
322 if ( !linesSink->addFeature( feat, QgsFeatureSink::FastInsert ) )
323 throw QgsProcessingException( writeFeatureError( linesSink.get(), parameters, QStringLiteral( "OUTPUT_LINES" ) ) );
324 }
325
326 feedback->setProgress( i * step );
327 } // snappedPoints
328
329 QVariantMap outputs;
330 if ( pointsSink )
331 {
332 pointsSink->finalize();
333 outputs.insert( QStringLiteral( "OUTPUT" ), pointsSinkId );
334 }
335 if ( linesSink )
336 {
337 linesSink->finalize();
338 outputs.insert( QStringLiteral( "OUTPUT_LINES" ), linesSinkId );
339 }
340 if ( nonRoutableSink )
341 {
342 nonRoutableSink->finalize();
343 outputs.insert( QStringLiteral( "OUTPUT_NON_ROUTABLE" ), nonRoutableSinkId );
344 }
345
346 return outputs;
347}
348
@ VectorPoint
Vector point layers.
@ VectorLine
Vector line layers.
@ MultiPoint
MultiPoint.
@ MultiLineString
MultiLineString.
@ Hidden
Parameter is hidden and should not be shown to users.
@ Advanced
Parameter is an advanced parameter which should be hidden from users by default.
A vector of attributes.
Custom exception class for Coordinate Reference System related exceptions.
@ FastInsert
Use faster inserts, at the cost of updating the passed features to reflect changes made at the provid...
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition qgsfeature.h:58
void setAttributes(const QgsAttributes &attrs)
Sets the feature's attributes.
void setGeometry(const QgsGeometry &geometry)
Set the feature's geometry.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition qgsfeedback.h:53
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition qgsfeedback.h:61
Encapsulate a field in an attribute table or data source.
Definition qgsfield.h:53
Container of fields for a vector layer.
Definition qgsfields.h:46
bool append(const QgsField &field, Qgis::FieldOrigin origin=Qgis::FieldOrigin::Provider, int originIndex=-1)
Appends a field.
Definition qgsfields.cpp:70
static QgsPointXY interpolatePointOnLineByValue(double x1, double y1, double v1, double x2, double y2, double v2, double value)
Interpolates the position of a point along the line from (x1, y1) to (x2, y2).
A geometry is the spatial representation of a feature.
static QgsGeometry fromMultiPolylineXY(const QgsMultiPolylineXY &multiline)
Creates a new geometry from a QgsMultiPolylineXY object.
static QgsGeometry fromMultiPointXY(const QgsMultiPointXY &multipoint)
Creates a new geometry from a QgsMultiPointXY object.
static QgsGeometry fromPointXY(const QgsPointXY &point)
Creates a new geometry from a QgsPointXY object.
static void dijkstra(const QgsGraph *source, int startVertexIdx, int criterionNum, QVector< int > *resultTree=nullptr, QVector< double > *resultCost=nullptr)
Solve shortest path problem using Dijkstra algorithm.
Represents an edge in a graph.
Definition qgsgraph.h:44
int toVertex() const
Returns the index of the vertex at the end of this edge.
Definition qgsgraph.cpp:184
QVariant cost(int strategyIndex) const
Returns edge cost calculated using specified strategy.
Definition qgsgraph.cpp:169
Represents a 2D point.
Definition qgspointxy.h:60
QString toString(int precision=-1) const
Returns a string representation of the point (x, y) with a preset precision.
double y
Definition qgspointxy.h:64
double x
Definition qgspointxy.h:63
Contains information about the context in which a processing algorithm is executed.
Custom exception class for processing related exceptions.
Base class for providing feedback from a processing algorithm.
virtual void pushInfo(const QString &info)
Pushes a general informational message from the algorithm.
virtual void pushWarning(const QString &warning)
Pushes a warning informational message from the algorithm.
An input feature source (such as vector layers) parameter for processing algorithms.
A numeric parameter for processing algorithms.
static QgsFields combineFields(const QgsFields &fieldsA, const QgsFields &fieldsB, const QString &fieldsBPrefix=QString())
Combines two field lists, avoiding duplicate field names (in a case-insensitive manner).
QVector< QgsPolylineXY > QgsMultiPolylineXY
A collection of QgsPolylines that share a common collection of attributes.
Definition qgsgeometry.h:84
QVector< QgsPointXY > QgsMultiPointXY
A collection of QgsPoints that share a common collection of attributes.
Definition qgsgeometry.h:80
QVector< QgsPointXY > QgsPolylineXY
Polyline as represented as a vector of two-dimensional points.
Definition qgsgeometry.h:62