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