QGIS API Documentation  2.99.0-Master (6a61179)
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 "Point3D.h"
26 #include "qgsfeature.h"
27 #include "qgsgeometry.h"
28 #include "qgsvectorlayer.h"
29 #include "qgswkbptr.h"
30 #include <QProgressDialog>
31 
32 QgsTINInterpolator::QgsTINInterpolator( const QList<LayerData>& inputData, TIN_INTERPOLATION interpolation, bool showProgressDialog )
33  : QgsInterpolator( inputData )
34  , mTriangulation( nullptr )
35  , mTriangleInterpolator( nullptr )
36  , mIsInitialized( false )
37  , mShowProgressDialog( showProgressDialog )
38  , mExportTriangulationToFile( false )
39  , mInterpolation( interpolation )
40 {
41 }
42 
44 {
45  delete mTriangulation;
46  delete mTriangleInterpolator;
47 }
48 
49 int QgsTINInterpolator::interpolatePoint( double x, double y, double& result )
50 {
51  if ( !mIsInitialized )
52  {
53  initialize();
54  }
55 
56  if ( !mTriangleInterpolator )
57  {
58  return 1;
59  }
60 
61  Point3D r;
62  if ( !mTriangleInterpolator->calcPoint( x, y, &r ) )
63  {
64  return 2;
65  }
66  result = r.getZ();
67  return 0;
68 }
69 
70 void QgsTINInterpolator::initialize()
71 {
72  DualEdgeTriangulation* theDualEdgeTriangulation = new DualEdgeTriangulation( 100000, nullptr );
73  if ( mInterpolation == CloughTocher )
74  {
76  dec->addTriangulation( theDualEdgeTriangulation );
77  mTriangulation = dec;
78  }
79  else
80  {
81  mTriangulation = theDualEdgeTriangulation;
82  }
83 
84  //get number of features if we use a progress bar
85  int nFeatures = 0;
86  int nProcessedFeatures = 0;
87  if ( mShowProgressDialog )
88  {
89  Q_FOREACH ( const LayerData& layer, mLayerData )
90  {
91  if ( layer.vectorLayer )
92  {
93  nFeatures += layer.vectorLayer->featureCount();
94  }
95  }
96  }
97 
98  QProgressDialog* theProgressDialog = nullptr;
99  if ( mShowProgressDialog )
100  {
101  theProgressDialog = new QProgressDialog( QObject::tr( "Building triangulation..." ), QObject::tr( "Abort" ), 0, nFeatures, nullptr );
102  theProgressDialog->setWindowModality( Qt::WindowModal );
103  }
104 
105 
106  QgsFeature f;
107  Q_FOREACH ( const LayerData& layer, mLayerData )
108  {
109  if ( layer.vectorLayer )
110  {
111  QgsAttributeList attList;
112  if ( !layer.zCoordInterpolation )
113  {
114  attList.push_back( layer.interpolationAttribute );
115  }
116 
117  QgsFeatureIterator fit = layer.vectorLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( attList ) );
118 
119  while ( fit.nextFeature( f ) )
120  {
121  if ( mShowProgressDialog )
122  {
123  if ( theProgressDialog->wasCanceled() )
124  {
125  break;
126  }
127  theProgressDialog->setValue( nProcessedFeatures );
128  }
129  insertData( &f, layer.zCoordInterpolation, layer.interpolationAttribute, layer.mInputType );
130  ++nProcessedFeatures;
131  }
132  }
133  }
134 
135  delete theProgressDialog;
136 
137  if ( mInterpolation == CloughTocher )
138  {
139  CloughTocherInterpolator* ctInterpolator = new CloughTocherInterpolator();
140  NormVecDecorator* dec = dynamic_cast<NormVecDecorator*>( mTriangulation );
141  if ( dec )
142  {
143  QProgressDialog* progressDialog = nullptr;
144  if ( mShowProgressDialog ) //show a progress dialog because it can take a long time...
145  {
146  progressDialog = new QProgressDialog();
147  progressDialog->setLabelText( QObject::tr( "Estimating normal derivatives..." ) );
148  }
149  dec->estimateFirstDerivatives( progressDialog );
150  delete progressDialog;
151  ctInterpolator->setTriangulation( dec );
152  dec->setTriangleInterpolator( ctInterpolator );
153  mTriangleInterpolator = ctInterpolator;
154  }
155  }
156  else //linear
157  {
158  mTriangleInterpolator = new LinTriangleInterpolator( theDualEdgeTriangulation );
159  }
160  mIsInitialized = true;
161 
162  //debug
163  if ( mExportTriangulationToFile )
164  {
165  theDualEdgeTriangulation->saveAsShapefile( mTriangulationFilePath );
166  }
167 }
168 
169 int QgsTINInterpolator::insertData( QgsFeature* f, bool zCoord, int attr, InputType type )
170 {
171  if ( !f )
172  {
173  return 1;
174  }
175 
176  QgsGeometry g = f->geometry();
177  {
178  if ( g.isEmpty() )
179  {
180  return 2;
181  }
182  }
183 
184  //check attribute value
185  double attributeValue = 0;
186  bool attributeConversionOk = false;
187  if ( !zCoord )
188  {
189  QVariant attributeVariant = f->attribute( attr );
190  if ( !attributeVariant.isValid() ) //attribute not found, something must be wrong (e.g. NULL value)
191  {
192  return 3;
193  }
194  attributeValue = attributeVariant.toDouble( &attributeConversionOk );
195  if ( !attributeConversionOk || qIsNaN( attributeValue ) ) //don't consider vertices with attributes like 'nan' for the interpolation
196  {
197  return 4;
198  }
199  }
200 
201  //parse WKB. It is ugly, but we cannot use the methods with QgsPoint because they don't contain z-values for 25D types
202  bool hasZValue = false;
203  double x, y, z;
204  QByteArray wkb( g.exportToWkb() );
205  QgsConstWkbPtr currentWkbPtr( wkb );
206  currentWkbPtr.readHeader();
207  //maybe a structure or break line
208  Line3D* line = nullptr;
209 
210  QgsWkbTypes::Type wkbType = g.wkbType();
211  switch ( wkbType )
212  {
214  hasZValue = true;
215  FALLTHROUGH;
216  case QgsWkbTypes::Point:
217  {
218  currentWkbPtr >> x >> y;
219  if ( zCoord && hasZValue )
220  {
221  currentWkbPtr >> z;
222  }
223  else
224  {
225  z = attributeValue;
226  }
227  Point3D* thePoint = new Point3D( x, y, z );
228  if ( mTriangulation->addPoint( thePoint ) == -100 )
229  {
230  return -1;
231  }
232  break;
233  }
235  hasZValue = true;
236  FALLTHROUGH;
238  {
239  int nPoints;
240  currentWkbPtr >> nPoints;
241  for ( int index = 0; index < nPoints; ++index )
242  {
243  currentWkbPtr.readHeader();
244  currentWkbPtr >> x >> y;
245  if ( hasZValue ) //skip z-coordinate for 25D geometries
246  {
247  currentWkbPtr >> z;
248  }
249  else
250  {
251  z = attributeValue;
252  }
253  }
254  break;
255  }
257  hasZValue = true;
258  FALLTHROUGH;
260  {
261  if ( type != POINTS )
262  {
263  line = new Line3D();
264  }
265  int nPoints;
266  currentWkbPtr >> nPoints;
267  for ( int index = 0; index < nPoints; ++index )
268  {
269  currentWkbPtr >> x >> y;
270  if ( zCoord && hasZValue ) //skip z-coordinate for 25D geometries
271  {
272  currentWkbPtr >> z;
273  }
274  else
275  {
276  z = attributeValue;
277  }
278 
279  if ( type == POINTS )
280  {
281  //todo: handle error code -100
282  mTriangulation->addPoint( new Point3D( x, y, z ) );
283  }
284  else
285  {
286  line->insertPoint( new Point3D( x, y, z ) );
287  }
288  }
289 
290  if ( type != POINTS )
291  {
292  mTriangulation->addLine( line, type == BREAK_LINES );
293  }
294  break;
295  }
297  hasZValue = true;
298  FALLTHROUGH;
300  {
301  int nLines;
302  currentWkbPtr >> nLines;
303  for ( int index = 0; index < nLines; ++index )
304  {
305  if ( type != POINTS )
306  {
307  line = new Line3D();
308  }
309  int nPoints;
310  currentWkbPtr >> nPoints;
311  for ( int index2 = 0; index2 < nPoints; ++index2 )
312  {
313  currentWkbPtr >> x >> y;
314  if ( hasZValue ) //skip z-coordinate for 25D geometries
315  {
316  currentWkbPtr >> z;
317  }
318  else
319  {
320  z = attributeValue;
321  }
322 
323  if ( type == POINTS )
324  {
325  //todo: handle error code -100
326  mTriangulation->addPoint( new Point3D( x, y, z ) );
327  }
328  else
329  {
330  line->insertPoint( new Point3D( x, y, z ) );
331  }
332  }
333  if ( type != POINTS )
334  {
335  mTriangulation->addLine( line, type == BREAK_LINES );
336  }
337  }
338  break;
339  }
341  hasZValue = true;
342  FALLTHROUGH;
344  {
345  int nRings;
346  currentWkbPtr >> nRings;
347  for ( int index = 0; index < nRings; ++index )
348  {
349  if ( type != POINTS )
350  {
351  line = new Line3D();
352  }
353 
354  int nPoints;
355  currentWkbPtr >> nPoints;
356  for ( int index2 = 0; index2 < nPoints; ++index2 )
357  {
358  currentWkbPtr >> x >> y;
359  if ( hasZValue ) //skip z-coordinate for 25D geometries
360  {
361  currentWkbPtr >> z;
362  }
363  else
364  {
365  z = attributeValue;
366  }
367  if ( type == POINTS )
368  {
369  //todo: handle error code -100
370  mTriangulation->addPoint( new Point3D( x, y, z ) );
371  }
372  else
373  {
374  line->insertPoint( new Point3D( x, y, z ) );
375  }
376  }
377 
378  if ( type != POINTS )
379  {
380  mTriangulation->addLine( line, type == BREAK_LINES );
381  }
382  }
383  break;
384  }
385 
387  hasZValue = true;
388  FALLTHROUGH;
390  {
391  int nPolys;
392  currentWkbPtr >> nPolys;
393  for ( int index = 0; index < nPolys; ++index )
394  {
395  currentWkbPtr.readHeader();
396  int nRings;
397  currentWkbPtr >> nRings;
398  for ( int index2 = 0; index2 < nRings; ++index2 )
399  {
400  if ( type != POINTS )
401  {
402  line = new Line3D();
403  }
404  int nPoints;
405  currentWkbPtr >> nPoints;
406  for ( int index3 = 0; index3 < nPoints; ++index3 )
407  {
408  currentWkbPtr >> x >> y;
409  if ( hasZValue ) //skip z-coordinate for 25D geometries
410  {
411  currentWkbPtr >> z;
412  }
413  else
414  {
415  z = attributeValue;
416  }
417  if ( type == POINTS )
418  {
419  //todo: handle error code -100
420  mTriangulation->addPoint( new Point3D( x, y, z ) );
421  }
422  else
423  {
424  line->insertPoint( new Point3D( x, y, z ) );
425  }
426  }
427  if ( type != POINTS )
428  {
429  mTriangulation->addLine( line, type == BREAK_LINES );
430  }
431  }
432  }
433  break;
434  }
435  default:
436  //should not happen...
437  break;
438  }
439 
440  return 0;
441 }
442 
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.
static unsigned index
QList< LayerData > mLayerData
virtual int addPoint(Point3D *p)=0
Adds a point to the triangulation Ownership is transferred to this class.
Interface class for interpolations.
QgsVectorLayer * vectorLayer
virtual void addTriangulation(Triangulation *t)
Adds an association to a triangulation.
Definition: TriDecorator.h:78
LinTriangleInterpolator is a class which interpolates linearly on a triangulation.
QgsWkbTypes::Type wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
double getZ() const
Returns the z-coordinate of the point.
Definition: Point3D.h:94
virtual void setTriangulation(NormVecDecorator *tin)
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:78
QgsTINInterpolator(const QList< LayerData > &inputData, TIN_INTERPOLATION interpolation=Linear, bool showProgressDialog=false)
This class represents a line.
Definition: Line3D.h:24
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:135
DualEdgeTriangulation is an implementation of a triangulation class based on the dual edge data struc...
InputType
Describes the type of input data.
Point3D is a class to represent a three dimensional point.
Definition: Point3D.h:24
#define FALLTHROUGH
Definition: qgis.h:375
QByteArray exportToWkb() const
Export the geometry to WKB.
long featureCount(const QString &legendKey) const
Number of features rendered with specified legend key.
QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const
Query the layer for features specified in request.
This class wraps a request for features to a vector layer (or directly its vector data provider)...
QList< int > QgsAttributeList
bool isEmpty() const
Returns true if the geometry is empty (ie, contains no underlying geometry accessible via geometry)...
virtual void addLine(Line3D *line, bool breakline)=0
Adds a line (e.g.
QgsGeometry geometry() const
Returns the geometry associated with this feature.
Definition: qgsfeature.cpp:113
void setTriangleInterpolator(TriangleInterpolator *inter) override
Sets an interpolator.
bool estimateFirstDerivatives(QProgressDialog *d=nullptr)
This method adds the functionality of estimating normals at the data points. Return true in the case ...
void insertPoint(Point3D *p)
Inserts a node behind the current position and sets the current position to this new node...
Definition: Line3D.cc:49
A layer together with the information about interpolation attribute / z-coordinate interpolation and ...
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.
virtual bool calcPoint(double x, double y, Point3D *result)=0
Performs a linear interpolation in a triangle and assigns the x-,y- and z-coordinates to point...
bool nextFeature(QgsFeature &f)
QVariant attribute(const QString &name) const
Lookup attribute value from attribute name.
Definition: qgsfeature.cpp:277
virtual bool saveAsShapefile(const QString &fileName) const override
Saves the triangulation as a (line) shapefile.