QGIS API Documentation 3.37.0-Master (fdefdf9c27f)
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"
22#include "NormVecDecorator.h"
24#include "qgspoint.h"
25#include "qgsfeature.h"
26#include "qgsgeometry.h"
27#include "qgsvectorlayer.h"
28#include "qgsvariantutils.h"
29#include "qgsfeedback.h"
30#include "qgscurve.h"
31#include "qgsmulticurve.h"
32#include "qgscurvepolygon.h"
33#include "qgsmultisurface.h"
34
35QgsTinInterpolator::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
49int 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
80void QgsTinInterpolator::initialize()
81{
82 QgsDualEdgeTriangulation *dualEdgeTriangulation = new QgsDualEdgeTriangulation( 100000 );
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 long long nFeatures = 0;
96 long long nProcessedFeatures = 0;
97 if ( mFeedback )
98 {
99 for ( const LayerData &layer : std::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 : std::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 NormVecDecorator *dec = dynamic_cast<NormVecDecorator *>( mTriangulation );
149 if ( dec )
150 {
151 auto ctInterpolator = std::make_unique<CloughTocherInterpolator>();
152 dec->estimateFirstDerivatives( mFeedback );
153 ctInterpolator->setTriangulation( dec );
154 mTriangleInterpolator = ctInterpolator.release();
155 dec->setTriangleInterpolator( mTriangleInterpolator );
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
171int 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 ( QgsVariantUtils::isNull( attributeVariant ) ) //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( mc->curveN( 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 QgsPointSequence linePoints;
291 curve->points( linePoints );
292 for ( QgsPoint &point : linePoints )
293 {
294 switch ( source )
295 {
296 case ValueAttribute:
297 if ( point.is3D() )
298 point.setZ( attributeValue );
299 else
300 point.addZValue( attributeValue );
301 break;
302
303 case ValueM:
304 if ( point.is3D() )
305 point.setZ( point.m() );
306 else
307 point.addZValue( point.m() );
308 break;
309
310 case ValueZ:
311 break;
312 }
313 }
314 mTriangulation->addLine( linePoints, type );
315 }
316 break;
317 }
320 break;
321 }
322 break;
323 }
324 }
325
326 return 0;
327}
328
329
330int QgsTinInterpolator::addPointsFromGeometry( const QgsGeometry &g, ValueSource source, double attributeValue )
331{
332 // loop through all vertices and add to triangulation
333 for ( auto point = g.vertices_begin(); point != g.vertices_end(); ++point )
334 {
335 QgsPoint p = *point;
336 double z = 0;
337 switch ( source )
338 {
339 case ValueAttribute:
340 z = attributeValue;
341 break;
342
343 case ValueZ:
344 z = p.z();
345 break;
346
347 case ValueM:
348 z = p.m();
349 break;
350 }
351 if ( mTriangulation->addPoint( QgsPoint( p.x(), p.y(), z ) ) == -100 )
352 {
353 return -1;
354 }
355 }
356 return 0;
357}
LinTriangleInterpolator is a class which interpolates linearly on a triangulation.
Decorator class which adds the functionality of estimating normals at the data points.
void setTriangleInterpolator(TriangleInterpolator *inter) override
Sets an interpolator.
bool estimateFirstDerivatives(QgsFeedback *feedback=nullptr)
This method adds the functionality of estimating normals at the data points. Return true in the case ...
@ Polygon
Polygons.
@ Unknown
Unknown types.
@ Null
No geometry.
bool isMeasure() const
Returns true if the geometry contains m values.
bool is3D() const
Returns true if the geometry is 3D and contains a z-value.
This class represents a coordinate reference system (CRS).
Curve polygon geometry type.
Abstract base class for curved geometry type.
Definition: qgscurve.h:35
DualEdgeTriangulation is an implementation of a triangulation class based on the dual edge data struc...
bool saveTriangulation(QgsFeatureSink *sink, QgsFeedback *feedback=nullptr) const override
Saves the triangulation features to a feature sink.
Wrapper for iterator of features from vector data provider or vector layer.
bool nextFeature(QgsFeature &f)
Fetch next feature and stores in f, returns true on success.
This class wraps a request for features to a vector layer (or directly its vector data provider).
An interface for objects which accept features via addFeature(s) methods.
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition: qgsfeature.h:56
QgsGeometry geometry
Definition: qgsfeature.h:67
QVariant attribute(const QString &name) const
Lookup attribute value by attribute name.
Definition: qgsfeature.cpp:335
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition: qgsfeedback.h:44
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:53
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:61
Container of fields for a vector layer.
Definition: qgsfields.h:45
int numGeometries() const
Returns the number of geometries within the collection.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:162
Q_GADGET bool isNull
Definition: qgsgeometry.h:164
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
bool isMultipart() const
Returns true if WKB of the geometry is of WKBMulti* type.
bool isEmpty() const
Returns true if the geometry is empty (eg a linestring with no vertices, or a collection with no geom...
QgsAbstractGeometry::vertex_iterator vertices_begin() const
Returns STL-style iterator pointing to the first vertex of the geometry.
Qgis::WkbType wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
QgsAbstractGeometry::vertex_iterator vertices_end() const
Returns STL-style iterator pointing to the imaginary vertex after the last vertex of the geometry.
Interface class for interpolations.
ValueSource
Source for interpolated values from features.
@ ValueM
Use feature's geometry M values for interpolation.
@ ValueAttribute
Take value from feature's attribute.
@ ValueZ
Use feature's geometry Z values for interpolation.
@ SourcePoints
Point source.
@ SourceStructureLines
Structure lines.
@ SourceBreakLines
Break lines.
QList< LayerData > mLayerData
Information about the input vector layers and the attributes (or z-values) that are used for interpol...
Multi curve geometry collection.
Definition: qgsmulticurve.h:29
QgsCurve * curveN(int index)
Returns the curve with the specified index.
Multi surface geometry collection.
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:49
Q_GADGET double x
Definition: qgspoint.h:52
double z
Definition: qgspoint.h:54
double m
Definition: qgspoint.h:55
double y
Definition: qgspoint.h:53
QgsTinInterpolator(const QList< QgsInterpolator::LayerData > &inputData, TinInterpolation interpolation=Linear, QgsFeedback *feedback=nullptr)
Constructor for QgsTinInterpolator.
TinInterpolation
Indicates the type of interpolation to be performed.
@ CloughTocher
Clough-Tocher interpolation.
int interpolatePoint(double x, double y, double &result, QgsFeedback *feedback) override
Calculates interpolation value for map coordinates x, y.
static QgsFields triangulationFields()
Returns the fields output by features when saving the triangulation.
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().
virtual void addLine(const QgsPointSequence &points, QgsInterpolator::SourceType lineType)=0
Adds a line (e.g.
virtual int addPoint(const QgsPoint &point)=0
Adds a point to the triangulation.
static bool isNull(const QVariant &variant, bool silenceNullWarnings=false)
Returns true if the specified variant should be considered a NULL value.
static Qgis::GeometryType geometryType(Qgis::WkbType type)
Returns the geometry type for a WKB type, e.g., both MultiPolygon and CurvePolygon would have a Polyg...
Definition: qgswkbtypes.h:862
virtual void addTriangulation(QgsTriangulation *t)
Adds an association to a triangulation.
Definition: TriDecorator.h:76
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.
QVector< QgsPoint > QgsPointSequence
QList< int > QgsAttributeList
Definition: qgsfield.h:27
const QgsCoordinateReferenceSystem & crs