QGIS API Documentation  2.99.0-Master (f867b65)
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 #include <QProgressDialog>
32 
33 QgsTINInterpolator::QgsTINInterpolator( const QList<LayerData> &inputData, TINInterpolation interpolation, QgsFeedback *feedback )
34  : QgsInterpolator( inputData )
35  , mTriangulation( nullptr )
36  , mTriangleInterpolator( nullptr )
37  , mIsInitialized( false )
38  , mFeedback( feedback )
39  , mExportTriangulationToFile( false )
40  , mInterpolation( interpolation )
41 {
42 }
43 
45 {
46  delete mTriangulation;
47  delete mTriangleInterpolator;
48 }
49 
50 int QgsTINInterpolator::interpolatePoint( double x, double y, double &result )
51 {
52  if ( !mIsInitialized )
53  {
54  initialize();
55  }
56 
57  if ( !mTriangleInterpolator )
58  {
59  return 1;
60  }
61 
62  QgsPoint r( 0, 0, 0 );
63  if ( !mTriangleInterpolator->calcPoint( x, y, &r ) )
64  {
65  return 2;
66  }
67  result = r.z();
68  return 0;
69 }
70 
71 void QgsTINInterpolator::initialize()
72 {
73  DualEdgeTriangulation *dualEdgeTriangulation = new DualEdgeTriangulation( 100000, nullptr );
74  if ( mInterpolation == CloughTocher )
75  {
77  dec->addTriangulation( dualEdgeTriangulation );
78  mTriangulation = dec;
79  }
80  else
81  {
82  mTriangulation = dualEdgeTriangulation;
83  }
84 
85  //get number of features if we use a progress bar
86  int nFeatures = 0;
87  int nProcessedFeatures = 0;
88  if ( mFeedback )
89  {
90  Q_FOREACH ( const LayerData &layer, mLayerData )
91  {
92  if ( layer.vectorLayer )
93  {
94  nFeatures += layer.vectorLayer->featureCount();
95  }
96  }
97  }
98 
99  QgsFeature f;
100  Q_FOREACH ( const LayerData &layer, mLayerData )
101  {
102  if ( layer.vectorLayer )
103  {
104  QgsAttributeList attList;
105  if ( !layer.zCoordInterpolation )
106  {
107  attList.push_back( layer.interpolationAttribute );
108  }
109 
110  QgsFeatureIterator fit = layer.vectorLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( attList ) );
111 
112  while ( fit.nextFeature( f ) )
113  {
114  if ( mFeedback )
115  {
116  if ( mFeedback->isCanceled() )
117  {
118  break;
119  }
120  mFeedback->setProgress( 100.0 * static_cast< double >( nProcessedFeatures ) / nFeatures );
121  }
122  insertData( &f, layer.zCoordInterpolation, layer.interpolationAttribute, layer.mInputType );
123  ++nProcessedFeatures;
124  }
125  }
126  }
127 
128  if ( mInterpolation == CloughTocher )
129  {
130  CloughTocherInterpolator *ctInterpolator = new CloughTocherInterpolator();
131  NormVecDecorator *dec = dynamic_cast<NormVecDecorator *>( mTriangulation );
132  if ( dec )
133  {
134  dec->estimateFirstDerivatives( mFeedback );
135  ctInterpolator->setTriangulation( dec );
136  dec->setTriangleInterpolator( ctInterpolator );
137  mTriangleInterpolator = ctInterpolator;
138  }
139  }
140  else //linear
141  {
142  mTriangleInterpolator = new LinTriangleInterpolator( dualEdgeTriangulation );
143  }
144  mIsInitialized = true;
145 
146  //debug
147  if ( mExportTriangulationToFile )
148  {
149  dualEdgeTriangulation->saveAsShapefile( mTriangulationFilePath );
150  }
151 }
152 
153 int QgsTINInterpolator::insertData( QgsFeature *f, bool zCoord, int attr, InputType type )
154 {
155  if ( !f )
156  {
157  return 1;
158  }
159 
160  QgsGeometry g = f->geometry();
161  {
162  if ( g.isNull() )
163  {
164  return 2;
165  }
166  }
167 
168  //check attribute value
169  double attributeValue = 0;
170  bool attributeConversionOk = false;
171  if ( !zCoord )
172  {
173  QVariant attributeVariant = f->attribute( attr );
174  if ( !attributeVariant.isValid() ) //attribute not found, something must be wrong (e.g. NULL value)
175  {
176  return 3;
177  }
178  attributeValue = attributeVariant.toDouble( &attributeConversionOk );
179  if ( !attributeConversionOk || qIsNaN( attributeValue ) ) //don't consider vertices with attributes like 'nan' for the interpolation
180  {
181  return 4;
182  }
183  }
184 
185  //parse WKB. It is ugly, but we cannot use the methods with QgsPointXY because they don't contain z-values for 25D types
186  bool hasZValue = false;
187  double x, y, z;
188  QByteArray wkb( g.exportToWkb() );
189  QgsConstWkbPtr currentWkbPtr( wkb );
190  currentWkbPtr.readHeader();
191  //maybe a structure or break line
192  Line3D *line = nullptr;
193 
194  QgsWkbTypes::Type wkbType = g.wkbType();
195  switch ( wkbType )
196  {
198  hasZValue = true;
199  FALLTHROUGH;
200  case QgsWkbTypes::Point:
201  {
202  currentWkbPtr >> x >> y;
203  if ( zCoord && hasZValue )
204  {
205  currentWkbPtr >> z;
206  }
207  else
208  {
209  z = attributeValue;
210  }
211  QgsPoint *point = new QgsPoint( x, y, z );
212  if ( mTriangulation->addPoint( point ) == -100 )
213  {
214  return -1;
215  }
216  break;
217  }
219  hasZValue = true;
220  FALLTHROUGH;
222  {
223  int nPoints;
224  currentWkbPtr >> nPoints;
225  for ( int index = 0; index < nPoints; ++index )
226  {
227  currentWkbPtr.readHeader();
228  currentWkbPtr >> x >> y;
229  if ( hasZValue ) //skip z-coordinate for 25D geometries
230  {
231  currentWkbPtr >> z;
232  }
233  else
234  {
235  z = attributeValue;
236  }
237  }
238  break;
239  }
241  hasZValue = true;
242  FALLTHROUGH;
244  {
245  if ( type != POINTS )
246  {
247  line = new Line3D();
248  }
249  int nPoints;
250  currentWkbPtr >> nPoints;
251  for ( int index = 0; index < nPoints; ++index )
252  {
253  currentWkbPtr >> x >> y;
254  if ( zCoord && hasZValue ) //skip z-coordinate for 25D geometries
255  {
256  currentWkbPtr >> z;
257  }
258  else
259  {
260  z = attributeValue;
261  }
262 
263  if ( type == POINTS )
264  {
265  //todo: handle error code -100
266  mTriangulation->addPoint( new QgsPoint( x, y, z ) );
267  }
268  else
269  {
270  line->insertPoint( new QgsPoint( x, y, z ) );
271  }
272  }
273 
274  if ( type != POINTS )
275  {
276  mTriangulation->addLine( line, type == BREAK_LINES );
277  }
278  break;
279  }
281  hasZValue = true;
282  FALLTHROUGH;
284  {
285  int nLines;
286  currentWkbPtr >> nLines;
287  for ( int index = 0; index < nLines; ++index )
288  {
289  if ( type != POINTS )
290  {
291  line = new Line3D();
292  }
293  int nPoints;
294  currentWkbPtr >> nPoints;
295  for ( int index2 = 0; index2 < nPoints; ++index2 )
296  {
297  currentWkbPtr >> x >> y;
298  if ( hasZValue ) //skip z-coordinate for 25D geometries
299  {
300  currentWkbPtr >> z;
301  }
302  else
303  {
304  z = attributeValue;
305  }
306 
307  if ( type == POINTS )
308  {
309  //todo: handle error code -100
310  mTriangulation->addPoint( new QgsPoint( x, y, z ) );
311  }
312  else
313  {
314  line->insertPoint( new QgsPoint( x, y, z ) );
315  }
316  }
317  if ( type != POINTS )
318  {
319  mTriangulation->addLine( line, type == BREAK_LINES );
320  }
321  }
322  break;
323  }
325  hasZValue = true;
326  FALLTHROUGH;
328  {
329  int nRings;
330  currentWkbPtr >> nRings;
331  for ( int index = 0; index < nRings; ++index )
332  {
333  if ( type != POINTS )
334  {
335  line = new Line3D();
336  }
337 
338  int nPoints;
339  currentWkbPtr >> nPoints;
340  for ( int index2 = 0; index2 < nPoints; ++index2 )
341  {
342  currentWkbPtr >> x >> y;
343  if ( hasZValue ) //skip z-coordinate for 25D geometries
344  {
345  currentWkbPtr >> z;
346  }
347  else
348  {
349  z = attributeValue;
350  }
351  if ( type == POINTS )
352  {
353  //todo: handle error code -100
354  mTriangulation->addPoint( new QgsPoint( x, y, z ) );
355  }
356  else
357  {
358  line->insertPoint( new QgsPoint( x, y, z ) );
359  }
360  }
361 
362  if ( type != POINTS )
363  {
364  mTriangulation->addLine( line, type == BREAK_LINES );
365  }
366  }
367  break;
368  }
369 
371  hasZValue = true;
372  FALLTHROUGH;
374  {
375  int nPolys;
376  currentWkbPtr >> nPolys;
377  for ( int index = 0; index < nPolys; ++index )
378  {
379  currentWkbPtr.readHeader();
380  int nRings;
381  currentWkbPtr >> nRings;
382  for ( int index2 = 0; index2 < nRings; ++index2 )
383  {
384  if ( type != POINTS )
385  {
386  line = new Line3D();
387  }
388  int nPoints;
389  currentWkbPtr >> nPoints;
390  for ( int index3 = 0; index3 < nPoints; ++index3 )
391  {
392  currentWkbPtr >> x >> y;
393  if ( hasZValue ) //skip z-coordinate for 25D geometries
394  {
395  currentWkbPtr >> z;
396  }
397  else
398  {
399  z = attributeValue;
400  }
401  if ( type == POINTS )
402  {
403  //todo: handle error code -100
404  mTriangulation->addPoint( new QgsPoint( x, y, z ) );
405  }
406  else
407  {
408  line->insertPoint( new QgsPoint( x, y, z ) );
409  }
410  }
411  if ( type != POINTS )
412  {
413  mTriangulation->addLine( line, type == BREAK_LINES );
414  }
415  }
416  }
417  break;
418  }
419  default:
420  //should not happen...
421  break;
422  }
423 
424  return 0;
425 }
426 
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:82
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.cc: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:72
virtual void setTriangulation(NormVecDecorator *tin)
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:96
This class represents a line.
Definition: Line3D.h:26
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:61
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:43
Type
The WKB type describes the number of dimensions a geometry has.
Definition: qgswkbtypes.h:66
InputType
Describes the type of input data.
#define FALLTHROUGH
Definition: qgis.h:415
QByteArray exportToWkb() const
Export the geometry to WKB.
long featureCount(const QString &legendKey) const
Number of features rendered with specified legend key.
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:36
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:63
int interpolatePoint(double x, double y, double &result) override
Calculates interpolation value for map coordinates x, y.
This is an implementation of a Clough-Tocher interpolator based on a triangular tessellation.
double z
Definition: qgspoint.h:42
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 ...
QVariant attribute(const QString &name) const
Lookup attribute value from attribute name.
Definition: qgsfeature.cpp:255
virtual bool saveAsShapefile(const QString &fileName) const override
Saves the triangulation as a (line) shapefile.