QGIS API Documentation  3.14.0-Pi (9f7028fd23)
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  const QgsCoordinateReferenceSystem crs = !mLayerData.empty() ? mLayerData.at( 0 ).source->sourceCrs() : QgsCoordinateReferenceSystem();
109 
110  QgsFeature f;
111  for ( const LayerData &layer : qgis::as_const( mLayerData ) )
112  {
113  if ( layer.source )
114  {
115  QgsAttributeList attList;
116  switch ( layer.valueSource )
117  {
119  attList.push_back( layer.interpolationAttribute );
120  break;
121 
124  break;
125  }
126 
127  QgsFeatureIterator fit = layer.source->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( attList ).setDestinationCrs( crs, layer.transformContext ) );
128 
129  while ( fit.nextFeature( f ) )
130  {
131  if ( mFeedback )
132  {
133  if ( mFeedback->isCanceled() )
134  {
135  break;
136  }
137  if ( nFeatures > 0 )
138  mFeedback->setProgress( 100.0 * static_cast< double >( nProcessedFeatures ) / nFeatures );
139  }
140  insertData( f, layer.valueSource, layer.interpolationAttribute, layer.sourceType );
141  ++nProcessedFeatures;
142  }
143  }
144  }
145 
146  if ( mInterpolation == CloughTocher )
147  {
148  CloughTocherInterpolator *ctInterpolator = new CloughTocherInterpolator();
149  NormVecDecorator *dec = dynamic_cast<NormVecDecorator *>( mTriangulation );
150  if ( dec )
151  {
152  dec->estimateFirstDerivatives( mFeedback );
153  ctInterpolator->setTriangulation( dec );
154  dec->setTriangleInterpolator( ctInterpolator );
155  mTriangleInterpolator = ctInterpolator;
156  }
157  }
158  else //linear
159  {
160  mTriangleInterpolator = new LinTriangleInterpolator( dualEdgeTriangulation );
161  }
162  mIsInitialized = true;
163 
164  //debug
165  if ( mTriangulationSink )
166  {
167  dualEdgeTriangulation->saveTriangulation( mTriangulationSink, mFeedback );
168  }
169 }
170 
171 int QgsTinInterpolator::insertData( const QgsFeature &f, QgsInterpolator::ValueSource source, int attr, SourceType type )
172 {
173  QgsGeometry g = f.geometry();
174  if ( g.isNull() || g.isEmpty() )
175  {
176  return 2;
177  }
178 
179  //check attribute value
180  double attributeValue = 0;
181  bool attributeConversionOk = false;
182  switch ( source )
183  {
184  case ValueAttribute:
185  {
186  QVariant attributeVariant = f.attribute( attr );
187  if ( !attributeVariant.isValid() ) //attribute not found, something must be wrong (e.g. NULL value)
188  {
189  return 3;
190  }
191  attributeValue = attributeVariant.toDouble( &attributeConversionOk );
192  if ( !attributeConversionOk || std::isnan( attributeValue ) ) //don't consider vertices with attributes like 'nan' for the interpolation
193  {
194  return 4;
195  }
196  break;
197  }
198 
199  case ValueM:
200  if ( !g.constGet()->isMeasure() )
201  return 3;
202  else
203  break;
204 
205  case ValueZ:
206  if ( !g.constGet()->is3D() )
207  return 3;
208  else
209  break;
210  }
211 
212 
213  switch ( type )
214  {
215  case SourcePoints:
216  {
217  if ( addPointsFromGeometry( g, source, attributeValue ) != 0 )
218  return -1;
219  break;
220  }
221 
222  case SourceBreakLines:
224  {
225  switch ( QgsWkbTypes::geometryType( g.wkbType() ) )
226  {
228  {
229  if ( addPointsFromGeometry( g, source, attributeValue ) != 0 )
230  return -1;
231  break;
232  }
233 
236  {
237  // need to extract all rings from input geometry
238  std::vector<const QgsCurve *> curves;
240  {
241  std::vector< const QgsCurvePolygon * > polygons;
242  if ( g.isMultipart() )
243  {
244  const QgsMultiSurface *ms = qgsgeometry_cast< const QgsMultiSurface * >( g.constGet() );
245  for ( int i = 0; i < ms->numGeometries(); ++i )
246  {
247  polygons.emplace_back( qgsgeometry_cast< const QgsCurvePolygon * >( ms->geometryN( i ) ) );
248  }
249  }
250  else
251  {
252  polygons.emplace_back( qgsgeometry_cast< const QgsCurvePolygon * >( g.constGet() ) );
253  }
254 
255  for ( const QgsCurvePolygon *polygon : polygons )
256  {
257  if ( !polygon )
258  continue;
259 
260  if ( polygon->exteriorRing() )
261  curves.emplace_back( polygon->exteriorRing() );
262 
263  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
264  {
265  curves.emplace_back( polygon->interiorRing( i ) );
266  }
267  }
268  }
269  else
270  {
271  if ( g.isMultipart() )
272  {
273  const QgsMultiCurve *mc = qgsgeometry_cast< const QgsMultiCurve * >( g.constGet() );
274  for ( int i = 0; i < mc->numGeometries(); ++i )
275  {
276  curves.emplace_back( qgsgeometry_cast< const QgsCurve * >( mc->geometryN( i ) ) );
277  }
278  }
279  else
280  {
281  curves.emplace_back( qgsgeometry_cast< const QgsCurve * >( g.constGet() ) );
282  }
283  }
284 
285  for ( const QgsCurve *curve : curves )
286  {
287  if ( !curve )
288  continue;
289 
290  QVector< QgsPoint > linePoints;
291  for ( auto point = g.vertices_begin(); point != g.vertices_end(); ++point )
292  {
293  QgsPoint p = *point;
294  double z = 0;
295  switch ( source )
296  {
297  case ValueAttribute:
298  z = attributeValue;
299  break;
300 
301  case ValueZ:
302  z = p.z();
303  break;
304 
305  case ValueM:
306  z = p.m();
307  break;
308  }
309 
310  linePoints.append( QgsPoint( p.x(), p.y(), z ) );
311  }
312  mTriangulation->addLine( linePoints, type );
313  }
314  break;
315  }
318  break;
319  }
320  break;
321  }
322  }
323 
324  return 0;
325 }
326 
327 
328 int QgsTinInterpolator::addPointsFromGeometry( const QgsGeometry &g, ValueSource source, double attributeValue )
329 {
330  // loop through all vertices and add to triangulation
331  for ( auto point = g.vertices_begin(); point != g.vertices_end(); ++point )
332  {
333  QgsPoint p = *point;
334  double z = 0;
335  switch ( source )
336  {
337  case ValueAttribute:
338  z = attributeValue;
339  break;
340 
341  case ValueZ:
342  z = p.z();
343  break;
344 
345  case ValueM:
346  z = p.m();
347  break;
348  }
349  if ( mTriangulation->addPoint( QgsPoint( p.x(), p.y(), z ) ) == -100 )
350  {
351  return -1;
352  }
353  }
354  return 0;
355 }
QgsCurve
Abstract base class for curved geometry type.
Definition: qgscurve.h:35
QgsTinInterpolator::triangulationFields
static QgsFields triangulationFields()
Returns the fields output by features when saving the triangulation.
Definition: qgstininterpolator.cpp:70
QgsFeedback::setProgress
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:75
QgsTinInterpolator::CloughTocher
@ CloughTocher
Clough-Tocher interpolation.
Definition: qgstininterpolator.h:45
Triangulation::addLine
virtual void addLine(const QVector< QgsPoint > &points, QgsInterpolator::SourceType lineType)=0
Adds a line (e.g.
QgsInterpolator::SourceStructureLines
@ SourceStructureLines
Structure lines.
Definition: qgsinterpolator.h:74
QgsInterpolator::ValueAttribute
@ ValueAttribute
Take value from feature's attribute.
Definition: qgsinterpolator.h:81
DualEdgeTriangulation::saveTriangulation
bool saveTriangulation(QgsFeatureSink *sink, QgsFeedback *feedback=nullptr) const override
Saves the triangulation features to a feature sink.
Definition: DualEdgeTriangulation.cpp:2973
QgsWkbTypes::NullGeometry
@ NullGeometry
Definition: qgswkbtypes.h:145
NormVecDecorator.h
QgsPoint
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
QgsTinInterpolator::~QgsTinInterpolator
~QgsTinInterpolator() override
Definition: qgstininterpolator.cpp:43
qgswkbptr.h
QgsGeometryCollection::numGeometries
int numGeometries() const
Returns the number of geometries within the collection.
Definition: qgsgeometrycollection.h:51
QgsGeometry::isMultipart
bool isMultipart() const
Returns true if WKB of the geometry is of WKBMulti* type.
Definition: qgsgeometry.cpp:377
QgsInterpolator::SourcePoints
@ SourcePoints
Point source.
Definition: qgsinterpolator.h:73
crs
const QgsCoordinateReferenceSystem & crs
Definition: qgswfsgetfeature.cpp:105
QgsCurvePolygon
Curve polygon geometry type.
Definition: qgscurvepolygon.h:34
qgsfeatureiterator.h
QgsFields
Definition: qgsfields.h:44
QgsInterpolator::ValueZ
@ ValueZ
Use feature's geometry Z values for interpolation.
Definition: qgsinterpolator.h:82
qgsfeature.h
QgsFeature::geometry
QgsGeometry geometry
Definition: qgsfeature.h:71
QgsGeometry::wkbType
QgsWkbTypes::Type wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
Definition: qgsgeometry.cpp:345
QgsPoint::z
double z
Definition: qgspoint.h:60
QgsInterpolator
Definition: qgsinterpolator.h:66
LinTriangleInterpolator.h
qgspoint.h
qgstininterpolator.h
QgsMultiSurface
Multi surface geometry collection.
Definition: qgsmultisurface.h:29
TriDecorator::addTriangulation
virtual void addTriangulation(Triangulation *t)
Adds an association to a triangulation.
Definition: TriDecorator.h:78
QgsInterpolator::mLayerData
QList< LayerData > mLayerData
Information about the input vector layers and the attributes (or z-values) that are used for interpol...
Definition: qgsinterpolator.h:152
QgsAttributeList
QList< int > QgsAttributeList
Definition: qgsfield.h:26
QgsWkbTypes::PolygonGeometry
@ PolygonGeometry
Definition: qgswkbtypes.h:143
Triangulation::triangulationFields
static QgsFields triangulationFields()
Returns the fields output by features when calling saveTriangulation().
Definition: Triangulation.cpp:18
DualEdgeTriangulation
Definition: DualEdgeTriangulation.h:47
QgsInterpolator::ValueSource
ValueSource
Source for interpolated values from features.
Definition: qgsinterpolator.h:79
QgsPoint::y
double y
Definition: qgspoint.h:59
QgsFeatureRequest
Definition: qgsfeaturerequest.h:75
QgsTinInterpolator::setTriangulationSink
void setTriangulationSink(QgsFeatureSink *sink)
Sets the optional sink for saving the triangulation features.
Definition: qgstininterpolator.cpp:75
LinTriangleInterpolator
Definition: LinTriangleInterpolator.h:31
QgsMultiCurve
Multi curve geometry collection.
Definition: qgsmulticurve.h:29
QgsTinInterpolator::TinInterpolation
TinInterpolation
Indicates the type of interpolation to be performed.
Definition: qgstininterpolator.h:42
QgsFeedback
Definition: qgsfeedback.h:43
CloughTocherInterpolator::setTriangulation
virtual void setTriangulation(NormVecDecorator *tin)
Definition: CloughTocherInterpolator.cpp:28
NormVecDecorator::estimateFirstDerivatives
bool estimateFirstDerivatives(QgsFeedback *feedback=nullptr)
This method adds the functionality of estimating normals at the data points. Return true in the case ...
Definition: NormVecDecorator.cpp:494
QgsWkbTypes::geometryType
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:937
QgsPoint::m
double m
Definition: qgspoint.h:61
QgsFeature::attribute
QVariant attribute(const QString &name) const
Lookup attribute value from attribute name.
Definition: qgsfeature.cpp:262
QgsAbstractGeometry::is3D
bool is3D() const
Returns true if the geometry is 3D and contains a z-value.
Definition: qgsabstractgeometry.h:202
QgsGeometry::isEmpty
bool isEmpty() const
Returns true if the geometry is empty (eg a linestring with no vertices, or a collection with no geom...
Definition: qgsgeometry.cpp:367
QgsGeometry::isNull
bool isNull
Definition: qgsgeometry.h:125
QgsGeometry::constGet
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
Definition: qgsgeometry.cpp:128
CloughTocherInterpolator.h
CloughTocherInterpolator
Definition: CloughTocherInterpolator.h:33
qgscurvepolygon.h
Triangulation::addPoint
virtual int addPoint(const QgsPoint &point)=0
Adds a point to the triangulation.
QgsTinInterpolator::QgsTinInterpolator
QgsTinInterpolator(const QList< QgsInterpolator::LayerData > &inputData, TinInterpolation interpolation=Linear, QgsFeedback *feedback=nullptr)
Constructor for QgsTinInterpolator.
Definition: qgstininterpolator.cpp:35
QgsCoordinateReferenceSystem
Definition: qgscoordinatereferencesystem.h:206
qgsvectorlayer.h
QgsFeedback::isCanceled
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:66
QgsInterpolator::SourceBreakLines
@ SourceBreakLines
Break lines.
Definition: qgsinterpolator.h:75
QgsWkbTypes::LineGeometry
@ LineGeometry
Definition: qgswkbtypes.h:142
QgsWkbTypes::PointGeometry
@ PointGeometry
Definition: qgswkbtypes.h:141
DualEdgeTriangulation.h
QgsGeometryCollection::geometryN
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
Definition: qgsgeometrycollection.h:79
qgsgeometry.h
QgsFeatureIterator::nextFeature
bool nextFeature(QgsFeature &f)
Definition: qgsfeatureiterator.h:373
qgscurve.h
QgsGeometry
Definition: qgsgeometry.h:122
TriangleInterpolator::calcPoint
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.
QgsWkbTypes::UnknownGeometry
@ UnknownGeometry
Definition: qgswkbtypes.h:144
QgsAbstractGeometry::isMeasure
bool isMeasure() const
Returns true if the geometry contains m values.
Definition: qgsabstractgeometry.h:211
NormVecDecorator::setTriangleInterpolator
void setTriangleInterpolator(TriangleInterpolator *inter) override
Sets an interpolator.
Definition: NormVecDecorator.h:106
QgsInterpolator::ValueM
@ ValueM
Use feature's geometry M values for interpolation.
Definition: qgsinterpolator.h:83
qgsmulticurve.h
QgsFeature
Definition: qgsfeature.h:55
QgsTinInterpolator::interpolatePoint
int interpolatePoint(double x, double y, double &result, QgsFeedback *feedback) override
Calculates interpolation value for map coordinates x, y.
Definition: qgstininterpolator.cpp:49
QgsGeometry::vertices_end
QgsAbstractGeometry::vertex_iterator vertices_end() const
Returns STL-style iterator pointing to the imaginary vertex after the last vertex of the geometry.
Definition: qgsgeometry.cpp:1858
QgsGeometry::vertices_begin
QgsAbstractGeometry::vertex_iterator vertices_begin() const
Returns STL-style iterator pointing to the first vertex of the geometry.
Definition: qgsgeometry.cpp:1851
QgsFeatureIterator
Definition: qgsfeatureiterator.h:263
qgsfeedback.h
QgsPoint::x
double x
Definition: qgspoint.h:58
qgsmultisurface.h
QgsFeatureSink
Definition: qgsfeaturesink.h:33
NormVecDecorator
Definition: NormVecDecorator.h:36