QGIS API Documentation 3.27.0-Master (75dc696944)
qgsvectorlayerdirector.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgslinevectorlayerdirector.cpp
3 --------------------------------------
4 Date : 2010-10-20
5 Copyright : (C) 2010 by Yakushev Sergey
7****************************************************************************
8* *
9* This program is free software; you can redistribute it and/or modify *
10* it under the terms of the GNU General Public License as published by *
11* the Free Software Foundation; either version 2 of the License, or *
12* (at your option) any later version. *
13* *
14***************************************************************************/
15
23
24#include "qgsfeatureiterator.h"
25#include "qgsfeaturesource.h"
27#include "qgspoint.h"
28#include "qgsgeometry.h"
29#include "qgsdistancearea.h"
30#include "qgswkbtypes.h"
31
32#include <QString>
33#include <QtAlgorithms>
34
35#include <spatialindex/SpatialIndex.h>
36
37using namespace SpatialIndex;
38
40{
41 TiePointInfo() = default;
42 TiePointInfo( int additionalPointId, QgsFeatureId featureId, const QgsPointXY &start, const QgsPointXY &end )
44 , mNetworkFeatureId( featureId )
45 , mFirstPoint( start )
46 , mLastPoint( end )
47 {}
48
51 double mLength = std::numeric_limits<double>::max();
55};
56
58 int directionFieldId,
59 const QString &directDirectionValue,
60 const QString &reverseDirectionValue,
61 const QString &bothDirectionValue,
62 const Direction defaultDirection )
63 : mSource( source )
64 , mDirectionFieldId( directionFieldId )
65 , mDirectDirectionValue( directDirectionValue )
66 , mReverseDirectionValue( reverseDirectionValue )
67 , mBothDirectionValue( bothDirectionValue )
68 , mDefaultDirection( defaultDirection )
69{
70}
71
73{
74 return QStringLiteral( "Vector line" );
75}
76
77QgsAttributeList QgsVectorLayerDirector::requiredAttributes() const
78{
79 QSet< int > attrs;
80
81 if ( mDirectionFieldId != -1 )
82 attrs.insert( mDirectionFieldId );
83
84 for ( const QgsNetworkStrategy *strategy : mStrategies )
85 {
86 attrs.unite( strategy->requiredAttributes() );
87 }
88 return qgis::setToList( attrs );
89}
90
91QgsVectorLayerDirector::Direction QgsVectorLayerDirector::directionForFeature( const QgsFeature &feature ) const
92{
93 if ( mDirectionFieldId < 0 )
94 return mDefaultDirection;
95
96 QString str = feature.attribute( mDirectionFieldId ).toString();
97 if ( str == mBothDirectionValue )
98 {
99 return Direction::DirectionBoth;
100 }
101 else if ( str == mDirectDirectionValue )
102 {
103 return Direction::DirectionForward;
104 }
105 else if ( str == mReverseDirectionValue )
106 {
107 return Direction::DirectionBackward;
108 }
109 else
110 {
111 return mDefaultDirection;
112 }
113}
114
116class QgsNetworkVisitor : public SpatialIndex::IVisitor
117{
118 public:
119 explicit QgsNetworkVisitor( QVector< int > &pointIndexes )
120 : mPoints( pointIndexes ) {}
121
122 void visitNode( const INode &n ) override
123 { Q_UNUSED( n ) }
124
125 void visitData( const IData &d ) override
126 {
127 mPoints.append( d.getIdentifier() );
128 }
129
130 void visitData( std::vector<const IData *> &v ) override
131 { Q_UNUSED( v ) }
132
133 private:
134 QVector< int > &mPoints;
135};
136
138
139std::unique_ptr< SpatialIndex::ISpatialIndex > createVertexSpatialIndex( SpatialIndex::IStorageManager &storageManager )
140{
141 // R-Tree parameters
142 double fillFactor = 0.7;
143 unsigned long indexCapacity = 10;
144 unsigned long leafCapacity = 10;
145 unsigned long dimension = 2;
146 RTree::RTreeVariant variant = RTree::RV_RSTAR;
147
148 // create R-tree
149 SpatialIndex::id_type indexId;
150 std::unique_ptr< SpatialIndex::ISpatialIndex > iRTree( RTree::createNewRTree( storageManager, fillFactor, indexCapacity,
151 leafCapacity, dimension, variant, indexId ) );
152 return iRTree;
153}
154
155int findClosestVertex( const QgsPointXY &point, SpatialIndex::ISpatialIndex *index, double tolerance )
156{
157 QVector< int > matching;
158 QgsNetworkVisitor visitor( matching );
159
160 double pt1[2] = { point.x() - tolerance, point.y() - tolerance },
161 pt2[2] = { point.x() + tolerance, point.y() + tolerance };
162 SpatialIndex::Region searchRegion( pt1, pt2, 2 );
163
164 index->intersectsWithQuery( searchRegion, visitor );
165
166 return matching.empty() ? -1 : matching.at( 0 );
167}
168
169void QgsVectorLayerDirector::makeGraph( QgsGraphBuilderInterface *builder, const QVector< QgsPointXY > &additionalPoints,
170 QVector< QgsPointXY > &snappedPoints, QgsFeedback *feedback ) const
171{
172 long featureCount = mSource->featureCount() * 2;
173 int step = 0;
174
176 ct.setSourceCrs( mSource->sourceCrs() );
177 if ( builder->coordinateTransformationEnabled() )
178 {
179 ct.setDestinationCrs( builder->destinationCrs() );
180 }
181
182 // clear existing snapped points list, and resize to length of provided additional points
183 snappedPoints = QVector< QgsPointXY >( additionalPoints.size(), QgsPointXY( 0.0, 0.0 ) );
184 // tie points = snapped location of specified additional points to network lines
185 QVector< TiePointInfo > additionalTiePoints( additionalPoints.size() );
186
187 // graph's vertices = all vertices in graph, with vertices within builder's tolerance collapsed together
188 QVector< QgsPointXY > graphVertices;
189
190 // spatial index for graph vertices
191 std::unique_ptr< SpatialIndex::IStorageManager > iStorage( StorageManager::createNewMemoryStorageManager() );
192 std::unique_ptr< SpatialIndex::ISpatialIndex > iRTree = createVertexSpatialIndex( *iStorage );
193
194 double tolerance = std::max( builder->topologyTolerance(), 1e-10 );
195 auto findPointWithinTolerance = [&iRTree, tolerance]( const QgsPointXY & point )->int
196 {
197 return findClosestVertex( point, iRTree.get(), tolerance );
198 };
199 auto addPointToIndex = [&iRTree]( const QgsPointXY & point, int index )
200 {
201 double coords[] = {point.x(), point.y()};
202 iRTree->insertData( 0, nullptr, SpatialIndex::Point( coords, 2 ), index );
203 };
204
205 // first iteration - get all nodes from network, and snap additional points to network
207 QgsFeature feature;
208 while ( fit.nextFeature( feature ) )
209 {
210 if ( feedback && feedback->isCanceled() )
211 return;
212
215 mpl = feature.geometry().asMultiPolyline();
217 mpl.push_back( feature.geometry().asPolyline() );
218
219 for ( const QgsPolylineXY &line : std::as_const( mpl ) )
220 {
221 QgsPointXY pt1, pt2;
222 bool isFirstPoint = true;
223 for ( const QgsPointXY &point : line )
224 {
225 pt2 = ct.transform( point );
226
227 int pt2Idx = findPointWithinTolerance( pt2 ) ;
228 if ( pt2Idx == -1 )
229 {
230 // no vertex already exists within tolerance - add to points, and index
231 addPointToIndex( pt2, graphVertices.count() );
232 graphVertices.push_back( pt2 );
233 }
234 else
235 {
236 // vertex already exists within tolerance - use that
237 pt2 = graphVertices.at( pt2Idx );
238 }
239
240 if ( !isFirstPoint )
241 {
242 // check if this line segment is a candidate for being closest to each additional point
243 int i = 0;
244 for ( const QgsPointXY &additionalPoint : additionalPoints )
245 {
246
247 QgsPointXY snappedPoint;
248 double thisSegmentClosestDist = std::numeric_limits<double>::max();
249 if ( pt1 == pt2 )
250 {
251 thisSegmentClosestDist = additionalPoint.sqrDist( pt1 );
252 snappedPoint = pt1;
253 }
254 else
255 {
256 thisSegmentClosestDist = additionalPoint.sqrDistToSegment( pt1.x(), pt1.y(),
257 pt2.x(), pt2.y(), snappedPoint, 0 );
258 }
259
260 if ( thisSegmentClosestDist < additionalTiePoints[ i ].mLength )
261 {
262 // found a closer segment for this additional point
263 TiePointInfo info( i, feature.id(), pt1, pt2 );
264 info.mLength = thisSegmentClosestDist;
265 info.mTiedPoint = snappedPoint;
266
267 additionalTiePoints[ i ] = info;
268 snappedPoints[ i ] = info.mTiedPoint;
269 }
270 i++;
271 }
272 }
273 pt1 = pt2;
274 isFirstPoint = false;
275 }
276 }
277 if ( feedback )
278 feedback->setProgress( 100.0 * static_cast< double >( ++step ) / featureCount );
279 }
280
281 // build a hash of feature ids to tie points which depend on this feature
282 QHash< QgsFeatureId, QList< int > > tiePointNetworkFeatures;
283 int i = 0;
284 for ( TiePointInfo &info : additionalTiePoints )
285 {
286 tiePointNetworkFeatures[ info.mNetworkFeatureId ] << i;
287 i++;
288 }
289
290 // add tied point to graph
291 for ( int i = 0; i < snappedPoints.size(); ++i )
292 {
293 // check index to see if vertex exists within tolerance of tie point
294 const QgsPointXY point = snappedPoints.at( i );
295 int ptIdx = findPointWithinTolerance( point );
296 if ( ptIdx == -1 )
297 {
298 // no vertex already within tolerance, add to index and network vertices
299 addPointToIndex( point, graphVertices.count() );
300 graphVertices.push_back( point );
301 }
302 else
303 {
304 // otherwise snap tie point to vertex
305 snappedPoints[ i ] = graphVertices.at( ptIdx );
306 }
307 }
308 // also need to update tie points - they need to be matched for snapped points
309 for ( int i = 0; i < additionalTiePoints.count(); ++i )
310 {
311 additionalTiePoints[ i ].mTiedPoint = snappedPoints.at( additionalTiePoints.at( i ).additionalPointId );
312 }
313
314
315 // begin graph construction
316
317 // add vertices to graph
318 {
319 int i = 0;
320 for ( const QgsPointXY &point : graphVertices )
321 {
322 builder->addVertex( i, point );
323 i++;
324 }
325 }
326
327 fit = mSource->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( requiredAttributes() ) );
328 while ( fit.nextFeature( feature ) )
329 {
330 if ( feedback && feedback->isCanceled() )
331 return;
332
333 Direction direction = directionForFeature( feature );
334
335 // begin features segments and add arc to the Graph;
338 mpl = feature.geometry().asMultiPolyline();
340 mpl.push_back( feature.geometry().asPolyline() );
341
342 for ( const QgsPolylineXY &line : std::as_const( mpl ) )
343 {
344 QgsPointXY pt1, pt2;
345
346 bool isFirstPoint = true;
347 for ( const QgsPointXY &point : line )
348 {
349 pt2 = ct.transform( point );
350 int pPt2idx = findPointWithinTolerance( pt2 );
351 Q_ASSERT_X( pPt2idx >= 0, "QgsVectorLayerDirectory::makeGraph", "encountered a vertex which was not present in graph" );
352 pt2 = graphVertices.at( pPt2idx );
353
354 if ( !isFirstPoint )
355 {
356 QMap< double, QgsPointXY > pointsOnArc;
357 pointsOnArc[ 0.0 ] = pt1;
358 pointsOnArc[ pt1.sqrDist( pt2 )] = pt2;
359
360 const QList< int > tiePointsForCurrentFeature = tiePointNetworkFeatures.value( feature.id() );
361 for ( int tiePointIdx : tiePointsForCurrentFeature )
362 {
363 const TiePointInfo &t = additionalTiePoints.at( tiePointIdx );
364 if ( t.mFirstPoint == pt1 && t.mLastPoint == pt2 )
365 {
366 pointsOnArc[ pt1.sqrDist( t.mTiedPoint )] = t.mTiedPoint;
367 }
368 }
369
370 QgsPointXY arcPt1;
371 QgsPointXY arcPt2;
372 int pt1idx = -1;
373 int pt2idx = -1;
374 bool isFirstPoint = true;
375 for ( auto arcPointIt = pointsOnArc.constBegin(); arcPointIt != pointsOnArc.constEnd(); ++arcPointIt )
376 {
377 arcPt2 = arcPointIt.value();
378
379 pt2idx = findPointWithinTolerance( arcPt2 );
380 Q_ASSERT_X( pt2idx >= 0, "QgsVectorLayerDirectory::makeGraph", "encountered a vertex which was not present in graph" );
381 arcPt2 = graphVertices.at( pt2idx );
382
383 if ( !isFirstPoint && arcPt1 != arcPt2 )
384 {
385 double distance = builder->distanceArea()->measureLine( arcPt1, arcPt2 );
386 QVector< QVariant > prop;
387 prop.reserve( mStrategies.size() );
388 for ( QgsNetworkStrategy *strategy : mStrategies )
389 {
390 prop.push_back( strategy->cost( distance, feature ) );
391 }
392
393 if ( direction == Direction::DirectionForward ||
394 direction == Direction::DirectionBoth )
395 {
396 builder->addEdge( pt1idx, arcPt1, pt2idx, arcPt2, prop );
397 }
398 if ( direction == Direction::DirectionBackward ||
399 direction == Direction::DirectionBoth )
400 {
401 builder->addEdge( pt2idx, arcPt2, pt1idx, arcPt1, prop );
402 }
403 }
404 pt1idx = pt2idx;
405 arcPt1 = arcPt2;
406 isFirstPoint = false;
407 }
408 }
409 pt1 = pt2;
410 isFirstPoint = false;
411 }
412 }
413 if ( feedback )
414 {
415 feedback->setProgress( 100.0 * static_cast< double >( ++step ) / featureCount );
416 }
417
418 }
419}
Class for doing transforms between two map coordinate systems.
void setSourceCrs(const QgsCoordinateReferenceSystem &crs)
Sets the source coordinate reference system.
void setDestinationCrs(const QgsCoordinateReferenceSystem &crs)
Sets the destination coordinate reference system.
QgsPointXY transform(const QgsPointXY &point, Qgis::TransformDirection direction=Qgis::TransformDirection::Forward) const SIP_THROW(QgsCsException)
Transform the point from the source CRS to the destination CRS.
double measureLine(const QVector< QgsPointXY > &points) const
Measures the length of a line with multiple segments.
Wrapper for iterator of features from vector data provider or vector layer.
bool nextFeature(QgsFeature &f)
This class wraps a request for features to a vector layer (or directly its vector data provider).
QgsFeatureRequest & setNoAttributes()
Set that no attributes will be fetched.
An interface for objects which provide features via a getFeatures method.
virtual QgsCoordinateReferenceSystem sourceCrs() const =0
Returns the coordinate reference system for features in the source.
virtual QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const =0
Returns an iterator for the features in the source.
virtual long long featureCount() const =0
Returns the number of features contained in the source, or -1 if the feature count is unknown.
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition: qgsfeature.h:56
QgsGeometry geometry
Definition: qgsfeature.h:67
QVariant attribute(const QString &name) const
Lookup attribute value by attribute name.
Definition: qgsfeature.cpp:335
Q_GADGET QgsFeatureId id
Definition: qgsfeature.h:64
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition: qgsfeedback.h:45
bool isCanceled() const SIP_HOLDGIL
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
QgsWkbTypes::Type wkbType() const SIP_HOLDGIL
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
QgsPolylineXY asPolyline() const
Returns the contents of the geometry as a polyline.
QgsMultiPolylineXY asMultiPolyline() const
Returns the contents of the geometry as a multi-linestring.
Determine interface for creating a graph.
virtual void addVertex(int id, const QgsPointXY &pt)
Add vertex to the graph.
QgsCoordinateReferenceSystem destinationCrs() const
Returns destinaltion CRS.
QgsDistanceArea * distanceArea()
Returns measurement tool.
virtual void addEdge(int pt1id, const QgsPointXY &pt1, int pt2id, const QgsPointXY &pt2, const QVector< QVariant > &strategies)
Add edge to the graph.
double topologyTolerance() const
Returns topology tolerance.
bool coordinateTransformationEnabled() const
Returns coordinate transformation enabled.
QList< QgsNetworkStrategy * > mStrategies
QgsNetworkStrategy defines strategy used for calculation of the edge cost.
A class to represent a 2D point.
Definition: qgspointxy.h:59
double sqrDist(double x, double y) const SIP_HOLDGIL
Returns the squared distance between this point a specified x, y coordinate.
Definition: qgspointxy.h:190
double y
Definition: qgspointxy.h:63
Q_GADGET double x
Definition: qgspointxy.h:62
double sqrDistToSegment(double x1, double y1, double x2, double y2, QgsPointXY &minDistPoint, double epsilon=DEFAULT_SEGMENT_EPSILON) const SIP_HOLDGIL
Returns the minimum distance between this point and a segment.
Definition: qgspointxy.cpp:95
QString name() const override
Returns director name.
void makeGraph(QgsGraphBuilderInterface *builder, const QVector< QgsPointXY > &additionalPoints, QVector< QgsPointXY > &snappedPoints, QgsFeedback *feedback=nullptr) const override
Make a graph using QgsGraphBuilder.
Direction
Edge direction Edge can be one-way with direct flow (one can move only from the start point to the en...
QgsVectorLayerDirector(QgsFeatureSource *source, int directionFieldId, const QString &directDirectionValue, const QString &reverseDirectionValue, const QString &bothDirectionValue, Direction defaultDirection)
Default constructor.
static Type flatType(Type type) SIP_HOLDGIL
Returns the flat type for a WKB type.
Definition: qgswkbtypes.h:732
#define str(x)
Definition: qgis.cpp:37
qint64 QgsFeatureId
64 bit feature ids negative numbers are used for uncommitted/newly added features
Definition: qgsfeatureid.h:28
QList< int > QgsAttributeList
Definition: qgsfield.h:26
QVector< QgsPolylineXY > QgsMultiPolylineXY
A collection of QgsPolylines that share a common collection of attributes.
Definition: qgsgeometry.h:86
QVector< QgsPointXY > QgsPolylineXY
Polyline as represented as a vector of two-dimensional points.
Definition: qgsgeometry.h:63
std::unique_ptr< SpatialIndex::ISpatialIndex > createVertexSpatialIndex(SpatialIndex::IStorageManager &storageManager)
int findClosestVertex(const QgsPointXY &point, SpatialIndex::ISpatialIndex *index, double tolerance)
TiePointInfo(int additionalPointId, QgsFeatureId featureId, const QgsPointXY &start, const QgsPointXY &end)
TiePointInfo()=default
QgsFeatureId mNetworkFeatureId