QGIS API Documentation  2.99.0-Master (ae4d26a)
qgstininterpolator.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgstininterpolator.cpp
3  ----------------------
4  begin : March 10, 2008
5  copyright : (C) 2008 by Marco Hugentobler
6  email : marco dot hugentobler at karto dot baug dot ethz dot ch
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #include "qgstininterpolator.h"
19 #include "qgsfeatureiterator.h"
21 #include "DualEdgeTriangulation.h"
22 #include "NormVecDecorator.h"
24 #include "Line3D.h"
25 #include "qgspoint.h"
26 #include "qgsfeature.h"
27 #include "qgsgeometry.h"
28 #include "qgsvectorlayer.h"
29 #include "qgswkbptr.h"
30 #include "qgsfeedback.h"
31 
32 QgsTINInterpolator::QgsTINInterpolator( const QList<LayerData> &inputData, TINInterpolation interpolation, QgsFeedback *feedback )
33  : QgsInterpolator( inputData )
34  , mIsInitialized( false )
35  , mFeedback( feedback )
36  , mInterpolation( interpolation )
37 {
38 }
39 
41 {
42  delete mTriangulation;
43  delete mTriangleInterpolator;
44 }
45 
46 int QgsTINInterpolator::interpolatePoint( double x, double y, double &result )
47 {
48  if ( !mIsInitialized )
49  {
50  initialize();
51  }
52 
53  if ( !mTriangleInterpolator )
54  {
55  return 1;
56  }
57 
58  QgsPoint r( 0, 0, 0 );
59  if ( !mTriangleInterpolator->calcPoint( x, y, &r ) )
60  {
61  return 2;
62  }
63  result = r.z();
64  return 0;
65 }
66 
68 {
70 }
71 
73 {
74  mTriangulationSink = sink;
75 }
76 
77 void QgsTINInterpolator::initialize()
78 {
79  DualEdgeTriangulation *dualEdgeTriangulation = new DualEdgeTriangulation( 100000, nullptr );
80  if ( mInterpolation == CloughTocher )
81  {
83  dec->addTriangulation( dualEdgeTriangulation );
84  mTriangulation = dec;
85  }
86  else
87  {
88  mTriangulation = dualEdgeTriangulation;
89  }
90 
91  //get number of features if we use a progress bar
92  int nFeatures = 0;
93  int nProcessedFeatures = 0;
94  if ( mFeedback )
95  {
96  Q_FOREACH ( const LayerData &layer, mLayerData )
97  {
98  if ( layer.vectorLayer )
99  {
100  nFeatures += layer.vectorLayer->featureCount();
101  }
102  }
103  }
104 
105  QgsFeature f;
106  Q_FOREACH ( const LayerData &layer, mLayerData )
107  {
108  if ( layer.vectorLayer )
109  {
110  QgsAttributeList attList;
111  if ( !layer.zCoordInterpolation )
112  {
113  attList.push_back( layer.interpolationAttribute );
114  }
115 
116  QgsFeatureIterator fit = layer.vectorLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( attList ) );
117 
118  while ( fit.nextFeature( f ) )
119  {
120  if ( mFeedback )
121  {
122  if ( mFeedback->isCanceled() )
123  {
124  break;
125  }
126  if ( nFeatures > 0 )
127  mFeedback->setProgress( 100.0 * static_cast< double >( nProcessedFeatures ) / nFeatures );
128  }
129  insertData( &f, layer.zCoordInterpolation, layer.interpolationAttribute, layer.mInputType );
130  ++nProcessedFeatures;
131  }
132  }
133  }
134 
135  if ( mInterpolation == CloughTocher )
136  {
137  CloughTocherInterpolator *ctInterpolator = new CloughTocherInterpolator();
138  NormVecDecorator *dec = dynamic_cast<NormVecDecorator *>( mTriangulation );
139  if ( dec )
140  {
141  dec->estimateFirstDerivatives( mFeedback );
142  ctInterpolator->setTriangulation( dec );
143  dec->setTriangleInterpolator( ctInterpolator );
144  mTriangleInterpolator = ctInterpolator;
145  }
146  }
147  else //linear
148  {
149  mTriangleInterpolator = new LinTriangleInterpolator( dualEdgeTriangulation );
150  }
151  mIsInitialized = true;
152 
153  //debug
154  if ( mTriangulationSink )
155  {
156  dualEdgeTriangulation->saveTriangulation( mTriangulationSink, mFeedback );
157  }
158 }
159 
160 int QgsTINInterpolator::insertData( QgsFeature *f, bool zCoord, int attr, InputType type )
161 {
162  if ( !f )
163  {
164  return 1;
165  }
166 
167  QgsGeometry g = f->geometry();
168  {
169  if ( g.isNull() )
170  {
171  return 2;
172  }
173  }
174 
175  //check attribute value
176  double attributeValue = 0;
177  bool attributeConversionOk = false;
178  if ( !zCoord )
179  {
180  QVariant attributeVariant = f->attribute( attr );
181  if ( !attributeVariant.isValid() ) //attribute not found, something must be wrong (e.g. NULL value)
182  {
183  return 3;
184  }
185  attributeValue = attributeVariant.toDouble( &attributeConversionOk );
186  if ( !attributeConversionOk || std::isnan( attributeValue ) ) //don't consider vertices with attributes like 'nan' for the interpolation
187  {
188  return 4;
189  }
190  }
191 
192  //parse WKB. It is ugly, but we cannot use the methods with QgsPointXY because they don't contain z-values for 25D types
193  bool hasZValue = false;
194  double x, y, z;
195  QByteArray wkb( g.exportToWkb() );
196  QgsConstWkbPtr currentWkbPtr( wkb );
197  currentWkbPtr.readHeader();
198  //maybe a structure or break line
199  Line3D *line = nullptr;
200 
201  QgsWkbTypes::Type wkbType = g.wkbType();
202  switch ( wkbType )
203  {
205  hasZValue = true;
206  FALLTHROUGH;
207  case QgsWkbTypes::Point:
208  {
209  currentWkbPtr >> x >> y;
210  if ( zCoord && hasZValue )
211  {
212  currentWkbPtr >> z;
213  }
214  else
215  {
216  z = attributeValue;
217  }
218  QgsPoint *point = new QgsPoint( x, y, z );
219  if ( mTriangulation->addPoint( point ) == -100 )
220  {
221  return -1;
222  }
223  break;
224  }
226  hasZValue = true;
227  FALLTHROUGH;
229  {
230  int nPoints;
231  currentWkbPtr >> nPoints;
232  for ( int index = 0; index < nPoints; ++index )
233  {
234  currentWkbPtr.readHeader();
235  currentWkbPtr >> x >> y;
236  if ( hasZValue ) //skip z-coordinate for 25D geometries
237  {
238  currentWkbPtr >> z;
239  }
240  else
241  {
242  z = attributeValue;
243  }
244  }
245  break;
246  }
248  hasZValue = true;
249  FALLTHROUGH;
251  {
252  if ( type != POINTS )
253  {
254  line = new Line3D();
255  }
256  int nPoints;
257  currentWkbPtr >> nPoints;
258  for ( int index = 0; index < nPoints; ++index )
259  {
260  currentWkbPtr >> x >> y;
261  if ( zCoord && hasZValue ) //skip z-coordinate for 25D geometries
262  {
263  currentWkbPtr >> z;
264  }
265  else
266  {
267  z = attributeValue;
268  }
269 
270  if ( type == POINTS )
271  {
272  //todo: handle error code -100
273  mTriangulation->addPoint( new QgsPoint( x, y, z ) );
274  }
275  else
276  {
277  line->insertPoint( new QgsPoint( x, y, z ) );
278  }
279  }
280 
281  if ( type != POINTS )
282  {
283  mTriangulation->addLine( line, type == BREAK_LINES );
284  }
285  break;
286  }
288  hasZValue = true;
289  FALLTHROUGH;
291  {
292  int nLines;
293  currentWkbPtr >> nLines;
294  for ( int index = 0; index < nLines; ++index )
295  {
296  if ( type != POINTS )
297  {
298  line = new Line3D();
299  }
300  int nPoints;
301  currentWkbPtr >> nPoints;
302  for ( int index2 = 0; index2 < nPoints; ++index2 )
303  {
304  currentWkbPtr >> x >> y;
305  if ( hasZValue ) //skip z-coordinate for 25D geometries
306  {
307  currentWkbPtr >> z;
308  }
309  else
310  {
311  z = attributeValue;
312  }
313 
314  if ( type == POINTS )
315  {
316  //todo: handle error code -100
317  mTriangulation->addPoint( new QgsPoint( x, y, z ) );
318  }
319  else
320  {
321  line->insertPoint( new QgsPoint( x, y, z ) );
322  }
323  }
324  if ( type != POINTS )
325  {
326  mTriangulation->addLine( line, type == BREAK_LINES );
327  }
328  }
329  break;
330  }
332  hasZValue = true;
333  FALLTHROUGH;
335  {
336  int nRings;
337  currentWkbPtr >> nRings;
338  for ( int index = 0; index < nRings; ++index )
339  {
340  if ( type != POINTS )
341  {
342  line = new Line3D();
343  }
344 
345  int nPoints;
346  currentWkbPtr >> nPoints;
347  for ( int index2 = 0; index2 < nPoints; ++index2 )
348  {
349  currentWkbPtr >> x >> y;
350  if ( hasZValue ) //skip z-coordinate for 25D geometries
351  {
352  currentWkbPtr >> z;
353  }
354  else
355  {
356  z = attributeValue;
357  }
358  if ( type == POINTS )
359  {
360  //todo: handle error code -100
361  mTriangulation->addPoint( new QgsPoint( x, y, z ) );
362  }
363  else
364  {
365  line->insertPoint( new QgsPoint( x, y, z ) );
366  }
367  }
368 
369  if ( type != POINTS )
370  {
371  mTriangulation->addLine( line, type == BREAK_LINES );
372  }
373  }
374  break;
375  }
376 
378  hasZValue = true;
379  FALLTHROUGH;
381  {
382  int nPolys;
383  currentWkbPtr >> nPolys;
384  for ( int index = 0; index < nPolys; ++index )
385  {
386  currentWkbPtr.readHeader();
387  int nRings;
388  currentWkbPtr >> nRings;
389  for ( int index2 = 0; index2 < nRings; ++index2 )
390  {
391  if ( type != POINTS )
392  {
393  line = new Line3D();
394  }
395  int nPoints;
396  currentWkbPtr >> nPoints;
397  for ( int index3 = 0; index3 < nPoints; ++index3 )
398  {
399  currentWkbPtr >> x >> y;
400  if ( hasZValue ) //skip z-coordinate for 25D geometries
401  {
402  currentWkbPtr >> z;
403  }
404  else
405  {
406  z = attributeValue;
407  }
408  if ( type == POINTS )
409  {
410  //todo: handle error code -100
411  mTriangulation->addPoint( new QgsPoint( x, y, z ) );
412  }
413  else
414  {
415  line->insertPoint( new QgsPoint( x, y, z ) );
416  }
417  }
418  if ( type != POINTS )
419  {
420  mTriangulation->addLine( line, type == BREAK_LINES );
421  }
422  }
423  }
424  break;
425  }
426  default:
427  //should not happen...
428  break;
429  }
430 
431  return 0;
432 }
433 
QgsTINInterpolator(const QList< QgsInterpolator::LayerData > &inputData, TINInterpolation interpolation=Linear, QgsFeedback *feedback=nullptr)
Constructor for QgsTINInterpolator.
Decorator class which adds the functionality of estimating normals at the data points.
Wrapper for iterator of features from vector data provider or vector layer.
QList< LayerData > mLayerData
virtual int addPoint(QgsPoint *p)=0
Adds a point to the triangulation Ownership is transferred to this class.
Interface class for interpolations.
QgsVectorLayer * vectorLayer
bool isNull() const
Returns true if the geometry is null (ie, contains no underlying geometry accessible via geometry() )...
virtual void addTriangulation(Triangulation *t)
Adds an association to a triangulation.
Definition: TriDecorator.h:88
LinTriangleInterpolator is a class which interpolates linearly on a triangulation.
void insertPoint(QgsPoint *p)
Inserts a node behind the current position and sets the current position to this new node...
Definition: Line3D.cpp:49
QgsWkbTypes::Type wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
An interface for objects which accept features via addFeature(s) methods.
Container of fields for a vector layer.
Definition: qgsfields.h:42
virtual void setTriangulation(NormVecDecorator *tin)
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:94
This class represents a line.
Definition: Line3D.h:31
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:62
DualEdgeTriangulation is an implementation of a triangulation class based on the dual edge data struc...
Base class for feedback objects to be used for cancelation of something running in a worker thread...
Definition: qgsfeedback.h:44
Type
The WKB type describes the number of dimensions a geometry has.
Definition: qgswkbtypes.h:67
InputType
Describes the type of input data.
#define FALLTHROUGH
Definition: qgis.h:512
QByteArray exportToWkb() const
Export the geometry to WKB.
long featureCount(const QString &legendKey) const
Number of features rendered with specified legend key.
virtual bool saveTriangulation(QgsFeatureSink *sink, QgsFeedback *feedback=nullptr) const override
Saves the triangulation features to a feature sink.
This class wraps a request for features to a vector layer (or directly its vector data provider)...
virtual void addLine(Line3D *line, bool breakline)=0
Adds a line (e.g.
QgsInterpolator::InputType mInputType
QgsGeometry geometry() const
Returns the geometry associated with this feature.
Definition: qgsfeature.cpp:101
void setTriangleInterpolator(TriangleInterpolator *inter) override
Sets an interpolator.
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
static QgsFields triangulationFields()
Returns the fields output by features when saving the triangulation.
QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const override
Query the layer for features specified in request.
A layer together with the information about interpolation attribute / z-coordinate interpolation and ...
virtual bool calcPoint(double x, double y, QgsPoint *result SIP_OUT)=0
Performs a linear interpolation in a triangle and assigns the x-,y- and z-coordinates to point...
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
int interpolatePoint(double x, double y, double &result) override
Calculates interpolation value for map coordinates x, y.
static QgsFields triangulationFields()
Returns the fields output by features when calling saveTriangulation().
This is an implementation of a Clough-Tocher interpolator based on a triangular tessellation.
double z
Definition: qgspoint.h:43
QList< int > QgsAttributeList
Definition: qgsfield.h:27
bool nextFeature(QgsFeature &f)
bool estimateFirstDerivatives(QgsFeedback *feedback=nullptr)
This method adds the functionality of estimating normals at the data points. Return true in the case ...
void setTriangulationSink(QgsFeatureSink *sink)
Sets the optional sink for saving the triangulation features.
QVariant attribute(const QString &name) const
Lookup attribute value from attribute name.
Definition: qgsfeature.cpp:255