QGIS API Documentation  3.8.0-Zanzibar (11aff65)
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 "qgspoint.h"
25 #include "qgsfeature.h"
26 #include "qgsgeometry.h"
27 #include "qgsvectorlayer.h"
28 #include "qgswkbptr.h"
29 #include "qgsfeedback.h"
30 #include "qgscurve.h"
31 #include "qgsmulticurve.h"
32 #include "qgscurvepolygon.h"
33 #include "qgsmultisurface.h"
34 
35 QgsTinInterpolator::QgsTinInterpolator( const QList<LayerData> &inputData, TinInterpolation interpolation, QgsFeedback *feedback )
36  : QgsInterpolator( inputData )
37  , mIsInitialized( false )
38  , mFeedback( feedback )
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, QgsFeedback * )
50 {
51  if ( !mIsInitialized )
52  {
53  initialize();
54  }
55 
56  if ( !mTriangleInterpolator )
57  {
58  return 1;
59  }
60 
61  QgsPoint r( 0, 0, 0 );
62  if ( !mTriangleInterpolator->calcPoint( x, y, r ) )
63  {
64  return 2;
65  }
66  result = r.z();
67  return 0;
68 }
69 
71 {
73 }
74 
76 {
77  mTriangulationSink = sink;
78 }
79 
80 void QgsTinInterpolator::initialize()
81 {
82  DualEdgeTriangulation *dualEdgeTriangulation = new DualEdgeTriangulation( 100000, nullptr );
83  if ( mInterpolation == CloughTocher )
84  {
86  dec->addTriangulation( dualEdgeTriangulation );
87  mTriangulation = dec;
88  }
89  else
90  {
91  mTriangulation = dualEdgeTriangulation;
92  }
93 
94  //get number of features if we use a progress bar
95  int nFeatures = 0;
96  int nProcessedFeatures = 0;
97  if ( mFeedback )
98  {
99  for ( const LayerData &layer : qgis::as_const( mLayerData ) )
100  {
101  if ( layer.source )
102  {
103  nFeatures += layer.source->featureCount();
104  }
105  }
106  }
107 
108  QgsFeature f;
109  for ( const LayerData &layer : qgis::as_const( mLayerData ) )
110  {
111  if ( layer.source )
112  {
113  QgsAttributeList attList;
114  switch ( layer.valueSource )
115  {
117  attList.push_back( layer.interpolationAttribute );
118  break;
119 
122  break;
123  }
124 
125  QgsFeatureIterator fit = layer.source->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( attList ) );
126 
127  while ( fit.nextFeature( f ) )
128  {
129  if ( mFeedback )
130  {
131  if ( mFeedback->isCanceled() )
132  {
133  break;
134  }
135  if ( nFeatures > 0 )
136  mFeedback->setProgress( 100.0 * static_cast< double >( nProcessedFeatures ) / nFeatures );
137  }
138  insertData( f, layer.valueSource, layer.interpolationAttribute, layer.sourceType );
139  ++nProcessedFeatures;
140  }
141  }
142  }
143 
144  if ( mInterpolation == CloughTocher )
145  {
146  CloughTocherInterpolator *ctInterpolator = new CloughTocherInterpolator();
147  NormVecDecorator *dec = dynamic_cast<NormVecDecorator *>( mTriangulation );
148  if ( dec )
149  {
150  dec->estimateFirstDerivatives( mFeedback );
151  ctInterpolator->setTriangulation( dec );
152  dec->setTriangleInterpolator( ctInterpolator );
153  mTriangleInterpolator = ctInterpolator;
154  }
155  }
156  else //linear
157  {
158  mTriangleInterpolator = new LinTriangleInterpolator( dualEdgeTriangulation );
159  }
160  mIsInitialized = true;
161 
162  //debug
163  if ( mTriangulationSink )
164  {
165  dualEdgeTriangulation->saveTriangulation( mTriangulationSink, mFeedback );
166  }
167 }
168 
169 int QgsTinInterpolator::insertData( const QgsFeature &f, QgsInterpolator::ValueSource source, int attr, SourceType type )
170 {
171  QgsGeometry g = f.geometry();
172  if ( g.isNull() || g.isEmpty() )
173  {
174  return 2;
175  }
176 
177  //check attribute value
178  double attributeValue = 0;
179  bool attributeConversionOk = false;
180  switch ( source )
181  {
182  case ValueAttribute:
183  {
184  QVariant attributeVariant = f.attribute( attr );
185  if ( !attributeVariant.isValid() ) //attribute not found, something must be wrong (e.g. NULL value)
186  {
187  return 3;
188  }
189  attributeValue = attributeVariant.toDouble( &attributeConversionOk );
190  if ( !attributeConversionOk || std::isnan( attributeValue ) ) //don't consider vertices with attributes like 'nan' for the interpolation
191  {
192  return 4;
193  }
194  break;
195  }
196 
197  case ValueM:
198  if ( !g.constGet()->isMeasure() )
199  return 3;
200  else
201  break;
202 
203  case ValueZ:
204  if ( !g.constGet()->is3D() )
205  return 3;
206  else
207  break;
208  }
209 
210 
211  switch ( type )
212  {
213  case SourcePoints:
214  {
215  if ( addPointsFromGeometry( g, source, attributeValue ) != 0 )
216  return -1;
217  break;
218  }
219 
220  case SourceBreakLines:
222  {
223  switch ( QgsWkbTypes::geometryType( g.wkbType() ) )
224  {
226  {
227  if ( addPointsFromGeometry( g, source, attributeValue ) != 0 )
228  return -1;
229  break;
230  }
231 
234  {
235  // need to extract all rings from input geometry
236  std::vector<const QgsCurve *> curves;
238  {
239  std::vector< const QgsCurvePolygon * > polygons;
240  if ( g.isMultipart() )
241  {
242  const QgsMultiSurface *ms = qgsgeometry_cast< const QgsMultiSurface * >( g.constGet() );
243  for ( int i = 0; i < ms->numGeometries(); ++i )
244  {
245  polygons.emplace_back( qgsgeometry_cast< const QgsCurvePolygon * >( ms->geometryN( i ) ) );
246  }
247  }
248  else
249  {
250  polygons.emplace_back( qgsgeometry_cast< const QgsCurvePolygon * >( g.constGet() ) );
251  }
252 
253  for ( const QgsCurvePolygon *polygon : polygons )
254  {
255  if ( !polygon )
256  continue;
257 
258  if ( polygon->exteriorRing() )
259  curves.emplace_back( polygon->exteriorRing() );
260 
261  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
262  {
263  curves.emplace_back( polygon->interiorRing( i ) );
264  }
265  }
266  }
267  else
268  {
269  if ( g.isMultipart() )
270  {
271  const QgsMultiCurve *mc = qgsgeometry_cast< const QgsMultiCurve * >( g.constGet() );
272  for ( int i = 0; i < mc->numGeometries(); ++i )
273  {
274  curves.emplace_back( qgsgeometry_cast< const QgsCurve * >( mc->geometryN( i ) ) );
275  }
276  }
277  else
278  {
279  curves.emplace_back( qgsgeometry_cast< const QgsCurve * >( g.constGet() ) );
280  }
281  }
282 
283  for ( const QgsCurve *curve : curves )
284  {
285  if ( !curve )
286  continue;
287 
288  QVector< QgsPoint > linePoints;
289  for ( auto point = g.vertices_begin(); point != g.vertices_end(); ++point )
290  {
291  QgsPoint p = *point;
292  double z = 0;
293  switch ( source )
294  {
295  case ValueAttribute:
296  z = attributeValue;
297  break;
298 
299  case ValueZ:
300  z = p.z();
301  break;
302 
303  case ValueM:
304  z = p.m();
305  break;
306  }
307 
308  linePoints.append( QgsPoint( p.x(), p.y(), z ) );
309  }
310  mTriangulation->addLine( linePoints, type );
311  }
312  break;
313  }
316  break;
317  }
318  break;
319  }
320  }
321 
322  return 0;
323 }
324 
325 
326 int QgsTinInterpolator::addPointsFromGeometry( const QgsGeometry &g, ValueSource source, double attributeValue )
327 {
328  // loop through all vertices and add to triangulation
329  for ( auto point = g.vertices_begin(); point != g.vertices_end(); ++point )
330  {
331  QgsPoint p = *point;
332  double z = 0;
333  switch ( source )
334  {
335  case ValueAttribute:
336  z = attributeValue;
337  break;
338 
339  case ValueZ:
340  z = p.z();
341  break;
342 
343  case ValueM:
344  z = p.m();
345  break;
346  }
347  if ( mTriangulation->addPoint( QgsPoint( p.x(), p.y(), z ) ) == -100 )
348  {
349  return -1;
350  }
351  }
352  return 0;
353 }
bool isMeasure() const
Returns true if the geometry contains m values.
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
Information about the input vector layers and the attributes (or z-values) that are used for interpol...
double y
Definition: qgspoint.h:42
Interface class for interpolations.
bool isMultipart() const
Returns true if WKB of the geometry is of WKBMulti* type.
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.)
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
An interface for objects which accept features via addFeature(s) methods.
Curve polygon geometry type.
Clough-Tocher interpolation.
Container of fields for a vector layer.
Definition: qgsfields.h:42
virtual void setTriangulation(NormVecDecorator *tin)
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:111
SourceType
Describes the type of input data.
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:55
Multi surface geometry collection.
DualEdgeTriangulation is an implementation of a triangulation class based on the dual edge data struc...
Base class for feedback objects to be used for cancellation of something running in a worker thread...
Definition: qgsfeedback.h:44
ValueSource
Source for interpolated values from features.
Take value from feature&#39;s attribute.
static QgsFields triangulationFields()
Returns the fields output by features when saving the triangulation.
static GeometryType geometryType(Type type)
Returns the geometry type for a WKB type, e.g., both MultiPolygon and CurvePolygon would have a Polyg...
Definition: qgswkbtypes.h:666
bool saveTriangulation(QgsFeatureSink *sink, QgsFeedback *feedback=nullptr) const override
Saves the triangulation features to a feature sink.
This class wraps a request for features to a vector layer (or directly its vector data provider)...
T qgsgeometry_cast(const QgsAbstractGeometry *geom)
bool isEmpty() const
Returns true if the geometry is empty (eg a linestring with no vertices, or a collection with no geom...
Multi curve geometry collection.
Definition: qgsmulticurve.h:29
Abstract base class for curved geometry type.
Definition: qgscurve.h:35
int interpolatePoint(double x, double y, double &result, QgsFeedback *feedback) override
Calculates interpolation value for map coordinates x, y.
void setTriangleInterpolator(TriangleInterpolator *inter) override
Sets an interpolator.
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
Use feature&#39;s geometry Z values for interpolation.
virtual bool calcPoint(double x, double y, QgsPoint &result)=0
Performs a linear interpolation in a triangle and assigns the x-,y- and z-coordinates to point...
int numGeometries() const
Returns the number of geometries within the collection.
virtual int addPoint(const QgsPoint &point)=0
Adds a point to the triangulation.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
A source together with the information about interpolation attribute / z-coordinate interpolation and...
TinInterpolation
Indicates the type of interpolation to be performed.
QgsTinInterpolator(const QList< QgsInterpolator::LayerData > &inputData, TinInterpolation interpolation=Linear, QgsFeedback *feedback=nullptr)
Constructor for QgsTinInterpolator.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
void setTriangulationSink(QgsFeatureSink *sink)
Sets the optional sink for saving the triangulation features.
static QgsFields triangulationFields()
Returns the fields output by features when calling saveTriangulation().
Use feature&#39;s geometry M values for interpolation.
This is an implementation of a Clough-Tocher interpolator based on a triangular tessellation.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
double z
Definition: qgspoint.h:43
QgsGeometry geometry
Definition: qgsfeature.h:67
virtual void addLine(const QVector< QgsPoint > &points, QgsInterpolator::SourceType lineType)=0
Adds a line (e.g.
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 ...
QgsAbstractGeometry::vertex_iterator vertices_end() const
Returns STL-style iterator pointing to the imaginary vertex after the last vertex of the geometry...
QVariant attribute(const QString &name) const
Lookup attribute value from attribute name.
Definition: qgsfeature.cpp:262
QgsAbstractGeometry::vertex_iterator vertices_begin() const
Returns STL-style iterator pointing to the first vertex of the geometry.
bool is3D() const
Returns true if the geometry is 3D and contains a z-value.
double m
Definition: qgspoint.h:44
double x
Definition: qgspoint.h:41