QGIS API Documentation  3.0.2-Girona (307d082)
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.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 = DBL_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().setSubsetOfAttributes( QgsAttributeList() ) );
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 = DBL_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 );
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  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 }
QgsFeatureId id
Definition: qgsfeature.h:71
Wrapper for iterator of features from vector data provider or vector layer.
QgsNetworkStrategy defines strategy used for calculation of the edge cost. For example it can take in...
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.
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
Determine interface for creating a graph.
TiePointInfo(int additionalPointId, QgsFeatureId featureId, const QgsPointXY &start, const QgsPointXY &end)
virtual void addEdge(int pt1id, const QgsPointXY &pt1, int pt2id, const QgsPointXY &pt2, const QVector< QVariant > &strategies)
Add edge to the graph.
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:62
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.cpp:68
QgsMultiPolylineXY asMultiPolyline() const
Returns contents of the geometry as a multi linestring if wkbType is WKBMultiLineString, otherwise an empty list.
Base class for feedback objects to be used for cancelation 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:136
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)...
QgsVectorLayerDirector(QgsFeatureSource *source, int directionFieldId, const QString &directDirectionValue, const QString &reverseDirectionValue, const QString &bothDirectionValue, const Direction defaultDirection)
Default constructor.
QgsGeometry geometry() const
Returns the geometry associated with this feature.
Definition: qgsfeature.cpp:101
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.
virtual void addVertex(int id, const QgsPointXY &pt)
Add vertex to the graph.
QgsPointXY transform(const QgsPointXY &point, TransformDirection direction=ForwardTransform) const
Transform the point from the source CRS to the destination CRS.
QList< QgsNetworkStrategy * > mStrategies
Class for doing transforms between two map coordinate systems.
qint64 QgsFeatureId
Definition: qgsfeature.h:37
double topologyTolerance()
Returns topology tolerance.
QgsPolylineXY asPolyline() const
Returns contents of the geometry as a polyline if wkbType is WKBLineString, otherwise an empty list...
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:255
static Type flatType(Type type)
Returns the flat type for a WKB type.
Definition: qgswkbtypes.h:427
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.