QGIS API Documentation  2.2.0-Valmiera
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
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 <QProgressDialog>
28 
29 QgsTINInterpolator::QgsTINInterpolator( const QList<LayerData>& inputData, TIN_INTERPOLATION interpolation, bool showProgressDialog )
30  : QgsInterpolator( inputData )
31  , mTriangulation( 0 )
32  , mTriangleInterpolator( 0 )
33  , mIsInitialized( false )
34  , mShowProgressDialog( showProgressDialog )
35  , mExportTriangulationToFile( false )
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  Point3D r;
59  if ( !mTriangleInterpolator->calcPoint( x, y, &r ) )
60  {
61  return 2;
62  }
63  result = r.getZ();
64  return 0;
65 }
66 
68 {
69  DualEdgeTriangulation* theDualEdgeTriangulation = new DualEdgeTriangulation( 100000, 0 );
71  {
73  dec->addTriangulation( theDualEdgeTriangulation );
74  mTriangulation = dec;
75  }
76  else
77  {
78  mTriangulation = theDualEdgeTriangulation;
79  }
80 
81  //get number of features if we use a progress bar
82  int nFeatures = 0;
83  int nProcessedFeatures = 0;
84  if ( mShowProgressDialog )
85  {
86  QList<LayerData>::iterator layerDataIt = mLayerData.begin();
87  for ( ; layerDataIt != mLayerData.end(); ++layerDataIt )
88  {
89  if ( layerDataIt->vectorLayer )
90  {
91  nFeatures += layerDataIt->vectorLayer->featureCount();
92  }
93  }
94  }
95 
96  QProgressDialog* theProgressDialog = 0;
97  if ( mShowProgressDialog )
98  {
99  theProgressDialog = new QProgressDialog( QObject::tr( "Building triangulation..." ), QObject::tr( "Abort" ), 0, nFeatures, 0 );
100  theProgressDialog->setWindowModality( Qt::WindowModal );
101  }
102 
103 
104  QgsFeature f;
105  QList<LayerData>::iterator layerDataIt = mLayerData.begin();
106  for ( ; layerDataIt != mLayerData.end(); ++layerDataIt )
107  {
108  if ( layerDataIt->vectorLayer )
109  {
110  QgsAttributeList attList;
111  if ( !layerDataIt->zCoordInterpolation )
112  {
113  attList.push_back( layerDataIt->interpolationAttribute );
114  }
115 
116  QgsFeatureIterator fit = layerDataIt->vectorLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( attList ) );
117 
118  while ( fit.nextFeature( f ) )
119  {
120  if ( mShowProgressDialog )
121  {
122  if ( theProgressDialog->wasCanceled() )
123  {
124  break;
125  }
126  theProgressDialog->setValue( nProcessedFeatures );
127  }
128  insertData( &f, layerDataIt->zCoordInterpolation, layerDataIt->interpolationAttribute, layerDataIt->mInputType );
129  ++nProcessedFeatures;
130  }
131  }
132  }
133 
134  delete theProgressDialog;
135 
136  if ( mInterpolation == CloughTocher )
137  {
138  CloughTocherInterpolator* ctInterpolator = new CloughTocherInterpolator();
139  NormVecDecorator* dec = dynamic_cast<NormVecDecorator*>( mTriangulation );
140  if ( dec )
141  {
142  QProgressDialog* progressDialog = 0;
143  if ( mShowProgressDialog ) //show a progress dialog because it can take a long time...
144  {
145  progressDialog = new QProgressDialog();
146  progressDialog->setLabelText( QObject::tr( "Estimating normal derivatives..." ) );
147  }
148  dec->estimateFirstDerivatives( progressDialog );
149  delete progressDialog;
150  ctInterpolator->setTriangulation( dec );
151  dec->setTriangleInterpolator( ctInterpolator );
152  mTriangleInterpolator = ctInterpolator;
153  }
154  }
155  else //linear
156  {
157  mTriangleInterpolator = new LinTriangleInterpolator( theDualEdgeTriangulation );
158  }
159  mIsInitialized = true;
160 
161  //debug
163  {
164  theDualEdgeTriangulation->saveAsShapefile( mTriangulationFilePath );
165  }
166 }
167 
168 int QgsTINInterpolator::insertData( QgsFeature* f, bool zCoord, int attr, InputType type )
169 {
170  if ( !f )
171  {
172  return 1;
173  }
174 
175  QgsGeometry* g = f->geometry();
176  {
177  if ( !g )
178  {
179  return 2;
180  }
181  }
182 
183  //check attribute value
184  double attributeValue = 0;
185  bool attributeConversionOk = false;
186  if ( !zCoord )
187  {
188  QVariant attributeVariant = f->attribute( attr );
189  if ( !attributeVariant.isValid() ) //attribute not found, something must be wrong (e.g. NULL value)
190  {
191  return 3;
192  }
193  attributeValue = attributeVariant.toDouble( &attributeConversionOk );
194  if ( !attributeConversionOk || qIsNaN( attributeValue ) ) //don't consider vertices with attributes like 'nan' for the interpolation
195  {
196  return 4;
197  }
198  }
199 
200  //parse WKB. It is ugly, but we cannot use the methods with QgsPoint because they don't contain z-values for 25D types
201  bool hasZValue = false;
202  double x, y, z;
203  QgsConstWkbPtr currentWkbPtr( g->asWkb() + 1 + sizeof( int ) );
204  //maybe a structure or break line
205  Line3D* line = 0;
206 
207  QGis::WkbType wkbType = g->wkbType();
208  switch ( wkbType )
209  {
210  case QGis::WKBPoint25D:
211  hasZValue = true;
212  case QGis::WKBPoint:
213  {
214  currentWkbPtr >> x >> y;
215  if ( zCoord && hasZValue )
216  {
217  currentWkbPtr >> z;
218  }
219  else
220  {
221  z = attributeValue;
222  }
223  Point3D* thePoint = new Point3D( x, y, z );
224  if ( mTriangulation->addPoint( thePoint ) == -100 )
225  {
226  return -1;
227  }
228  break;
229  }
231  hasZValue = true;
232  case QGis::WKBMultiPoint:
233  {
234  int nPoints;
235  currentWkbPtr >> nPoints;
236  for ( int index = 0; index < nPoints; ++index )
237  {
238  currentWkbPtr += 1 + sizeof( int );
239  currentWkbPtr >> x >> y;
240  if ( hasZValue ) //skip z-coordinate for 25D geometries
241  {
242  currentWkbPtr >> z;
243  }
244  else
245  {
246  z = attributeValue;
247  }
248  }
249  break;
250  }
252  hasZValue = true;
253  case QGis::WKBLineString:
254  {
255  if ( type != POINTS )
256  {
257  line = new Line3D();
258  }
259  int nPoints;
260  currentWkbPtr >> nPoints;
261  for ( int index = 0; index < nPoints; ++index )
262  {
263  currentWkbPtr >> x >> y;
264  if ( zCoord && hasZValue ) //skip z-coordinate for 25D geometries
265  {
266  currentWkbPtr >> z;
267  }
268  else
269  {
270  z = attributeValue;
271  }
272 
273  if ( type == POINTS )
274  {
275  //todo: handle error code -100
276  mTriangulation->addPoint( new Point3D( x, y, z ) );
277  }
278  else
279  {
280  line->insertPoint( new Point3D( x, y, z ) );
281  }
282  }
283 
284  if ( type != POINTS )
285  {
286  mTriangulation->addLine( line, type == BREAK_LINES );
287  }
288  break;
289  }
291  hasZValue = true;
293  {
294  int nLines;
295  currentWkbPtr >> nLines;
296  for ( int index = 0; index < nLines; ++index )
297  {
298  if ( type != POINTS )
299  {
300  line = new Line3D();
301  }
302  int nPoints;
303  currentWkbPtr >> nPoints;
304  for ( int index2 = 0; index2 < nPoints; ++index2 )
305  {
306  currentWkbPtr >> x >> y;
307  if ( hasZValue ) //skip z-coordinate for 25D geometries
308  {
309  currentWkbPtr >> z;
310  }
311  else
312  {
313  z = attributeValue;
314  }
315 
316  if ( type == POINTS )
317  {
318  //todo: handle error code -100
319  mTriangulation->addPoint( new Point3D( x, y, z ) );
320  }
321  else
322  {
323  line->insertPoint( new Point3D( x, y, z ) );
324  }
325  }
326  if ( type != POINTS )
327  {
328  mTriangulation->addLine( line, type == BREAK_LINES );
329  }
330  }
331  break;
332  }
333  case QGis::WKBPolygon25D:
334  hasZValue = true;
335  case QGis::WKBPolygon:
336  {
337  int nRings;
338  currentWkbPtr >> nRings;
339  for ( int index = 0; index < nRings; ++index )
340  {
341  if ( type != POINTS )
342  {
343  line = new Line3D();
344  }
345 
346  int nPoints;
347  currentWkbPtr >> nPoints;
348  for ( int index2 = 0; index2 < nPoints; ++index2 )
349  {
350  currentWkbPtr >> x >> y;
351  if ( hasZValue ) //skip z-coordinate for 25D geometries
352  {
353  currentWkbPtr >> z;
354  }
355  else
356  {
357  z = attributeValue;
358  }
359  if ( type == POINTS )
360  {
361  //todo: handle error code -100
362  mTriangulation->addPoint( new Point3D( x, y, z ) );
363  }
364  else
365  {
366  line->insertPoint( new Point3D( x, y, z ) );
367  }
368  }
369 
370  if ( type != POINTS )
371  {
372  mTriangulation->addLine( line, type == BREAK_LINES );
373  }
374  }
375  break;
376  }
377 
379  hasZValue = true;
381  {
382  int nPolys;
383  currentWkbPtr >> nPolys;
384  for ( int index = 0; index < nPolys; ++index )
385  {
386  currentWkbPtr += 1 + sizeof( int );
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 Point3D( x, y, z ) );
412  }
413  else
414  {
415  line->insertPoint( new Point3D( 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