QGIS API Documentation  2.11.0-Master
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( 0 )
33  , mTriangleInterpolator( 0 )
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, 0 );
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  QList<LayerData>::iterator layerDataIt = mLayerData.begin();
88  for ( ; layerDataIt != mLayerData.end(); ++layerDataIt )
89  {
90  if ( layerDataIt->vectorLayer )
91  {
92  nFeatures += layerDataIt->vectorLayer->featureCount();
93  }
94  }
95  }
96 
97  QProgressDialog* theProgressDialog = 0;
98  if ( mShowProgressDialog )
99  {
100  theProgressDialog = new QProgressDialog( QObject::tr( "Building triangulation..." ), QObject::tr( "Abort" ), 0, nFeatures, 0 );
101  theProgressDialog->setWindowModality( Qt::WindowModal );
102  }
103 
104 
105  QgsFeature f;
106  QList<LayerData>::iterator layerDataIt = mLayerData.begin();
107  for ( ; layerDataIt != mLayerData.end(); ++layerDataIt )
108  {
109  if ( layerDataIt->vectorLayer )
110  {
111  QgsAttributeList attList;
112  if ( !layerDataIt->zCoordInterpolation )
113  {
114  attList.push_back( layerDataIt->interpolationAttribute );
115  }
116 
117  QgsFeatureIterator fit = layerDataIt->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, layerDataIt->zCoordInterpolation, layerDataIt->interpolationAttribute, layerDataIt->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 = 0;
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  const QgsGeometry* g = f->constGeometry();
177  {
178  if ( !g )
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  QgsConstWkbPtr currentWkbPtr( g->asWkb() + 1 + sizeof( int ) );
205  //maybe a structure or break line
206  Line3D* line = 0;
207 
208  QGis::WkbType wkbType = g->wkbType();
209  switch ( wkbType )
210  {
211  case QGis::WKBPoint25D:
212  hasZValue = true;
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  case QGis::WKBMultiPoint:
234  {
235  int nPoints;
236  currentWkbPtr >> nPoints;
237  for ( int index = 0; index < nPoints; ++index )
238  {
239  currentWkbPtr += 1 + sizeof( int );
240  currentWkbPtr >> x >> y;
241  if ( hasZValue ) //skip z-coordinate for 25D geometries
242  {
243  currentWkbPtr >> z;
244  }
245  else
246  {
247  z = attributeValue;
248  }
249  }
250  break;
251  }
253  hasZValue = true;
254  case QGis::WKBLineString:
255  {
256  if ( type != POINTS )
257  {
258  line = new Line3D();
259  }
260  int nPoints;
261  currentWkbPtr >> nPoints;
262  for ( int index = 0; index < nPoints; ++index )
263  {
264  currentWkbPtr >> x >> y;
265  if ( zCoord && hasZValue ) //skip z-coordinate for 25D geometries
266  {
267  currentWkbPtr >> z;
268  }
269  else
270  {
271  z = attributeValue;
272  }
273 
274  if ( type == POINTS )
275  {
276  //todo: handle error code -100
277  mTriangulation->addPoint( new Point3D( x, y, z ) );
278  }
279  else
280  {
281  line->insertPoint( new Point3D( x, y, z ) );
282  }
283  }
284 
285  if ( type != POINTS )
286  {
287  mTriangulation->addLine( line, type == BREAK_LINES );
288  }
289  break;
290  }
292  hasZValue = true;
294  {
295  int nLines;
296  currentWkbPtr >> nLines;
297  for ( int index = 0; index < nLines; ++index )
298  {
299  if ( type != POINTS )
300  {
301  line = new Line3D();
302  }
303  int nPoints;
304  currentWkbPtr >> nPoints;
305  for ( int index2 = 0; index2 < nPoints; ++index2 )
306  {
307  currentWkbPtr >> x >> y;
308  if ( hasZValue ) //skip z-coordinate for 25D geometries
309  {
310  currentWkbPtr >> z;
311  }
312  else
313  {
314  z = attributeValue;
315  }
316 
317  if ( type == POINTS )
318  {
319  //todo: handle error code -100
320  mTriangulation->addPoint( new Point3D( x, y, z ) );
321  }
322  else
323  {
324  line->insertPoint( new Point3D( x, y, z ) );
325  }
326  }
327  if ( type != POINTS )
328  {
329  mTriangulation->addLine( line, type == BREAK_LINES );
330  }
331  }
332  break;
333  }
334  case QGis::WKBPolygon25D:
335  hasZValue = true;
336  case QGis::WKBPolygon:
337  {
338  int nRings;
339  currentWkbPtr >> nRings;
340  for ( int index = 0; index < nRings; ++index )
341  {
342  if ( type != POINTS )
343  {
344  line = new Line3D();
345  }
346 
347  int nPoints;
348  currentWkbPtr >> nPoints;
349  for ( int index2 = 0; index2 < nPoints; ++index2 )
350  {
351  currentWkbPtr >> x >> y;
352  if ( hasZValue ) //skip z-coordinate for 25D geometries
353  {
354  currentWkbPtr >> z;
355  }
356  else
357  {
358  z = attributeValue;
359  }
360  if ( type == POINTS )
361  {
362  //todo: handle error code -100
363  mTriangulation->addPoint( new Point3D( x, y, z ) );
364  }
365  else
366  {
367  line->insertPoint( new Point3D( x, y, z ) );
368  }
369  }
370 
371  if ( type != POINTS )
372  {
373  mTriangulation->addLine( line, type == BREAK_LINES );
374  }
375  }
376  break;
377  }
378 
380  hasZValue = true;
382  {
383  int nPolys;
384  currentWkbPtr >> nPolys;
385  for ( int index = 0; index < nPolys; ++index )
386  {
387  currentWkbPtr += 1 + sizeof( int );
388  int nRings;
389  currentWkbPtr >> nRings;
390  for ( int index2 = 0; index2 < nRings; ++index2 )
391  {
392  if ( type != POINTS )
393  {
394  line = new Line3D();
395  }
396  int nPoints;
397  currentWkbPtr >> nPoints;
398  for ( int index3 = 0; index3 < nPoints; ++index3 )
399  {
400  currentWkbPtr >> x >> y;
401  if ( hasZValue ) //skip z-coordinate for 25D geometries
402  {
403  currentWkbPtr >> z;
404  }
405  else
406  {
407  z = attributeValue;
408  }
409  if ( type == POINTS )
410  {
411  //todo: handle error code -100
412  mTriangulation->addPoint( new Point3D( x, y, z ) );
413  }
414  else
415  {
416  line->insertPoint( new Point3D( x, y, z ) );
417  }
418  }
419  if ( type != POINTS )
420  {
421  mTriangulation->addLine( line, type == BREAK_LINES );
422  }
423  }
424  }
425  break;
426  }
427  default:
428  //should not happen...
429  break;
430  }
431 
432  return 0;
433 }
434 
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.
virtual void addTriangulation(Triangulation *t)
Adds an association to a triangulation.
Definition: TriDecorator.h:78
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)
void setLabelText(const QString &text)
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:75
QgsTINInterpolator(const QList< LayerData > &inputData, TIN_INTERPOLATION interpolation=Linear, bool showProgressDialog=false)
WkbType
Used for symbology operations.
Definition: qgis.h:53
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:162
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...
void setValue(int progress)
Point3D is a class to represent a three dimensional point.
Definition: Point3D.h:23
bool estimateFirstDerivatives(QProgressDialog *d=0)
This method adds the functionality of estimating normals at the data points.
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.
void insertPoint(Point3D *p)
Inserts a node behind the current position and sets the current position to this new node...
QVariant attribute(const QString &name) const
Lookup attribute value from attribute name.
Definition: qgsfeature.cpp:236
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:68
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.