QGIS API Documentation 3.40.0-Bratislava (b56115d8743)
Loading...
Searching...
No Matches
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#include "qgslogger.h"
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 {
99 }
100 else if ( str == mDirectDirectionValue )
101 {
103 }
104 else if ( str == mReverseDirectionValue )
105 {
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 = 0;
385 try
386 {
387 distance = builder->distanceArea()->measureLine( arcPt1, arcPt2 );
388 }
389 catch ( QgsCsException & )
390 {
391 // TODO report errors to user
392 QgsDebugError( QStringLiteral( "An error occurred while calculating length" ) );
393 }
394
395 QVector< QVariant > prop;
396 prop.reserve( mStrategies.size() );
397 for ( QgsNetworkStrategy *strategy : mStrategies )
398 {
399 prop.push_back( strategy->cost( distance, feature ) );
400 }
401
402 if ( direction == Direction::DirectionForward ||
403 direction == Direction::DirectionBoth )
404 {
405 builder->addEdge( pt1idx, arcPt1, pt2idx, arcPt2, prop );
406 }
407 if ( direction == Direction::DirectionBackward ||
408 direction == Direction::DirectionBoth )
409 {
410 builder->addEdge( pt2idx, arcPt2, pt1idx, arcPt1, prop );
411 }
412 }
413 pt1idx = pt2idx;
414 arcPt1 = arcPt2;
415 isFirstPoint = false;
416 }
417 }
418 pt1 = pt2;
419 isFirstPoint = false;
420 }
421 }
422 if ( feedback )
423 {
424 feedback->setProgress( 100.0 * static_cast< double >( ++step ) / featureCount );
425 }
426
427 }
428}
@ 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
Transform the point from the source CRS to the destination CRS.
Custom exception class for Coordinate Reference System related exceptions.
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)
Fetch next feature and stores in f, returns true on success.
This class wraps a request for features to a vector layer (or directly its vector data provider).
QgsFeatureRequest & setSubsetOfAttributes(const QgsAttributeList &attrs)
Set a subset of attributes that will be fetched.
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:58
QgsFeatureId id
Definition qgsfeature.h:66
QgsGeometry geometry
Definition qgsfeature.h:69
Q_INVOKABLE QVariant attribute(const QString &name) const
Lookup attribute value by attribute name.
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition qgsfeedback.h:44
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
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.
Qgis::WkbType wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
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:60
double sqrDist(double x, double y) const
Returns the squared distance between this point a specified x, y coordinate.
Definition qgspointxy.h:186
double y
Definition qgspointxy.h:64
double sqrDistToSegment(double x1, double y1, double x2, double y2, QgsPointXY &minDistPoint, double epsilon=DEFAULT_SEGMENT_EPSILON) const
Returns the minimum distance between this point and a segment.
double x
Definition qgspointxy.h:63
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...
@ DirectionBackward
One-way reversed.
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)
Returns the flat type for a WKB type.
#define str(x)
Definition qgis.cpp:38
qint64 QgsFeatureId
64 bit feature ids negative numbers are used for uncommitted/newly added features
QList< int > QgsAttributeList
Definition qgsfield.h:27
QVector< QgsPolylineXY > QgsMultiPolylineXY
A collection of QgsPolylines that share a common collection of attributes.
Definition qgsgeometry.h:84
QVector< QgsPointXY > QgsPolylineXY
Polyline as represented as a vector of two-dimensional points.
Definition qgsgeometry.h:62
#define QgsDebugError(str)
Definition qgslogger.h:38
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