QGIS API Documentation  3.8.0-Zanzibar (11aff65)
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
6  Email : [email protected]
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 
21 #include "qgsvectorlayerdirector.h"
23 
24 #include "qgsfeatureiterator.h"
25 #include "qgsfeaturesource.h"
26 #include "qgsvectordataprovider.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 
37 using namespace SpatialIndex;
38 
40 {
41  TiePointInfo() = default;
42  TiePointInfo( int additionalPointId, QgsFeatureId featureId, const QgsPointXY &start, const QgsPointXY &end )
43  : additionalPointId( additionalPointId )
44  , mNetworkFeatureId( featureId )
45  , mFirstPoint( start )
46  , mLastPoint( end )
47  {}
48 
49  int additionalPointId = -1;
51  double mLength = std::numeric_limits<double>::max();
52  QgsFeatureId mNetworkFeatureId = -1;
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 
77 QgsAttributeList 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 attrs.toList();
89 }
90 
91 QgsVectorLayerDirector::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 
116 class 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 
139 std::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 
155 int 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 
169 void 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
206  QgsFeatureIterator fit = mSource->getFeatures( QgsFeatureRequest().setNoAttributes() );
207  QgsFeature feature;
208  while ( fit.nextFeature( feature ) )
209  {
210  if ( feedback && feedback->isCanceled() )
211  return;
212 
213  QgsMultiPolylineXY mpl;
215  mpl = feature.geometry().asMultiPolyline();
216  else if ( QgsWkbTypes::flatType( feature.geometry().wkbType() ) == QgsWkbTypes::LineString )
217  mpl.push_back( feature.geometry().asPolyline() );
218 
219  for ( const QgsPolylineXY &line : qgis::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;
336  QgsMultiPolylineXY mpl;
338  mpl = feature.geometry().asMultiPolyline();
339  else if ( QgsWkbTypes::flatType( feature.geometry().wkbType() ) == QgsWkbTypes::LineString )
340  mpl.push_back( feature.geometry().asPolyline() );
341 
342  for ( const QgsPolylineXY &line : qgis::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 }
QgsFeatureId id
Definition: qgsfeature.h:64
Wrapper for iterator of features from vector data provider or vector layer.
QgsNetworkStrategy defines strategy used for calculation of the edge cost.
QgsVectorLayerDirector(QgsFeatureSource *source, int directionFieldId, const QString &directDirectionValue, const QString &reverseDirectionValue, const QString &bothDirectionValue, Direction defaultDirection)
Default constructor.
Direction
Edge direction Edge can be one-way with direct flow (one can move only from the start point to the en...
int findClosestVertex(const QgsPointXY &point, SpatialIndex::ISpatialIndex *index, double tolerance)
void makeGraph(QgsGraphBuilderInterface *builder, const QVector< QgsPointXY > &additionalPoints, QVector< QgsPointXY > &snappedPoints, QgsFeedback *feedback=nullptr) const override
Make a graph using QgsGraphBuilder.
QgsPointXY transform(const QgsPointXY &point, TransformDirection direction=ForwardTransform) const SIP_THROW(QgsCsException)
Transform the point from the source CRS to the destination CRS.
QgsWkbTypes::Type wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
double y
Definition: qgspointxy.h:48
bool coordinateTransformationEnabled()
Returns coordinate transformation enabled.
A class to represent a 2D point.
Definition: qgspointxy.h:43
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
qint64 QgsFeatureId
Definition: qgsfeatureid.h:25
Determine interface for creating a graph.
TiePointInfo(int additionalPointId, QgsFeatureId featureId, const QgsPointXY &start, const QgsPointXY &end)
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:55
QVector< QgsPolylineXY > QgsMultiPolylineXY
A collection of QgsPolylines that share a common collection of attributes.
Definition: qgsgeometry.h:83
double sqrDist(double x, double y) const
Returns the squared distance between this point a specified x, y coordinate.
Definition: qgspointxy.h:169
virtual void addVertex(int id, const QgsPointXY &pt)
Add vertex to the graph.
QgsMultiPolylineXY asMultiPolyline() const
Returns the contents of the geometry as a multi-linestring.
Base class for feedback objects to be used for cancellation of something running in a worker thread...
Definition: qgsfeedback.h:44
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.
Definition: qgspointxy.cpp:78
void setDestinationCrs(const QgsCoordinateReferenceSystem &crs)
Sets the destination coordinate reference system.
This class wraps a request for features to a vector layer (or directly its vector data provider)...
virtual void addEdge(int pt1id, const QgsPointXY &pt1, int pt2id, const QgsPointXY &pt2, const QVector< QVariant > &strategies)
Add edge to the graph.
double x
Definition: qgspointxy.h:47
std::unique_ptr< SpatialIndex::ISpatialIndex > createVertexSpatialIndex(SpatialIndex::IStorageManager &storageManager)
QVector< QgsPointXY > QgsPolylineXY
Polyline as represented as a vector of two-dimensional points.
Definition: qgsgeometry.h:49
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
An interface for objects which provide features via a getFeatures method.
QList< QgsNetworkStrategy * > mStrategies
Class for doing transforms between two map coordinate systems.
double topologyTolerance()
Returns topology tolerance.
QgsPolylineXY asPolyline() const
Returns the contents of the geometry as a polyline.
QgsGeometry geometry
Definition: qgsfeature.h:67
QList< int > QgsAttributeList
Definition: qgsfield.h:27
bool nextFeature(QgsFeature &f)
QgsCoordinateReferenceSystem destinationCrs() const
Returns destinaltion CRS.
QVariant attribute(const QString &name) const
Lookup attribute value from attribute name.
Definition: qgsfeature.cpp:262
static Type flatType(Type type)
Returns the flat type for a WKB type.
Definition: qgswkbtypes.h:430
void setSourceCrs(const QgsCoordinateReferenceSystem &crs)
Sets the source coordinate reference system.
QgsDistanceArea * distanceArea()
Returns measurement tool.
double measureLine(const QVector< QgsPointXY > &points) const
Measures the length of a line with multiple segments.
QString name() const override
Returns director name.