QGIS API Documentation 3.30.0-'s-Hertogenbosch (f186b8efe0)
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 "qgsgeometry.h"
28#include "qgsdistancearea.h"
29#include "qgswkbtypes.h"
30
31#include <QString>
32#include <QtAlgorithms>
33
34#include <spatialindex/SpatialIndex.h>
35
36using namespace SpatialIndex;
37
39{
40 TiePointInfo() = default;
41 TiePointInfo( int additionalPointId, QgsFeatureId featureId, const QgsPointXY &start, const QgsPointXY &end )
43 , mNetworkFeatureId( featureId )
44 , mFirstPoint( start )
45 , mLastPoint( end )
46 {}
47
50 double mLength = std::numeric_limits<double>::max();
54};
55
57 int directionFieldId,
58 const QString &directDirectionValue,
59 const QString &reverseDirectionValue,
60 const QString &bothDirectionValue,
61 const Direction defaultDirection )
62 : mSource( source )
63 , mDirectionFieldId( directionFieldId )
64 , mDirectDirectionValue( directDirectionValue )
65 , mReverseDirectionValue( reverseDirectionValue )
66 , mBothDirectionValue( bothDirectionValue )
67 , mDefaultDirection( defaultDirection )
68{
69}
70
72{
73 return QStringLiteral( "Vector line" );
74}
75
76QgsAttributeList QgsVectorLayerDirector::requiredAttributes() const
77{
78 QSet< int > attrs;
79
80 if ( mDirectionFieldId != -1 )
81 attrs.insert( mDirectionFieldId );
82
83 for ( const QgsNetworkStrategy *strategy : mStrategies )
84 {
85 attrs.unite( strategy->requiredAttributes() );
86 }
87 return qgis::setToList( attrs );
88}
89
90QgsVectorLayerDirector::Direction QgsVectorLayerDirector::directionForFeature( const QgsFeature &feature ) const
91{
92 if ( mDirectionFieldId < 0 )
93 return mDefaultDirection;
94
95 QString str = feature.attribute( mDirectionFieldId ).toString();
96 if ( str == mBothDirectionValue )
97 {
98 return Direction::DirectionBoth;
99 }
100 else if ( str == mDirectDirectionValue )
101 {
102 return Direction::DirectionForward;
103 }
104 else if ( str == mReverseDirectionValue )
105 {
106 return Direction::DirectionBackward;
107 }
108 else
109 {
110 return mDefaultDirection;
111 }
112}
113
115class QgsNetworkVisitor : public SpatialIndex::IVisitor
116{
117 public:
118 explicit QgsNetworkVisitor( QVector< int > &pointIndexes )
119 : mPoints( pointIndexes ) {}
120
121 void visitNode( const INode &n ) override
122 { Q_UNUSED( n ) }
123
124 void visitData( const IData &d ) override
125 {
126 mPoints.append( d.getIdentifier() );
127 }
128
129 void visitData( std::vector<const IData *> &v ) override
130 { Q_UNUSED( v ) }
131
132 private:
133 QVector< int > &mPoints;
134};
135
137
138std::unique_ptr< SpatialIndex::ISpatialIndex > createVertexSpatialIndex( SpatialIndex::IStorageManager &storageManager )
139{
140 // R-Tree parameters
141 double fillFactor = 0.7;
142 unsigned long indexCapacity = 10;
143 unsigned long leafCapacity = 10;
144 unsigned long dimension = 2;
145 RTree::RTreeVariant variant = RTree::RV_RSTAR;
146
147 // create R-tree
148 SpatialIndex::id_type indexId;
149 std::unique_ptr< SpatialIndex::ISpatialIndex > iRTree( RTree::createNewRTree( storageManager, fillFactor, indexCapacity,
150 leafCapacity, dimension, variant, indexId ) );
151 return iRTree;
152}
153
154int findClosestVertex( const QgsPointXY &point, SpatialIndex::ISpatialIndex *index, double tolerance )
155{
156 QVector< int > matching;
157 QgsNetworkVisitor visitor( matching );
158
159 double pt1[2] = { point.x() - tolerance, point.y() - tolerance },
160 pt2[2] = { point.x() + tolerance, point.y() + tolerance };
161 SpatialIndex::Region searchRegion( pt1, pt2, 2 );
162
163 index->intersectsWithQuery( searchRegion, visitor );
164
165 return matching.empty() ? -1 : matching.at( 0 );
166}
167
168void QgsVectorLayerDirector::makeGraph( QgsGraphBuilderInterface *builder, const QVector< QgsPointXY > &additionalPoints,
169 QVector< QgsPointXY > &snappedPoints, QgsFeedback *feedback ) const
170{
171 long featureCount = mSource->featureCount() * 2;
172 int step = 0;
173
175 ct.setSourceCrs( mSource->sourceCrs() );
176 if ( builder->coordinateTransformationEnabled() )
177 {
178 ct.setDestinationCrs( builder->destinationCrs() );
179 }
180
181 // clear existing snapped points list, and resize to length of provided additional points
182 snappedPoints = QVector< QgsPointXY >( additionalPoints.size(), QgsPointXY( 0.0, 0.0 ) );
183 // tie points = snapped location of specified additional points to network lines
184 QVector< TiePointInfo > additionalTiePoints( additionalPoints.size() );
185
186 // graph's vertices = all vertices in graph, with vertices within builder's tolerance collapsed together
187 QVector< QgsPointXY > graphVertices;
188
189 // spatial index for graph vertices
190 std::unique_ptr< SpatialIndex::IStorageManager > iStorage( StorageManager::createNewMemoryStorageManager() );
191 std::unique_ptr< SpatialIndex::ISpatialIndex > iRTree = createVertexSpatialIndex( *iStorage );
192
193 double tolerance = std::max( builder->topologyTolerance(), 1e-10 );
194 auto findPointWithinTolerance = [&iRTree, tolerance]( const QgsPointXY & point )->int
195 {
196 return findClosestVertex( point, iRTree.get(), tolerance );
197 };
198 auto addPointToIndex = [&iRTree]( const QgsPointXY & point, int index )
199 {
200 double coords[] = {point.x(), point.y()};
201 iRTree->insertData( 0, nullptr, SpatialIndex::Point( coords, 2 ), index );
202 };
203
204 // first iteration - get all nodes from network, and snap additional points to network
206 QgsFeature feature;
207 while ( fit.nextFeature( feature ) )
208 {
209 if ( feedback && feedback->isCanceled() )
210 return;
211
214 mpl = feature.geometry().asMultiPolyline();
216 mpl.push_back( feature.geometry().asPolyline() );
217
218 for ( const QgsPolylineXY &line : std::as_const( mpl ) )
219 {
220 QgsPointXY pt1, pt2;
221 bool isFirstPoint = true;
222 for ( const QgsPointXY &point : line )
223 {
224 pt2 = ct.transform( point );
225
226 int pt2Idx = findPointWithinTolerance( pt2 ) ;
227 if ( pt2Idx == -1 )
228 {
229 // no vertex already exists within tolerance - add to points, and index
230 addPointToIndex( pt2, graphVertices.count() );
231 graphVertices.push_back( pt2 );
232 }
233 else
234 {
235 // vertex already exists within tolerance - use that
236 pt2 = graphVertices.at( pt2Idx );
237 }
238
239 if ( !isFirstPoint )
240 {
241 // check if this line segment is a candidate for being closest to each additional point
242 int i = 0;
243 for ( const QgsPointXY &additionalPoint : additionalPoints )
244 {
245
246 QgsPointXY snappedPoint;
247 double thisSegmentClosestDist = std::numeric_limits<double>::max();
248 if ( pt1 == pt2 )
249 {
250 thisSegmentClosestDist = additionalPoint.sqrDist( pt1 );
251 snappedPoint = pt1;
252 }
253 else
254 {
255 thisSegmentClosestDist = additionalPoint.sqrDistToSegment( pt1.x(), pt1.y(),
256 pt2.x(), pt2.y(), snappedPoint, 0 );
257 }
258
259 if ( thisSegmentClosestDist < additionalTiePoints[ i ].mLength )
260 {
261 // found a closer segment for this additional point
262 TiePointInfo info( i, feature.id(), pt1, pt2 );
263 info.mLength = thisSegmentClosestDist;
264 info.mTiedPoint = snappedPoint;
265
266 additionalTiePoints[ i ] = info;
267 snappedPoints[ i ] = info.mTiedPoint;
268 }
269 i++;
270 }
271 }
272 pt1 = pt2;
273 isFirstPoint = false;
274 }
275 }
276 if ( feedback )
277 feedback->setProgress( 100.0 * static_cast< double >( ++step ) / featureCount );
278 }
279
280 // build a hash of feature ids to tie points which depend on this feature
281 QHash< QgsFeatureId, QList< int > > tiePointNetworkFeatures;
282 int i = 0;
283 for ( TiePointInfo &info : additionalTiePoints )
284 {
285 tiePointNetworkFeatures[ info.mNetworkFeatureId ] << i;
286 i++;
287 }
288
289 // add tied point to graph
290 for ( int i = 0; i < snappedPoints.size(); ++i )
291 {
292 // check index to see if vertex exists within tolerance of tie point
293 const QgsPointXY point = snappedPoints.at( i );
294 int ptIdx = findPointWithinTolerance( point );
295 if ( ptIdx == -1 )
296 {
297 // no vertex already within tolerance, add to index and network vertices
298 addPointToIndex( point, graphVertices.count() );
299 graphVertices.push_back( point );
300 }
301 else
302 {
303 // otherwise snap tie point to vertex
304 snappedPoints[ i ] = graphVertices.at( ptIdx );
305 }
306 }
307 // also need to update tie points - they need to be matched for snapped points
308 for ( int i = 0; i < additionalTiePoints.count(); ++i )
309 {
310 additionalTiePoints[ i ].mTiedPoint = snappedPoints.at( additionalTiePoints.at( i ).additionalPointId );
311 }
312
313
314 // begin graph construction
315
316 // add vertices to graph
317 {
318 int i = 0;
319 for ( const QgsPointXY &point : graphVertices )
320 {
321 builder->addVertex( i, point );
322 i++;
323 }
324 }
325
326 fit = mSource->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( requiredAttributes() ) );
327 while ( fit.nextFeature( feature ) )
328 {
329 if ( feedback && feedback->isCanceled() )
330 return;
331
332 Direction direction = directionForFeature( feature );
333
334 // begin features segments and add arc to the Graph;
337 mpl = feature.geometry().asMultiPolyline();
339 mpl.push_back( feature.geometry().asPolyline() );
340
341 for ( const QgsPolylineXY &line : std::as_const( mpl ) )
342 {
343 QgsPointXY pt1, pt2;
344
345 bool isFirstPoint = true;
346 for ( const QgsPointXY &point : line )
347 {
348 pt2 = ct.transform( point );
349 int pPt2idx = findPointWithinTolerance( pt2 );
350 Q_ASSERT_X( pPt2idx >= 0, "QgsVectorLayerDirectory::makeGraph", "encountered a vertex which was not present in graph" );
351 pt2 = graphVertices.at( pPt2idx );
352
353 if ( !isFirstPoint )
354 {
355 QMap< double, QgsPointXY > pointsOnArc;
356 pointsOnArc[ 0.0 ] = pt1;
357 pointsOnArc[ pt1.sqrDist( pt2 )] = pt2;
358
359 const QList< int > tiePointsForCurrentFeature = tiePointNetworkFeatures.value( feature.id() );
360 for ( int tiePointIdx : tiePointsForCurrentFeature )
361 {
362 const TiePointInfo &t = additionalTiePoints.at( tiePointIdx );
363 if ( t.mFirstPoint == pt1 && t.mLastPoint == pt2 )
364 {
365 pointsOnArc[ pt1.sqrDist( t.mTiedPoint )] = t.mTiedPoint;
366 }
367 }
368
369 QgsPointXY arcPt1;
370 QgsPointXY arcPt2;
371 int pt1idx = -1;
372 int pt2idx = -1;
373 bool isFirstPoint = true;
374 for ( auto arcPointIt = pointsOnArc.constBegin(); arcPointIt != pointsOnArc.constEnd(); ++arcPointIt )
375 {
376 arcPt2 = arcPointIt.value();
377
378 pt2idx = findPointWithinTolerance( arcPt2 );
379 Q_ASSERT_X( pt2idx >= 0, "QgsVectorLayerDirectory::makeGraph", "encountered a vertex which was not present in graph" );
380 arcPt2 = graphVertices.at( pt2idx );
381
382 if ( !isFirstPoint && arcPt1 != arcPt2 )
383 {
384 double distance = builder->distanceArea()->measureLine( arcPt1, arcPt2 );
385 QVector< QVariant > prop;
386 prop.reserve( mStrategies.size() );
387 for ( QgsNetworkStrategy *strategy : mStrategies )
388 {
389 prop.push_back( strategy->cost( distance, feature ) );
390 }
391
392 if ( direction == Direction::DirectionForward ||
393 direction == Direction::DirectionBoth )
394 {
395 builder->addEdge( pt1idx, arcPt1, pt2idx, arcPt2, prop );
396 }
397 if ( direction == Direction::DirectionBackward ||
398 direction == Direction::DirectionBoth )
399 {
400 builder->addEdge( pt2idx, arcPt2, pt1idx, arcPt1, prop );
401 }
402 }
403 pt1idx = pt2idx;
404 arcPt1 = arcPt2;
405 isFirstPoint = false;
406 }
407 }
408 pt1 = pt2;
409 isFirstPoint = false;
410 }
411 }
412 if ( feedback )
413 {
414 feedback->setProgress( 100.0 * static_cast< double >( ++step ) / featureCount );
415 }
416
417 }
418}
@ LineString
LineString.
@ MultiLineString
MultiLineString.
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:338
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
Qgis::WkbType 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 Qgis::WkbType flatType(Qgis::WkbType type) SIP_HOLDGIL
Returns the flat type for a WKB type.
Definition: qgswkbtypes.h:629
#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