QGIS API Documentation 3.40.0-Bratislava (b56115d8743)
Loading...
Searching...
No Matches
qgsvectorlayerprofilegenerator.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsvectorlayerprofilegenerator.cpp
3 ---------------
4 begin : March 2022
5 copyright : (C) 2022 by Nyall Dawson
6 email : nyall dot dawson at gmail dot com
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 ***************************************************************************/
18#include "qgsabstractgeometry.h"
20#include "qgsprofilerequest.h"
21#include "qgscurve.h"
22#include "qgsvectorlayer.h"
25#include "qgsgeos.h"
27#include "qgsterrainprovider.h"
28#include "qgspolygon.h"
29#include "qgstessellator.h"
30#include "qgsmultipolygon.h"
31#include "qgsmeshlayerutils.h"
32#include "qgsmultipoint.h"
33#include "qgsmultilinestring.h"
34#include "qgslinesymbol.h"
35#include "qgsfillsymbol.h"
36#include "qgsmarkersymbol.h"
37#include "qgsprofilepoint.h"
38#include "qgsprofilesnapping.h"
40#include <QPolygonF>
41
42//
43// QgsVectorLayerProfileResults
44//
45
47{
48 return QStringLiteral( "vector" );
49}
50
52{
53 QVector<QgsGeometry> res;
54 res.reserve( features.size() );
55 for ( auto it = features.constBegin(); it != features.constEnd(); ++it )
56 {
57 for ( const Feature &feature : it.value() )
58 {
59 res.append( feature.geometry );
60 }
61 }
62 return res;
63}
64
65QVector<QgsAbstractProfileResults::Feature> QgsVectorLayerProfileResults::asFeatures( Qgis::ProfileExportType type, QgsFeedback *feedback ) const
66{
67 switch ( profileType )
68 {
71 return asIndividualFeatures( type, feedback );
72 // distance vs elevation table results are always handled like a continuous surface
73 [[fallthrough]];
74
77 }
79}
80
82{
83 switch ( profileType )
84 {
86 return snapPointToIndividualFeatures( point, context );
88 return QgsAbstractProfileSurfaceResults::snapPoint( point, context );
89 }
91}
92
93
94QVector<QgsProfileIdentifyResults> QgsVectorLayerProfileResults::identify( const QgsDoubleRange &distanceRange, const QgsDoubleRange &elevationRange, const QgsProfileIdentifyContext & )
95{
96 QgsFeatureIds ids;
97 auto visitFeature = [&ids]( QgsFeatureId featureId )
98 {
99 ids << featureId;
100 };
101
102 visitFeaturesInRange( distanceRange, elevationRange, visitFeature );
103 if ( ids.empty() )
104 return {};
105
106 QVector< QVariantMap> idsList;
107 for ( auto it = ids.constBegin(); it != ids.constEnd(); ++it )
108 idsList.append( QVariantMap( {{QStringLiteral( "id" ), *it}} ) );
109
110 return { QgsProfileIdentifyResults( mLayer, idsList ) };
111}
112
113QVector<QgsProfileIdentifyResults> QgsVectorLayerProfileResults::identify( const QgsProfilePoint &point, const QgsProfileIdentifyContext &context )
114{
115 QHash< QgsFeatureId, QVariantMap > features;
116 auto visitFeature = [&features]( QgsFeatureId featureId, double delta, double distance, double elevation )
117 {
118 auto it = features.find( featureId );
119 if ( it == features.end() )
120 {
121 features[ featureId ] = QVariantMap( {{QStringLiteral( "id" ), featureId },
122 {QStringLiteral( "delta" ), delta },
123 {QStringLiteral( "distance" ), distance },
124 {QStringLiteral( "elevation" ), elevation }
125 } );
126 }
127 else
128 {
129 const double currentDelta = it.value().value( QStringLiteral( "delta" ) ).toDouble();
130 if ( delta < currentDelta )
131 {
132 *it = QVariantMap( {{QStringLiteral( "id" ), featureId },
133 {QStringLiteral( "delta" ), delta },
134 {QStringLiteral( "distance" ), distance },
135 {QStringLiteral( "elevation" ), elevation }
136 } );
137 }
138 }
139 };
140
141 visitFeaturesAtPoint( point, context.maximumPointDistanceDelta, context.maximumPointElevationDelta, context.maximumSurfaceElevationDelta, visitFeature, true );
142
143 QVector< QVariantMap> attributes;
144 for ( auto it = features.constBegin(); it != features.constEnd(); ++it )
145 attributes.append( *it );
146
147 QVector<QgsProfileIdentifyResults> res;
148
149 if ( !attributes.empty() )
150 res.append( QgsProfileIdentifyResults( mLayer, attributes ) );
151
153 {
154 const QVector<QgsProfileIdentifyResults> surfaceResults = QgsAbstractProfileSurfaceResults::identify( point, context );
155 res.reserve( surfaceResults.size() );
156 for ( const QgsProfileIdentifyResults &surfaceResult : surfaceResults )
157 {
158 res.append( QgsProfileIdentifyResults( mLayer, surfaceResult.results() ) );
159 }
160 }
161
162 return res;
163}
164
165QgsProfileSnapResult QgsVectorLayerProfileResults::snapPointToIndividualFeatures( const QgsProfilePoint &point, const QgsProfileSnapContext &context )
166{
168 double bestSnapDistance = std::numeric_limits< double >::max();
169
170 auto visitFeature = [&bestSnapDistance, &res]( QgsFeatureId, double delta, double distance, double elevation )
171 {
172 if ( distance < bestSnapDistance )
173 {
174 bestSnapDistance = delta;
175 res.snappedPoint = QgsProfilePoint( distance, elevation );
176 }
177 };
178
179 visitFeaturesAtPoint( point, context.maximumPointDistanceDelta, context.maximumPointElevationDelta, context.maximumSurfaceElevationDelta, visitFeature, false );
180
181 return res;
182}
183
184void QgsVectorLayerProfileResults::visitFeaturesAtPoint( const QgsProfilePoint &point, double maximumPointDistanceDelta, double maximumPointElevationDelta, double maximumSurfaceElevationDelta, const std::function< void( QgsFeatureId, double delta, double distance, double elevation ) > &visitor, bool visitWithin )
185{
186 // TODO -- add spatial index if performance is an issue
187
188 const QgsPoint targetPoint( point.distance(), point.elevation() );
189
190 for ( auto it = features.constBegin(); it != features.constEnd(); ++it )
191 {
192 for ( const Feature &feature : it.value() )
193 {
194 const QgsRectangle featureBounds = feature.crossSectionGeometry.boundingBox();
195 if ( ( featureBounds.xMinimum() - maximumPointDistanceDelta <= point.distance() ) && ( featureBounds.xMaximum() + maximumPointDistanceDelta >= point.distance() ) )
196 {
197 switch ( feature.crossSectionGeometry.type() )
198 {
200 {
201 for ( auto partIt = feature.crossSectionGeometry.const_parts_begin(); partIt != feature.crossSectionGeometry.const_parts_end(); ++partIt )
202 {
203 if ( const QgsPoint *candidatePoint = qgsgeometry_cast< const QgsPoint * >( *partIt ) )
204 {
205 const double snapDistanceDelta = std::fabs( point.distance() - candidatePoint->x() );
206 if ( snapDistanceDelta > maximumPointDistanceDelta )
207 continue;
208
209 const double snapHeightDelta = std::fabs( point.elevation() - candidatePoint->y() );
210 if ( snapHeightDelta > maximumPointElevationDelta )
211 continue;
212
213 const double snapDistance = candidatePoint->distance( targetPoint );
214 visitor( feature.featureId, snapDistance, candidatePoint->x(), candidatePoint->y() );
215 }
216 }
217 break;
218 }
219
221 {
222 for ( auto partIt = feature.crossSectionGeometry.const_parts_begin(); partIt != feature.crossSectionGeometry.const_parts_end(); ++partIt )
223 {
224 if ( const QgsCurve *line = qgsgeometry_cast< const QgsCurve * >( *partIt ) )
225 {
226 // might be a vertical line
227 if ( const QgsLineString *lineString = qgsgeometry_cast< const QgsLineString * >( line ) )
228 {
229 if ( lineString->numPoints() == 2 && qgsDoubleNear( lineString->pointN( 0 ).x(), lineString->pointN( 1 ).x() ) )
230 {
231 const double snapDistanceDelta = std::fabs( point.distance() - lineString->pointN( 0 ).x() );
232 if ( snapDistanceDelta > maximumPointDistanceDelta )
233 continue;
234
235 const double snapHeightDelta = std::fabs( point.elevation() - lineString->pointN( 0 ).y() );
236 if ( snapHeightDelta <= maximumPointElevationDelta )
237 {
238 const double snapDistanceP1 = lineString->pointN( 0 ).distance( targetPoint );
239 visitor( feature.featureId, snapDistanceP1, lineString->pointN( 0 ).x(), lineString->pointN( 0 ).y() );
240 }
241
242 const double snapHeightDelta2 = std::fabs( point.elevation() - lineString->pointN( 1 ).y() );
243 if ( snapHeightDelta2 <= maximumPointElevationDelta )
244 {
245 const double snapDistanceP2 = lineString->pointN( 1 ).distance( targetPoint );
246 visitor( feature.featureId, snapDistanceP2, lineString->pointN( 1 ).x(), lineString->pointN( 1 ).y() );
247 }
248
249 if ( visitWithin )
250 {
251 double elevation1 = lineString->pointN( 0 ).y();
252 double elevation2 = lineString->pointN( 1 ).y();
253 if ( elevation1 > elevation2 )
254 std::swap( elevation1, elevation2 );
255
256 if ( point.elevation() > elevation1 && point.elevation() < elevation2 )
257 {
258 const double snapDistance = std::fabs( lineString->pointN( 0 ).x() - point.distance() );
259 visitor( feature.featureId, snapDistance, lineString->pointN( 0 ).x(), point.elevation() );
260 }
261 }
262 continue;
263 }
264 }
265
266 const QgsRectangle partBounds = ( *partIt )->boundingBox();
267 if ( point.distance() < partBounds.xMinimum() - maximumPointDistanceDelta || point.distance() > partBounds.xMaximum() + maximumPointDistanceDelta )
268 continue;
269
270 const double snappedDistance = point.distance() < partBounds.xMinimum() ? partBounds.xMinimum()
271 : point.distance() > partBounds.xMaximum() ? partBounds.xMaximum() : point.distance();
272
273 const QgsGeometry cutLine( new QgsLineString( QgsPoint( snappedDistance, qgsDoubleNear( minZ, maxZ ) ? minZ - 1 : minZ ), QgsPoint( snappedDistance, maxZ ) ) );
274 QgsGeos cutLineGeos( cutLine.constGet() );
275
276 const QgsGeometry points( cutLineGeos.intersection( line ) );
277
278 for ( auto vertexIt = points.vertices_begin(); vertexIt != points.vertices_end(); ++vertexIt )
279 {
280 const double snapHeightDelta = std::fabs( point.elevation() - ( *vertexIt ).y() );
281 if ( snapHeightDelta > maximumSurfaceElevationDelta )
282 continue;
283
284 const double snapDistance = ( *vertexIt ).distance( targetPoint );
285 visitor( feature.featureId, snapDistance, ( *vertexIt ).x(), ( *vertexIt ).y() );
286 }
287 }
288 }
289 break;
290 }
291
293 {
294 if ( visitWithin )
295 {
296 if ( feature.crossSectionGeometry.intersects( QgsGeometry::fromPointXY( QgsPointXY( point.distance(), point.elevation() ) ) ) )
297 {
298 visitor( feature.featureId, 0, point.distance(), point.elevation() );
299 break;
300 }
301 }
302 for ( auto partIt = feature.crossSectionGeometry.const_parts_begin(); partIt != feature.crossSectionGeometry.const_parts_end(); ++partIt )
303 {
304 if ( const QgsCurve *exterior = qgsgeometry_cast< const QgsPolygon * >( *partIt )->exteriorRing() )
305 {
306 const QgsRectangle partBounds = ( *partIt )->boundingBox();
307 if ( point.distance() < partBounds.xMinimum() - maximumPointDistanceDelta || point.distance() > partBounds.xMaximum() + maximumPointDistanceDelta )
308 continue;
309
310 const double snappedDistance = point.distance() < partBounds.xMinimum() ? partBounds.xMinimum()
311 : point.distance() > partBounds.xMaximum() ? partBounds.xMaximum() : point.distance();
312
313 const QgsGeometry cutLine( new QgsLineString( QgsPoint( snappedDistance, qgsDoubleNear( minZ, maxZ ) ? minZ - 1 : minZ ), QgsPoint( snappedDistance, maxZ ) ) );
314 QgsGeos cutLineGeos( cutLine.constGet() );
315
316 const QgsGeometry points( cutLineGeos.intersection( exterior ) );
317 for ( auto vertexIt = points.vertices_begin(); vertexIt != points.vertices_end(); ++vertexIt )
318 {
319 const double snapHeightDelta = std::fabs( point.elevation() - ( *vertexIt ).y() );
320 if ( snapHeightDelta > maximumSurfaceElevationDelta )
321 continue;
322
323 const double snapDistance = ( *vertexIt ).distance( targetPoint );
324 visitor( feature.featureId, snapDistance, ( *vertexIt ).x(), ( *vertexIt ).y() );
325 }
326 }
327 }
328 break;
329 }
332 break;
333 }
334 }
335 }
336 }
337}
338
339void QgsVectorLayerProfileResults::visitFeaturesInRange( const QgsDoubleRange &distanceRange, const QgsDoubleRange &elevationRange, const std::function<void ( QgsFeatureId )> &visitor )
340{
341 // TODO -- add spatial index if performance is an issue
342 const QgsRectangle profileRange( distanceRange.lower(), elevationRange.lower(), distanceRange.upper(), elevationRange.upper() );
343 const QgsGeometry profileRangeGeometry = QgsGeometry::fromRect( profileRange );
344 QgsGeos profileRangeGeos( profileRangeGeometry.constGet() );
345 profileRangeGeos.prepareGeometry();
346
347 for ( auto it = features.constBegin(); it != features.constEnd(); ++it )
348 {
349 for ( const Feature &feature : it.value() )
350 {
351 if ( feature.crossSectionGeometry.boundingBoxIntersects( profileRange ) )
352 {
353 switch ( feature.crossSectionGeometry.type() )
354 {
356 {
357 for ( auto partIt = feature.crossSectionGeometry.const_parts_begin(); partIt != feature.crossSectionGeometry.const_parts_end(); ++partIt )
358 {
359 if ( const QgsPoint *candidatePoint = qgsgeometry_cast< const QgsPoint * >( *partIt ) )
360 {
361 if ( profileRange.contains( candidatePoint->x(), candidatePoint->y() ) )
362 {
363 visitor( feature.featureId );
364 }
365 }
366 }
367 break;
368 }
369
372 {
373 if ( profileRangeGeos.intersects( feature.crossSectionGeometry.constGet() ) )
374 {
375 visitor( feature.featureId );
376 }
377 break;
378 }
379
382 break;
383 }
384 }
385 }
386 }
387}
388
390{
391 const QgsExpressionContextScopePopper scopePopper( context.renderContext().expressionContext(), mLayer ? mLayer->createExpressionContextScope() : nullptr );
392 switch ( profileType )
393 {
395 renderResultsAsIndividualFeatures( context );
396 break;
400 renderMarkersOverContinuousSurfacePlot( context );
401 break;
402 }
403}
404
405void QgsVectorLayerProfileResults::renderResultsAsIndividualFeatures( QgsProfileRenderContext &context )
406{
407 QPainter *painter = context.renderContext().painter();
408 if ( !painter )
409 return;
410
411 const QgsScopedQPainterState painterState( painter );
412
413 painter->setBrush( Qt::NoBrush );
414 painter->setPen( Qt::NoPen );
415
416 const double minDistance = context.distanceRange().lower();
417 const double maxDistance = context.distanceRange().upper();
418 const double minZ = context.elevationRange().lower();
419 const double maxZ = context.elevationRange().upper();
420
421 const QRectF visibleRegion( minDistance, minZ, maxDistance - minDistance, maxZ - minZ );
422 QPainterPath clipPath;
423 clipPath.addPolygon( context.worldTransform().map( visibleRegion ) );
424 painter->setClipPath( clipPath, Qt::ClipOperation::IntersectClip );
425
426 const QgsRectangle clipPathRect( clipPath.boundingRect() );
427
428 auto renderResult = [&context, &clipPathRect]( const Feature & profileFeature, QgsMarkerSymbol * markerSymbol, QgsLineSymbol * lineSymbol, QgsFillSymbol * fillSymbol )
429 {
430 if ( profileFeature.crossSectionGeometry.isEmpty() )
431 return;
432
433 QgsGeometry transformed = profileFeature.crossSectionGeometry;
434 transformed.transform( context.worldTransform() );
435
436 if ( !transformed.boundingBoxIntersects( clipPathRect ) )
437 return;
438
439 // we can take some shortcuts here, because we know that the geometry will already be segmentized and can't be a curved type
440 switch ( transformed.type() )
441 {
443 {
444 if ( const QgsPoint *point = qgsgeometry_cast< const QgsPoint * >( transformed.constGet() ) )
445 {
446 markerSymbol->renderPoint( QPointF( point->x(), point->y() ), nullptr, context.renderContext() );
447 }
448 else if ( const QgsMultiPoint *multipoint = qgsgeometry_cast< const QgsMultiPoint * >( transformed.constGet() ) )
449 {
450 const int numGeometries = multipoint->numGeometries();
451 for ( int i = 0; i < numGeometries; ++i )
452 {
453 markerSymbol->renderPoint( QPointF( multipoint->pointN( i )->x(), multipoint->pointN( i )->y() ), nullptr, context.renderContext() );
454 }
455 }
456 break;
457 }
458
460 {
461 if ( const QgsLineString *line = qgsgeometry_cast< const QgsLineString * >( transformed.constGet() ) )
462 {
463 lineSymbol->renderPolyline( line->asQPolygonF(), nullptr, context.renderContext() );
464 }
465 else if ( const QgsMultiLineString *multiLinestring = qgsgeometry_cast< const QgsMultiLineString * >( transformed.constGet() ) )
466 {
467 const int numGeometries = multiLinestring->numGeometries();
468 for ( int i = 0; i < numGeometries; ++i )
469 {
470 lineSymbol->renderPolyline( multiLinestring->lineStringN( i )->asQPolygonF(), nullptr, context.renderContext() );
471 }
472 }
473 break;
474 }
475
477 {
478 if ( const QgsPolygon *polygon = qgsgeometry_cast< const QgsPolygon * >( transformed.constGet() ) )
479 {
480 if ( const QgsCurve *exterior = polygon->exteriorRing() )
481 fillSymbol->renderPolygon( exterior->asQPolygonF(), nullptr, nullptr, context.renderContext() );
482 }
483 else if ( const QgsMultiPolygon *multiPolygon = qgsgeometry_cast< const QgsMultiPolygon * >( transformed.constGet() ) )
484 {
485 const int numGeometries = multiPolygon->numGeometries();
486 for ( int i = 0; i < numGeometries; ++i )
487 {
488 fillSymbol->renderPolygon( multiPolygon->polygonN( i )->exteriorRing()->asQPolygonF(), nullptr, nullptr, context.renderContext() );
489 }
490 }
491 break;
492 }
493
496 return;
497 }
498 };
499
501 req.setFilterFids( qgis::listToSet( features.keys() ) );
502
503 if ( respectLayerSymbology && mLayer && mLayer->renderer() )
504 {
505 std::unique_ptr< QgsFeatureRenderer > renderer( mLayer->renderer()->clone() );
506 renderer->startRender( context.renderContext(), mLayer->fields() );
507
508 // if we are respecting the layer's symbology then we'll fire off a feature request and iterate through
509 // features from the profile, rendering each in turn
510 QSet<QString> attributes = renderer->usedAttributes( context.renderContext() );
511
512 std::unique_ptr< QgsMarkerSymbol > marker( mMarkerSymbol->clone() );
513 std::unique_ptr< QgsLineSymbol > line( mLineSymbol->clone() );
514 std::unique_ptr< QgsFillSymbol > fill( mFillSymbol->clone() );
515 attributes.unite( marker->usedAttributes( context.renderContext() ) );
516 attributes.unite( line->usedAttributes( context.renderContext() ) );
517 attributes.unite( fill->usedAttributes( context.renderContext() ) );
518
519 req.setSubsetOfAttributes( attributes, mLayer->fields() );
520
521 QgsFeature feature;
522 QgsFeatureIterator it = mLayer->getFeatures( req );
523 while ( it.nextFeature( feature ) )
524 {
525 context.renderContext().expressionContext().setFeature( feature );
526 QgsSymbol *rendererSymbol = renderer->symbolForFeature( feature, context.renderContext() );
527 if ( !rendererSymbol )
528 continue;
529
530 marker->setColor( rendererSymbol->color() );
531 marker->setOpacity( rendererSymbol->opacity() );
532 line->setColor( rendererSymbol->color() );
533 line->setOpacity( rendererSymbol->opacity() );
534 fill->setColor( rendererSymbol->color() );
535 fill->setOpacity( rendererSymbol->opacity() );
536
537 marker->startRender( context.renderContext() );
538 line->startRender( context.renderContext() );
539 fill->startRender( context.renderContext() );
540
541 const QVector< Feature > profileFeatures = features.value( feature.id() );
542 for ( const Feature &profileFeature : profileFeatures )
543 {
544 renderResult( profileFeature,
545 rendererSymbol->type() == Qgis::SymbolType::Marker ? qgis::down_cast< QgsMarkerSymbol * >( rendererSymbol ) : marker.get(),
546 rendererSymbol->type() == Qgis::SymbolType::Line ? qgis::down_cast< QgsLineSymbol * >( rendererSymbol ) : line.get(),
547 rendererSymbol->type() == Qgis::SymbolType::Fill ? qgis::down_cast< QgsFillSymbol * >( rendererSymbol ) : fill.get() );
548 }
549
550 marker->stopRender( context.renderContext() );
551 line->stopRender( context.renderContext() );
552 fill->stopRender( context.renderContext() );
553 }
554
555 renderer->stopRender( context.renderContext() );
556 }
557 else if ( mLayer )
558 {
559 QSet<QString> attributes;
560 attributes.unite( mMarkerSymbol->usedAttributes( context.renderContext() ) );
561 attributes.unite( mFillSymbol->usedAttributes( context.renderContext() ) );
562 attributes.unite( mLineSymbol->usedAttributes( context.renderContext() ) );
563
564 mMarkerSymbol->startRender( context.renderContext() );
565 mFillSymbol->startRender( context.renderContext() );
566 mLineSymbol->startRender( context.renderContext() );
567 req.setSubsetOfAttributes( attributes, mLayer->fields() );
568
569 QgsFeature feature;
570 QgsFeatureIterator it = mLayer->getFeatures( req );
571 while ( it.nextFeature( feature ) )
572 {
573 context.renderContext().expressionContext().setFeature( feature );
574 const QVector< Feature > profileFeatures = features.value( feature.id() );
575 for ( const Feature &profileFeature : profileFeatures )
576 {
577 renderResult( profileFeature, mMarkerSymbol.get(), mLineSymbol.get(), mFillSymbol.get() );
578 }
579 }
580 mMarkerSymbol->stopRender( context.renderContext() );
581 mFillSymbol->stopRender( context.renderContext() );
582 mLineSymbol->stopRender( context.renderContext() );
583 }
584}
585
586void QgsVectorLayerProfileResults::renderMarkersOverContinuousSurfacePlot( QgsProfileRenderContext &context )
587{
588 QPainter *painter = context.renderContext().painter();
589 if ( !painter )
590 return;
591
592 const QgsScopedQPainterState painterState( painter );
593
594 painter->setBrush( Qt::NoBrush );
595 painter->setPen( Qt::NoPen );
596
597 const double minDistance = context.distanceRange().lower();
598 const double maxDistance = context.distanceRange().upper();
599 const double minZ = context.elevationRange().lower();
600 const double maxZ = context.elevationRange().upper();
601
602 const QRectF visibleRegion( minDistance, minZ, maxDistance - minDistance, maxZ - minZ );
603 QPainterPath clipPath;
604 clipPath.addPolygon( context.worldTransform().map( visibleRegion ) );
605 painter->setClipPath( clipPath, Qt::ClipOperation::IntersectClip );
606
607 mMarkerSymbol->startRender( context.renderContext() );
608
609 for ( auto pointIt = mDistanceToHeightMap.constBegin(); pointIt != mDistanceToHeightMap.constEnd(); ++pointIt )
610 {
611 if ( std::isnan( pointIt.value() ) )
612 continue;
613
614 mMarkerSymbol->renderPoint( context.worldTransform().map( QPointF( pointIt.key(), pointIt.value() ) ), nullptr, context.renderContext() );
615 }
616 mMarkerSymbol->stopRender( context.renderContext() );
617}
618
619QVector<QgsAbstractProfileResults::Feature> QgsVectorLayerProfileResults::asIndividualFeatures( Qgis::ProfileExportType type, QgsFeedback *feedback ) const
620{
621 QVector<QgsAbstractProfileResults::Feature> res;
622 res.reserve( features.size() );
623 for ( auto it = features.constBegin(); it != features.constEnd(); ++it )
624 {
625 if ( feedback && feedback->isCanceled() )
626 break;
627
628 for ( const Feature &feature : it.value() )
629 {
630 if ( feedback && feedback->isCanceled() )
631 break;
632
634 outFeature.layerIdentifier = mId;
635 outFeature.attributes = {{QStringLiteral( "id" ), feature.featureId }};
636 switch ( type )
637 {
639 outFeature.geometry = feature.geometry;
640 break;
641
643 outFeature.geometry = feature.crossSectionGeometry;
644 break;
645
647 break; // unreachable
648 }
649 res << outFeature;
650 }
651 }
652 return res;
653}
654
656{
658 const QgsVectorLayerProfileGenerator *vlGenerator = qgis::down_cast< const QgsVectorLayerProfileGenerator * >( generator );
659
660 mId = vlGenerator->mId;
661 profileType = vlGenerator->mType;
662 respectLayerSymbology = vlGenerator->mRespectLayerSymbology;
663 mMarkerSymbol.reset( vlGenerator->mProfileMarkerSymbol->clone() );
664 mShowMarkerSymbolInSurfacePlots = vlGenerator->mShowMarkerSymbolInSurfacePlots;
665}
666
667//
668// QgsVectorLayerProfileGenerator
669//
670
673 , mId( layer->id() )
674 , mFeedback( std::make_unique< QgsFeedback >() )
675 , mProfileCurve( request.profileCurve() ? request.profileCurve()->clone() : nullptr )
676 , mTerrainProvider( request.terrainProvider() ? request.terrainProvider()->clone() : nullptr )
677 , mTolerance( request.tolerance() )
678 , mSourceCrs( layer->crs3D() )
679 , mTargetCrs( request.crs() )
680 , mTransformContext( request.transformContext() )
681 , mExtent( layer->extent() )
682 , mSource( std::make_unique< QgsVectorLayerFeatureSource >( layer ) )
683 , mOffset( layer->elevationProperties()->zOffset() )
684 , mScale( layer->elevationProperties()->zScale() )
685 , mType( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->type() )
686 , mClamping( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->clamping() )
687 , mBinding( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->binding() )
688 , mExtrusionEnabled( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->extrusionEnabled() )
689 , mExtrusionHeight( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->extrusionHeight() )
690 , mExpressionContext( request.expressionContext() )
691 , mFields( layer->fields() )
692 , mDataDefinedProperties( layer->elevationProperties()->dataDefinedProperties() )
693 , mWkbType( layer->wkbType() )
694 , mRespectLayerSymbology( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->respectLayerSymbology() )
695 , mProfileMarkerSymbol( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->profileMarkerSymbol()->clone() )
696 , mShowMarkerSymbolInSurfacePlots( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->showMarkerSymbolInSurfacePlots() )
697 , mLayer( layer )
698{
699 if ( mTerrainProvider )
700 mTerrainProvider->prepare(); // must be done on main thread
701
702 // make sure profile curve is always 2d, or we may get unwanted z value averaging for intersections from GEOS
703 if ( mProfileCurve )
704 mProfileCurve->dropZValue();
705
706 mSymbology = qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->profileSymbology();
707 mElevationLimit = qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->elevationLimit();
708
709 mLineSymbol.reset( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->profileLineSymbol()->clone() );
710 mFillSymbol.reset( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->profileFillSymbol()->clone() );
711}
712
714{
715 return mId;
716}
717
719
721{
722 if ( !mProfileCurve || mFeedback->isCanceled() )
723 return false;
724
725 if ( QgsLineString *profileLine =
726 qgsgeometry_cast<QgsLineString *>( mProfileCurve.get() ) )
727 {
728 // The profile generation code can't deal with curves that enter a single
729 // point multiple times. We handle this for line strings by splitting them
730 // into multiple parts, each with no repeated points, and computing the
731 // profile for each by itself.
732 std::unique_ptr< QgsCurve > origCurve = std::move( mProfileCurve );
733 std::unique_ptr< QgsVectorLayerProfileResults > totalResults;
734 double distanceProcessed = 0;
735
736 QVector<QgsLineString *> disjointParts = profileLine->splitToDisjointXYParts();
737 for ( int i = 0; i < disjointParts.size(); i++ )
738 {
739 mProfileCurve.reset( disjointParts[i] );
740 if ( !generateProfileInner() )
741 {
742 mProfileCurve = std::move( origCurve );
743
744 // Free the rest of the parts
745 for ( int j = i + 1; j < disjointParts.size(); j++ )
746 delete disjointParts[j];
747
748 return false;
749 }
750
751 if ( !totalResults )
752 // Use the first result set as a base
753 totalResults.reset( mResults.release() );
754 else
755 {
756 // Merge the results, shifting them by distanceProcessed
757 totalResults->mRawPoints.append( mResults->mRawPoints );
758 totalResults->minZ = std::min( totalResults->minZ, mResults->minZ );
759 totalResults->maxZ = std::max( totalResults->maxZ, mResults->maxZ );
760 for ( auto it = mResults->mDistanceToHeightMap.constKeyValueBegin();
761 it != mResults->mDistanceToHeightMap.constKeyValueEnd();
762 ++it )
763 {
764 totalResults->mDistanceToHeightMap[it->first + distanceProcessed] = it->second;
765 }
766 for ( auto it = mResults->features.constKeyValueBegin();
767 it != mResults->features.constKeyValueEnd();
768 ++it )
769 {
770 for ( QgsVectorLayerProfileResults::Feature feature : it->second )
771 {
772 feature.crossSectionGeometry.translate( distanceProcessed, 0 );
773 totalResults->features[it->first].push_back( feature );
774 }
775 }
776 }
777
778 distanceProcessed += mProfileCurve->length();
779 }
780
781 mProfileCurve = std::move( origCurve );
782 mResults.reset( totalResults.release() );
783 return true;
784 }
785
786 return generateProfileInner();
787}
788
789bool QgsVectorLayerProfileGenerator::generateProfileInner( const QgsProfileGenerationContext & )
790{
791 // we need to transform the profile curve to the vector's CRS
792 mTransformedCurve.reset( mProfileCurve->clone() );
793 mLayerToTargetTransform = QgsCoordinateTransform( mSourceCrs, mTargetCrs, mTransformContext );
794 if ( mTerrainProvider )
795 mTargetToTerrainProviderTransform = QgsCoordinateTransform( mTargetCrs, mTerrainProvider->crs(), mTransformContext );
796
797 try
798 {
799 mTransformedCurve->transform( mLayerToTargetTransform, Qgis::TransformDirection::Reverse );
800 }
801 catch ( QgsCsException & )
802 {
803 QgsDebugError( QStringLiteral( "Error transforming profile line to vector CRS" ) );
804 return false;
805 }
806
807 const QgsRectangle profileCurveBoundingBox = mTransformedCurve->boundingBox();
808 if ( !profileCurveBoundingBox.intersects( mExtent ) )
809 return false;
810
811 if ( mFeedback->isCanceled() )
812 return false;
813
814 mResults = std::make_unique< QgsVectorLayerProfileResults >();
815 mResults->mLayer = mLayer;
816 mResults->copyPropertiesFromGenerator( this );
817
818 mProfileCurveEngine.reset( new QgsGeos( mProfileCurve.get() ) );
819 mProfileCurveEngine->prepareGeometry();
820
821 if ( mTolerance == 0.0 ) // geos does not handle very well buffer with 0 size
822 {
823 mProfileBufferedCurve = std::unique_ptr<QgsAbstractGeometry>( mProfileCurve->clone() );
824 }
825 else
826 {
827 mProfileBufferedCurve = std::unique_ptr<QgsAbstractGeometry>( mProfileCurveEngine->buffer( mTolerance, 8, Qgis::EndCapStyle::Flat, Qgis::JoinStyle::Round, 2 ) );
828 }
829
830 mProfileBufferedCurveEngine.reset( new QgsGeos( mProfileBufferedCurve.get() ) );
831 mProfileBufferedCurveEngine->prepareGeometry();
832
833 mDataDefinedProperties.prepare( mExpressionContext );
834
835 if ( mFeedback->isCanceled() )
836 return false;
837
838 switch ( QgsWkbTypes::geometryType( mWkbType ) )
839 {
841 if ( !generateProfileForPoints() )
842 return false;
843 break;
844
846 if ( !generateProfileForLines() )
847 return false;
848 break;
849
851 if ( !generateProfileForPolygons() )
852 return false;
853 break;
854
857 return false;
858 }
859
860 return true;
861}
862
867
869{
870 return mFeedback.get();
871}
872
873bool QgsVectorLayerProfileGenerator::generateProfileForPoints()
874{
875 // get features from layer
876 QgsFeatureRequest request;
877 request.setCoordinateTransform( QgsCoordinateTransform( mSourceCrs, mTargetCrs, mTransformContext ) );
878 request.setDistanceWithin( QgsGeometry( mProfileCurve->clone() ), mTolerance );
879 request.setSubsetOfAttributes( mDataDefinedProperties.referencedFields( mExpressionContext ), mFields );
880 request.setFeedback( mFeedback.get() );
881
882 // our feature request is using the optimised distance within check (allowing use of spatial index)
883 // BUT this will also include points which are within the tolerance distance before/after the end of line.
884 // So we also need to double check that they fall within the flat buffered curve too.
885
886 QgsFeature feature;
887 QgsFeatureIterator it = mSource->getFeatures( request );
888 while ( !mFeedback->isCanceled() && it.nextFeature( feature ) )
889 {
890 mExpressionContext.setFeature( feature );
891
892 const QgsGeometry g = feature.geometry();
893 for ( auto it = g.const_parts_begin(); !mFeedback->isCanceled() && it != g.const_parts_end(); ++it )
894 {
895 if ( mProfileBufferedCurveEngine->intersects( *it ) )
896 {
897 processIntersectionPoint( qgsgeometry_cast< const QgsPoint * >( *it ), feature );
898 }
899 }
900 }
901 return !mFeedback->isCanceled();
902}
903
904void QgsVectorLayerProfileGenerator::processIntersectionPoint( const QgsPoint *point, const QgsFeature &feature )
905{
906 QString error;
907 const double offset = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::Property::ZOffset, mExpressionContext, mOffset );
908
909 const double height = featureZToHeight( point->x(), point->y(), point->z(), offset );
910 mResults->mRawPoints.append( QgsPoint( point->x(), point->y(), height ) );
911 mResults->minZ = std::min( mResults->minZ, height );
912 mResults->maxZ = std::max( mResults->maxZ, height );
913
914 const double distanceAlongProfileCurve = mProfileCurveEngine->lineLocatePoint( *point, &error );
915 mResults->mDistanceToHeightMap.insert( distanceAlongProfileCurve, height );
916
918 resultFeature.featureId = feature.id();
919 if ( mExtrusionEnabled )
920 {
921 const double extrusion = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::Property::ExtrusionHeight, mExpressionContext, mExtrusionHeight );
922
923 resultFeature.geometry = QgsGeometry( new QgsLineString( QgsPoint( point->x(), point->y(), height ),
924 QgsPoint( point->x(), point->y(), height + extrusion ) ) );
925 resultFeature.crossSectionGeometry = QgsGeometry( new QgsLineString( QgsPoint( distanceAlongProfileCurve, height ),
926 QgsPoint( distanceAlongProfileCurve, height + extrusion ) ) );
927 mResults->minZ = std::min( mResults->minZ, height + extrusion );
928 mResults->maxZ = std::max( mResults->maxZ, height + extrusion );
929 }
930 else
931 {
932 resultFeature.geometry = QgsGeometry( new QgsPoint( point->x(), point->y(), height ) );
933 resultFeature.crossSectionGeometry = QgsGeometry( new QgsPoint( distanceAlongProfileCurve, height ) );
934 }
935
936 mResults->features[resultFeature.featureId].append( resultFeature );
937}
938
939void QgsVectorLayerProfileGenerator::processIntersectionCurve( const QgsLineString *intersectionCurve, const QgsFeature &feature )
940{
941 QString error;
942
944 resultFeature.featureId = feature.id();
945 double maxDistanceAlongProfileCurve = std::numeric_limits<double>::lowest();
946
947 const double offset = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::Property::ZOffset, mExpressionContext, mOffset );
948 const double extrusion = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::Property::ExtrusionHeight, mExpressionContext, mExtrusionHeight );
949
950 const int numPoints = intersectionCurve->numPoints();
951 QVector< double > newX( numPoints );
952 QVector< double > newY( numPoints );
953 QVector< double > newZ( numPoints );
954 QVector< double > newDistance( numPoints );
955
956 const double *inX = intersectionCurve->xData();
957 const double *inY = intersectionCurve->yData();
958 const double *inZ = intersectionCurve->is3D() ? intersectionCurve->zData() : nullptr;
959 double *outX = newX.data();
960 double *outY = newY.data();
961 double *outZ = newZ.data();
962 double *outDistance = newDistance.data();
963
964 QVector< double > extrudedZ;
965 double *extZOut = nullptr;
966 if ( mExtrusionEnabled )
967 {
968 extrudedZ.resize( numPoints );
969 extZOut = extrudedZ.data();
970 }
971
972 for ( int i = 0 ; ! mFeedback->isCanceled() && i < numPoints; ++i )
973 {
974 QgsPoint intersectionPoint( *inX, *inY, ( inZ ? *inZ : std::numeric_limits<double>::quiet_NaN() ) );
975
976 const double height = featureZToHeight( intersectionPoint.x(), intersectionPoint.y(), intersectionPoint.z(), offset );
977 const double distanceAlongProfileCurve = mProfileCurveEngine->lineLocatePoint( intersectionPoint, &error );
978
979 maxDistanceAlongProfileCurve = std::max( maxDistanceAlongProfileCurve, distanceAlongProfileCurve );
980
981 mResults->mRawPoints.append( QgsPoint( intersectionPoint.x(), intersectionPoint.y(), height ) );
982 mResults->minZ = std::min( mResults->minZ, height );
983 mResults->maxZ = std::max( mResults->maxZ, height );
984
985 mResults->mDistanceToHeightMap.insert( distanceAlongProfileCurve, height );
986 *outDistance++ = distanceAlongProfileCurve;
987
988 *outX++ = intersectionPoint.x();
989 *outY++ = intersectionPoint.y();
990 *outZ++ = height;
991 if ( extZOut )
992 *extZOut++ = height + extrusion;
993
994 if ( mExtrusionEnabled )
995 {
996 mResults->minZ = std::min( mResults->minZ, height + extrusion );
997 mResults->maxZ = std::max( mResults->maxZ, height + extrusion );
998 }
999 inX++;
1000 inY++;
1001 if ( inZ )
1002 inZ++;
1003 }
1004
1005 mResults->mDistanceToHeightMap.insert( maxDistanceAlongProfileCurve + 0.000001, std::numeric_limits<double>::quiet_NaN() );
1006
1007 if ( mFeedback->isCanceled() )
1008 return;
1009
1010 // create geometries from vector data
1011 if ( mExtrusionEnabled )
1012 {
1013 std::unique_ptr< QgsLineString > ring = std::make_unique< QgsLineString >( newX, newY, newZ );
1014 std::unique_ptr< QgsLineString > extrudedRing = std::make_unique< QgsLineString >( newX, newY, extrudedZ );
1015 std::unique_ptr< QgsLineString > reversedExtrusion( extrudedRing->reversed() );
1016 ring->append( reversedExtrusion.get() );
1017 ring->close();
1018 resultFeature.geometry = QgsGeometry( new QgsPolygon( ring.release() ) );
1019
1020 std::unique_ptr< QgsLineString > distanceVHeightRing = std::make_unique< QgsLineString >( newDistance, newZ );
1021 std::unique_ptr< QgsLineString > extrudedDistanceVHeightRing = std::make_unique< QgsLineString >( newDistance, extrudedZ );
1022 std::unique_ptr< QgsLineString > reversedDistanceVHeightExtrusion( extrudedDistanceVHeightRing->reversed() );
1023 distanceVHeightRing->append( reversedDistanceVHeightExtrusion.get() );
1024 distanceVHeightRing->close();
1025 resultFeature.crossSectionGeometry = QgsGeometry( new QgsPolygon( distanceVHeightRing.release() ) );
1026 }
1027 else
1028 {
1029 resultFeature.geometry = QgsGeometry( new QgsLineString( newX, newY, newZ ) ) ;
1030 resultFeature.crossSectionGeometry = QgsGeometry( new QgsLineString( newDistance, newZ ) );
1031 }
1032
1033 mResults->features[resultFeature.featureId].append( resultFeature );
1034}
1035
1036bool QgsVectorLayerProfileGenerator::generateProfileForLines()
1037{
1038 // get features from layer
1039 QgsFeatureRequest request;
1040 request.setDestinationCrs( mTargetCrs, mTransformContext );
1041 if ( mTolerance > 0 )
1042 {
1043 request.setDistanceWithin( QgsGeometry( mProfileCurve->clone() ), mTolerance );
1044 }
1045 else
1046 {
1047 request.setFilterRect( mProfileCurve->boundingBox() );
1048 }
1049 request.setSubsetOfAttributes( mDataDefinedProperties.referencedFields( mExpressionContext ), mFields );
1050 request.setFeedback( mFeedback.get() );
1051
1052 auto processCurve = [this]( const QgsFeature & feature, const QgsCurve * featGeomPart )
1053 {
1054 QString error;
1055 std::unique_ptr< QgsAbstractGeometry > intersection( mProfileBufferedCurveEngine->intersection( featGeomPart, &error ) );
1056 if ( !intersection )
1057 return;
1058
1059 if ( mFeedback->isCanceled() )
1060 return;
1061
1062
1063 // Intersection is empty : GEOS issue for vertical intersection : use feature geometry as intersection
1064 if ( intersection->isEmpty() )
1065 {
1066 intersection.reset( featGeomPart->clone() );
1067 }
1068
1069 QgsGeos featGeomPartGeos( featGeomPart );
1070 featGeomPartGeos.prepareGeometry();
1071
1072 for ( auto it = intersection->const_parts_begin();
1073 !mFeedback->isCanceled() && it != intersection->const_parts_end();
1074 ++it )
1075 {
1076 if ( const QgsPoint *intersectionPoint = qgsgeometry_cast< const QgsPoint * >( *it ) )
1077 {
1078 // unfortunately we need to do some work to interpolate the z value for the line -- GEOS doesn't give us this
1079 QString error;
1080 const double distance = featGeomPartGeos.lineLocatePoint( *intersectionPoint, &error );
1081 std::unique_ptr< QgsPoint > interpolatedPoint( featGeomPart->interpolatePoint( distance ) );
1082
1083 processIntersectionPoint( interpolatedPoint.get(), feature );
1084 }
1085 else if ( const QgsLineString *intersectionCurve = qgsgeometry_cast< const QgsLineString * >( *it ) )
1086 {
1087 processIntersectionCurve( intersectionCurve, feature );
1088 }
1089 }
1090 };
1091
1092 QgsFeature feature;
1093 QgsFeatureIterator it = mSource->getFeatures( request );
1094 while ( !mFeedback->isCanceled() && it.nextFeature( feature ) )
1095 {
1096 mExpressionContext.setFeature( feature );
1097
1098 const QgsGeometry g = feature.geometry();
1099 for ( auto it = g.const_parts_begin(); !mFeedback->isCanceled() && it != g.const_parts_end(); ++it )
1100 {
1101 if ( mProfileBufferedCurveEngine->intersects( *it ) )
1102 {
1103 processCurve( feature, qgsgeometry_cast< const QgsCurve * >( *it ) );
1104 }
1105 }
1106 }
1107
1108 return !mFeedback->isCanceled();
1109}
1110
1111QgsPoint QgsVectorLayerProfileGenerator::interpolatePointOnTriangle( const QgsPolygon *triangle, double x, double y ) const
1112{
1113 QgsPoint p1, p2, p3;
1115 triangle->exteriorRing()->pointAt( 0, p1, vt );
1116 triangle->exteriorRing()->pointAt( 1, p2, vt );
1117 triangle->exteriorRing()->pointAt( 2, p3, vt );
1118 const double z = QgsMeshLayerUtils::interpolateFromVerticesData( p1, p2, p3, p1.z(), p2.z(), p3.z(), QgsPointXY( x, y ) );
1119 return QgsPoint( x, y, z );
1120};
1121
1122void QgsVectorLayerProfileGenerator::processTriangleIntersectForPoint( const QgsPolygon *triangle, const QgsPoint *p, QVector< QgsGeometry > &transformedParts, QVector< QgsGeometry > &crossSectionParts )
1123{
1124 const QgsPoint interpolatedPoint = interpolatePointOnTriangle( triangle, p->x(), p->y() );
1125 mResults->mRawPoints.append( interpolatedPoint );
1126 mResults->minZ = std::min( mResults->minZ, interpolatedPoint.z() );
1127 mResults->maxZ = std::max( mResults->maxZ, interpolatedPoint.z() );
1128
1129 QString lastError;
1130 const double distance = mProfileCurveEngine->lineLocatePoint( *p, &lastError );
1131 mResults->mDistanceToHeightMap.insert( distance, interpolatedPoint.z() );
1132
1133 if ( mExtrusionEnabled )
1134 {
1135 const double extrusion = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::Property::ExtrusionHeight, mExpressionContext, mExtrusionHeight );
1136
1137 transformedParts.append( QgsGeometry( new QgsLineString( interpolatedPoint,
1138 QgsPoint( interpolatedPoint.x(), interpolatedPoint.y(), interpolatedPoint.z() + extrusion ) ) ) );
1139 crossSectionParts.append( QgsGeometry( new QgsLineString( QgsPoint( distance, interpolatedPoint.z() ),
1140 QgsPoint( distance, interpolatedPoint.z() + extrusion ) ) ) );
1141 mResults->minZ = std::min( mResults->minZ, interpolatedPoint.z() + extrusion );
1142 mResults->maxZ = std::max( mResults->maxZ, interpolatedPoint.z() + extrusion );
1143 }
1144 else
1145 {
1146 transformedParts.append( QgsGeometry( new QgsPoint( interpolatedPoint ) ) );
1147 crossSectionParts.append( QgsGeometry( new QgsPoint( distance, interpolatedPoint.z() ) ) );
1148 }
1149}
1150
1151void QgsVectorLayerProfileGenerator::processTriangleIntersectForLine( const QgsPolygon *triangle, const QgsLineString *intersectionLine, QVector< QgsGeometry > &transformedParts, QVector< QgsGeometry > &crossSectionParts )
1152{
1153 if ( triangle->exteriorRing()->numPoints() < 4 ) // not a polygon
1154 return;
1155
1156 int numPoints = intersectionLine->numPoints();
1157 QVector< double > newX( numPoints );
1158 QVector< double > newY( numPoints );
1159 QVector< double > newZ( numPoints );
1160 QVector< double > newDistance( numPoints );
1161
1162 const double *inX = intersectionLine->xData();
1163 const double *inY = intersectionLine->yData();
1164 const double *inZ = intersectionLine->is3D() ? intersectionLine->zData() : nullptr;
1165 double *outX = newX.data();
1166 double *outY = newY.data();
1167 double *outZ = newZ.data();
1168 double *outDistance = newDistance.data();
1169
1170 double lastDistanceAlongProfileCurve = 0.0;
1171 QVector< double > extrudedZ;
1172 double *extZOut = nullptr;
1173 double extrusion = 0;
1174
1175 if ( mExtrusionEnabled )
1176 {
1177 extrudedZ.resize( numPoints );
1178 extZOut = extrudedZ.data();
1179
1180 extrusion = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::Property::ExtrusionHeight, mExpressionContext, mExtrusionHeight );
1181 }
1182
1183 QString lastError;
1184 for ( int i = 0 ; ! mFeedback->isCanceled() && i < numPoints; ++i )
1185 {
1186 double x = *inX++;
1187 double y = *inY++;
1188 double z = inZ ? *inZ++ : 0;
1189
1190 QgsPoint interpolatedPoint( x, y, z ); // general case (not a triangle)
1191
1192 *outX++ = x;
1193 *outY++ = y;
1194 if ( triangle->exteriorRing()->numPoints() == 4 ) // triangle case
1195 {
1196 interpolatedPoint = interpolatePointOnTriangle( triangle, x, y );
1197 }
1198 double tempOutZ = std::isnan( interpolatedPoint.z() ) ? 0.0 : interpolatedPoint.z();
1199 *outZ++ = tempOutZ;
1200
1201 if ( mExtrusionEnabled )
1202 *extZOut++ = tempOutZ + extrusion;
1203
1204 mResults->mRawPoints.append( interpolatedPoint );
1205 mResults->minZ = std::min( mResults->minZ, interpolatedPoint.z() );
1206 mResults->maxZ = std::max( mResults->maxZ, interpolatedPoint.z() );
1207 if ( mExtrusionEnabled )
1208 {
1209 mResults->minZ = std::min( mResults->minZ, interpolatedPoint.z() + extrusion );
1210 mResults->maxZ = std::max( mResults->maxZ, interpolatedPoint.z() + extrusion );
1211 }
1212
1213 const double distance = mProfileCurveEngine->lineLocatePoint( interpolatedPoint, &lastError );
1214 *outDistance++ = distance;
1215
1216 mResults->mDistanceToHeightMap.insert( distance, interpolatedPoint.z() );
1217 lastDistanceAlongProfileCurve = distance;
1218 }
1219
1220 // insert nan point to end the line
1221 mResults->mDistanceToHeightMap.insert( lastDistanceAlongProfileCurve + 0.000001, std::numeric_limits<double>::quiet_NaN() );
1222
1223 if ( mFeedback->isCanceled() )
1224 return;
1225
1226 if ( mExtrusionEnabled )
1227 {
1228 std::unique_ptr< QgsLineString > ring = std::make_unique< QgsLineString >( newX, newY, newZ );
1229 std::unique_ptr< QgsLineString > extrudedRing = std::make_unique< QgsLineString >( newX, newY, extrudedZ );
1230 std::unique_ptr< QgsLineString > reversedExtrusion( extrudedRing->reversed() );
1231 ring->append( reversedExtrusion.get() );
1232 ring->close();
1233 transformedParts.append( QgsGeometry( new QgsPolygon( ring.release() ) ) );
1234
1235 std::unique_ptr< QgsLineString > distanceVHeightRing = std::make_unique< QgsLineString >( newDistance, newZ );
1236 std::unique_ptr< QgsLineString > extrudedDistanceVHeightRing = std::make_unique< QgsLineString >( newDistance, extrudedZ );
1237 std::unique_ptr< QgsLineString > reversedDistanceVHeightExtrusion( extrudedDistanceVHeightRing->reversed() );
1238 distanceVHeightRing->append( reversedDistanceVHeightExtrusion.get() );
1239 distanceVHeightRing->close();
1240 crossSectionParts.append( QgsGeometry( new QgsPolygon( distanceVHeightRing.release() ) ) );
1241 }
1242 else
1243 {
1244 transformedParts.append( QgsGeometry( new QgsLineString( newX, newY, newZ ) ) );
1245 crossSectionParts.append( QgsGeometry( new QgsLineString( newDistance, newZ ) ) );
1246 }
1247};
1248
1249void QgsVectorLayerProfileGenerator::processTriangleIntersectForPolygon( const QgsPolygon *sourcePolygon, const QgsPolygon *intersectionPolygon, QVector< QgsGeometry > &transformedParts, QVector< QgsGeometry > &crossSectionParts )
1250{
1251 bool oldExtrusion = mExtrusionEnabled;
1252
1253 /* Polyone extrusion produces I or C or inverted C shapes because the starting and ending points are the same.
1254 We observe the same case with linestrings if the starting and ending points are not at the ends.
1255 In the case below, the Z polygon projected onto the curve produces a shape that cannot be used to represent the extrusion ==> we would obtain a 3D volume.
1256 In order to avoid having strange shapes that cannot be understood by the end user, extrusion is deactivated in the case of polygons.
1257
1258 .^..
1259 ./ | \..
1260 ../ | \...
1261 ../ | \...
1262 ../ | \.. ....^..
1263 ../ | ........\.../ \... ^
1264 ../ ......|......./ \... \.... .../ \
1265 /,........../ | \.. \... / \
1266 v | \... ..../ \... \
1267 | \ ./ \... \
1268 | v \.. \
1269 | `v
1270 |
1271 .^..
1272 ./ \..
1273 ../ \...
1274 ../ \...
1275 ../ \.. ....^..
1276 ../ ........\.../ \... ^
1277 ../ ............../ \... \.... .../ \
1278 /,........../ \.. \... / \
1279 v \... ..../ \... \
1280 \ ./ \... \
1281 v \.. \
1282 `v
1283 */
1284 mExtrusionEnabled = false;
1285 if ( mProfileBufferedCurveEngine->contains( sourcePolygon ) ) // sourcePolygon is entirely inside curve buffer, we keep it as whole
1286 {
1287 if ( const QgsCurve *exterior = sourcePolygon->exteriorRing() )
1288 {
1289 QgsLineString *exteriorLine = qgsgeometry_cast<QgsLineString *>( exterior );
1290 processTriangleIntersectForLine( sourcePolygon, exteriorLine, transformedParts, crossSectionParts );
1291 }
1292 for ( int i = 0; i < sourcePolygon->numInteriorRings(); ++i )
1293 {
1294 QgsLineString *interiorLine = qgsgeometry_cast<QgsLineString *>( sourcePolygon->interiorRing( i ) );
1295 processTriangleIntersectForLine( sourcePolygon, interiorLine, transformedParts, crossSectionParts );
1296 }
1297 }
1298 else // sourcePolygon is partially inside curve buffer, the intersectionPolygon is closed due to the intersection operation then
1299 // it must be 'reopened'
1300 {
1301 if ( const QgsCurve *exterior = intersectionPolygon->exteriorRing() )
1302 {
1303 QgsLineString *exteriorLine = qgsgeometry_cast<QgsLineString *>( exterior )->clone();
1304 exteriorLine->deleteVertex( QgsVertexId( 0, 0, exteriorLine->numPoints() - 1 ) ); // open linestring
1305 processTriangleIntersectForLine( sourcePolygon, exteriorLine, transformedParts, crossSectionParts );
1306 delete exteriorLine;
1307 }
1308 for ( int i = 0; i < intersectionPolygon->numInteriorRings(); ++i )
1309 {
1310 QgsLineString *interiorLine = qgsgeometry_cast<QgsLineString *>( intersectionPolygon->interiorRing( i ) );
1311 if ( mProfileBufferedCurveEngine->contains( interiorLine ) ) // interiorLine is entirely inside curve buffer
1312 {
1313 processTriangleIntersectForLine( sourcePolygon, interiorLine, transformedParts, crossSectionParts );
1314 }
1315 else
1316 {
1317 interiorLine = qgsgeometry_cast<QgsLineString *>( intersectionPolygon->interiorRing( i ) )->clone();
1318 interiorLine->deleteVertex( QgsVertexId( 0, 0, interiorLine->numPoints() - 1 ) ); // open linestring
1319 processTriangleIntersectForLine( sourcePolygon, interiorLine, transformedParts, crossSectionParts );
1320 delete interiorLine;
1321 }
1322 }
1323 }
1324
1325 mExtrusionEnabled = oldExtrusion;
1326};
1327
1328bool QgsVectorLayerProfileGenerator::generateProfileForPolygons()
1329{
1330 // get features from layer
1331 QgsFeatureRequest request;
1332 request.setDestinationCrs( mTargetCrs, mTransformContext );
1333 if ( mTolerance > 0 )
1334 {
1335 request.setDistanceWithin( QgsGeometry( mProfileCurve->clone() ), mTolerance );
1336 }
1337 else
1338 {
1339 request.setFilterRect( mProfileCurve->boundingBox() );
1340 }
1341 request.setSubsetOfAttributes( mDataDefinedProperties.referencedFields( mExpressionContext ), mFields );
1342 request.setFeedback( mFeedback.get() );
1343
1344 std::function< void( const QgsPolygon *triangle, const QgsAbstractGeometry *intersect, QVector< QgsGeometry > &, QVector< QgsGeometry > & ) > processTriangleLineIntersect;
1345 processTriangleLineIntersect = [this]( const QgsPolygon * triangle, const QgsAbstractGeometry * intersection, QVector< QgsGeometry > &transformedParts, QVector< QgsGeometry > &crossSectionParts )
1346 {
1347 for ( auto it = intersection->const_parts_begin();
1348 ! mFeedback->isCanceled() && it != intersection->const_parts_end();
1349 ++it )
1350 {
1351 // intersect may be a (multi)point or (multi)linestring
1352 switch ( QgsWkbTypes::geometryType( ( *it )->wkbType() ) )
1353 {
1355 if ( const QgsPoint *p = qgsgeometry_cast< const QgsPoint * >( *it ) )
1356 {
1357 processTriangleIntersectForPoint( triangle, p, transformedParts, crossSectionParts );
1358 }
1359 break;
1360
1362 if ( const QgsLineString *intersectionLine = qgsgeometry_cast< const QgsLineString * >( *it ) )
1363 {
1364 processTriangleIntersectForLine( triangle, intersectionLine, transformedParts, crossSectionParts );
1365 }
1366 break;
1367
1369 if ( const QgsPolygon *poly = qgsgeometry_cast< const QgsPolygon * >( *it ) )
1370 {
1371 processTriangleIntersectForPolygon( triangle, poly, transformedParts, crossSectionParts );
1372 }
1373 break;
1374
1377 return;
1378 }
1379 }
1380 };
1381
1382 auto triangleIsCollinearInXYPlane = []( const QgsPolygon * polygon )-> bool
1383 {
1384 const QgsLineString *ring = qgsgeometry_cast< const QgsLineString * >( polygon->exteriorRing() );
1385 return QgsGeometryUtilsBase::pointsAreCollinear( ring->xAt( 0 ), ring->yAt( 0 ),
1386 ring->xAt( 1 ), ring->yAt( 1 ),
1387 ring->xAt( 2 ), ring->yAt( 2 ), 0.005 );
1388 };
1389
1390 auto processPolygon = [this, &processTriangleLineIntersect, &triangleIsCollinearInXYPlane]( const QgsCurvePolygon * polygon, QVector< QgsGeometry > &transformedParts, QVector< QgsGeometry > &crossSectionParts, double offset, bool & wasCollinear )
1391 {
1392 std::unique_ptr< QgsPolygon > clampedPolygon;
1393 if ( const QgsPolygon *p = qgsgeometry_cast< const QgsPolygon * >( polygon ) )
1394 {
1395 clampedPolygon.reset( p->clone() );
1396 }
1397 else
1398 {
1399 clampedPolygon.reset( qgsgeometry_cast< QgsPolygon * >( polygon->segmentize() ) );
1400 }
1401 clampAltitudes( clampedPolygon.get(), offset );
1402
1403 if ( mFeedback->isCanceled() )
1404 return;
1405
1406 if ( mTolerance > 0.0 ) // if the tolerance is not 0.0 we will have a polygon / polygon intersection, we do not need tessellation
1407 {
1408 QString error;
1409 if ( mProfileBufferedCurveEngine->intersects( clampedPolygon.get(), &error ) )
1410 {
1411 std::unique_ptr< QgsAbstractGeometry > intersection;
1412 intersection.reset( mProfileBufferedCurveEngine->intersection( clampedPolygon.get(), &error ) );
1413 if ( error.isEmpty() )
1414 {
1415 processTriangleLineIntersect( clampedPolygon.get(), intersection.get(), transformedParts, crossSectionParts );
1416 }
1417 else
1418 {
1419 // this case may occur with vertical object as geos does not handle very well 3D data.
1420 // Geos works in 2D from the 3D coordinates then re-add the Z values, but when 2D-from-3D objects are vertical, they are topologically incorrects!
1421 // This piece of code is just a fix to handle this case, a better and real 3D capable library is needed (like SFCGAL).
1422 QgsLineString *ring = qgsgeometry_cast< QgsLineString * >( clampedPolygon->exteriorRing() );
1423 int numPoints = ring->numPoints();
1424 QVector< double > newX( numPoints );
1425 QVector< double > newY( numPoints );
1426 QVector< double > newZ( numPoints );
1427 double *outX = newX.data();
1428 double *outY = newY.data();
1429 double *outZ = newZ.data();
1430
1431 const double *inX = ring->xData();
1432 const double *inY = ring->yData();
1433 const double *inZ = ring->zData();
1434 for ( int i = 0 ; ! mFeedback->isCanceled() && i < ring->numPoints() - 1; ++i )
1435 {
1436 *outX++ = inX[i] + i * 1.0e-9;
1437 *outY++ = inY[i] + i * 1.0e-9;
1438 *outZ++ = inZ[i];
1439 }
1440 std::unique_ptr< QgsPolygon > shiftedPoly;
1441 shiftedPoly.reset( new QgsPolygon( new QgsLineString( newX, newY, newZ ) ) );
1442
1443 intersection.reset( mProfileBufferedCurveEngine->intersection( shiftedPoly.get(), &error ) );
1444 if ( intersection.get() )
1445 processTriangleLineIntersect( clampedPolygon.get(), intersection.get(), transformedParts, crossSectionParts );
1446 else
1447 QgsDebugMsgLevel( QStringLiteral( "processPolygon after shift bad geom! error: %1" ).arg( error ), 0 );
1448 }
1449 }
1450
1451 }
1452 else // ie. polygon / line intersection ==> need tessellation
1453 {
1454 QgsGeometry tessellation;
1455 if ( clampedPolygon->numInteriorRings() == 0 && clampedPolygon->exteriorRing() && clampedPolygon->exteriorRing()->numPoints() == 4 && clampedPolygon->exteriorRing()->isClosed() )
1456 {
1457 // special case -- polygon is already a triangle, so no need to tessellate
1458 std::unique_ptr< QgsMultiPolygon > multiPolygon = std::make_unique< QgsMultiPolygon >();
1459 multiPolygon->addGeometry( clampedPolygon.release() );
1460 tessellation = QgsGeometry( std::move( multiPolygon ) );
1461 }
1462 else
1463 {
1464 const QgsRectangle bounds = clampedPolygon->boundingBox();
1465 QgsTessellator t( bounds, false, false, false, false );
1466 t.addPolygon( *clampedPolygon, 0 );
1467
1468 tessellation = QgsGeometry( t.asMultiPolygon() );
1469 if ( mFeedback->isCanceled() )
1470 return;
1471
1472 tessellation.translate( bounds.xMinimum(), bounds.yMinimum() );
1473 }
1474
1475 // iterate through the tessellation, finding triangles which intersect the line
1476 const int numTriangles = qgsgeometry_cast< const QgsMultiPolygon * >( tessellation.constGet() )->numGeometries();
1477 for ( int i = 0; ! mFeedback->isCanceled() && i < numTriangles; ++i )
1478 {
1479 const QgsPolygon *triangle = qgsgeometry_cast< const QgsPolygon * >( qgsgeometry_cast< const QgsMultiPolygon * >( tessellation.constGet() )->geometryN( i ) );
1480
1481 if ( triangleIsCollinearInXYPlane( triangle ) )
1482 {
1483 wasCollinear = true;
1484 const QgsLineString *ring = qgsgeometry_cast< const QgsLineString * >( polygon->exteriorRing() );
1485
1486 QString lastError;
1487 if ( const QgsLineString *ls = qgsgeometry_cast< const QgsLineString * >( mProfileCurve.get() ) )
1488 {
1489 for ( int curveSegmentIndex = 0; curveSegmentIndex < mProfileCurve->numPoints() - 1; ++curveSegmentIndex )
1490 {
1491 const QgsPoint p1 = ls->pointN( curveSegmentIndex );
1492 const QgsPoint p2 = ls->pointN( curveSegmentIndex + 1 );
1493
1494 QgsPoint intersectionPoint;
1495 double minZ = std::numeric_limits< double >::max();
1496 double maxZ = std::numeric_limits< double >::lowest();
1497
1498 for ( auto vertexPair : std::array<std::pair<int, int>, 3> {{ { 0, 1}, {1, 2}, {2, 0} }} )
1499 {
1500 bool isIntersection = false;
1501 if ( QgsGeometryUtils::segmentIntersection( ring->pointN( vertexPair.first ), ring->pointN( vertexPair.second ), p1, p2, intersectionPoint, isIntersection ) )
1502 {
1503 const double fraction = QgsGeometryUtilsBase::pointFractionAlongLine( ring->xAt( vertexPair.first ), ring->yAt( vertexPair.first ), ring->xAt( vertexPair.second ), ring->yAt( vertexPair.second ), intersectionPoint.x(), intersectionPoint.y() );
1504 const double intersectionZ = ring->zAt( vertexPair.first ) + ( ring->zAt( vertexPair.second ) - ring->zAt( vertexPair.first ) ) * fraction;
1505 minZ = std::min( minZ, intersectionZ );
1506 maxZ = std::max( maxZ, intersectionZ );
1507 }
1508 }
1509
1510 if ( !intersectionPoint.isEmpty() )
1511 {
1512 // need z?
1513 mResults->mRawPoints.append( intersectionPoint );
1514 mResults->minZ = std::min( mResults->minZ, minZ );
1515 mResults->maxZ = std::max( mResults->maxZ, maxZ );
1516
1517 const double distance = mProfileCurveEngine->lineLocatePoint( intersectionPoint, &lastError );
1518
1519 crossSectionParts.append( QgsGeometry( new QgsLineString( QVector< double > {distance, distance}, QVector< double > {minZ, maxZ} ) ) );
1520
1521 mResults->mDistanceToHeightMap.insert( distance, minZ );
1522 mResults->mDistanceToHeightMap.insert( distance, maxZ );
1523 }
1524 }
1525 }
1526 else
1527 {
1528 // curved geometries, not supported yet, but not possible through the GUI anyway
1529 QgsDebugError( QStringLiteral( "Collinear triangles with curved profile lines are not supported yet" ) );
1530 }
1531 }
1532 else // not collinear
1533 {
1534 QString error;
1535 if ( mProfileBufferedCurveEngine->intersects( triangle, &error ) )
1536 {
1537 std::unique_ptr< QgsAbstractGeometry > intersection( mProfileBufferedCurveEngine->intersection( triangle, &error ) );
1538 processTriangleLineIntersect( triangle, intersection.get(), transformedParts, crossSectionParts );
1539 }
1540 }
1541 }
1542 }
1543 };
1544
1545 // ========= MAIN JOB
1546 QgsFeature feature;
1547 QgsFeatureIterator it = mSource->getFeatures( request );
1548 while ( ! mFeedback->isCanceled() && it.nextFeature( feature ) )
1549 {
1550 if ( !mProfileBufferedCurveEngine->intersects( feature.geometry().constGet() ) )
1551 continue;
1552
1553 mExpressionContext.setFeature( feature );
1554
1555 const double offset = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::Property::ZOffset, mExpressionContext, mOffset );
1556 const QgsGeometry g = feature.geometry();
1557 QVector< QgsGeometry > transformedParts;
1558 QVector< QgsGeometry > crossSectionParts;
1559 bool wasCollinear = false;
1560
1561 // === process intersection of geometry feature parts with the mProfileBoxEngine
1562 for ( auto it = g.const_parts_begin(); ! mFeedback->isCanceled() && it != g.const_parts_end(); ++it )
1563 {
1564 if ( mProfileBufferedCurveEngine->intersects( *it ) )
1565 {
1566 if ( const QgsCurvePolygon *curvePolygon = qgsgeometry_cast< const QgsCurvePolygon * >( *it ) )
1567 {
1568 processPolygon( curvePolygon, transformedParts, crossSectionParts, offset, wasCollinear );
1569 }
1570 else if ( const QgsPolyhedralSurface *polySurface = qgsgeometry_cast< const QgsPolyhedralSurface * >( *it ) )
1571 {
1572 for ( int i = 0; i < polySurface->numPatches(); ++i )
1573 {
1574 const QgsPolygon *polygon = polySurface->patchN( i );
1575 if ( mProfileBufferedCurveEngine->intersects( polygon ) )
1576 {
1577 processPolygon( polygon, transformedParts, crossSectionParts, offset, wasCollinear );
1578 }
1579 }
1580 }
1581 else
1582 {
1583 QgsDebugError( QStringLiteral( "Unhandled Geometry type: %1" ).arg( ( *it )->wktTypeStr() ) );
1584 }
1585 }
1586 }
1587
1588 if ( mFeedback->isCanceled() )
1589 return false;
1590
1591 // === aggregate results for this feature
1593 resultFeature.featureId = feature.id();
1594 resultFeature.geometry = transformedParts.size() > 1 ? QgsGeometry::collectGeometry( transformedParts ) : transformedParts.value( 0 );
1595 if ( !crossSectionParts.empty() )
1596 {
1597 if ( !wasCollinear )
1598 {
1599 QgsGeometry unioned = QgsGeometry::unaryUnion( crossSectionParts );
1600 if ( unioned.isEmpty() )
1601 {
1602 resultFeature.crossSectionGeometry = QgsGeometry::collectGeometry( crossSectionParts );
1603 }
1604 else
1605 {
1606 if ( unioned.type() == Qgis::GeometryType::Line )
1607 {
1608 unioned = unioned.mergeLines();
1609 }
1610 resultFeature.crossSectionGeometry = unioned;
1611 }
1612 }
1613 else
1614 {
1615 resultFeature.crossSectionGeometry = QgsGeometry::collectGeometry( crossSectionParts );
1616 }
1617 }
1618 mResults->features[resultFeature.featureId].append( resultFeature );
1619 }
1620 return true;
1621}
1622
1623double QgsVectorLayerProfileGenerator::terrainHeight( double x, double y )
1624{
1625 if ( !mTerrainProvider )
1626 return std::numeric_limits<double>::quiet_NaN();
1627
1628 // transform feature point to terrain provider crs
1629 try
1630 {
1631 double dummyZ = 0;
1632 mTargetToTerrainProviderTransform.transformInPlace( x, y, dummyZ );
1633 }
1634 catch ( QgsCsException & )
1635 {
1636 return std::numeric_limits<double>::quiet_NaN();
1637 }
1638
1639 return mTerrainProvider->heightAt( x, y );
1640}
1641
1642double QgsVectorLayerProfileGenerator::featureZToHeight( double x, double y, double z, double offset )
1643{
1644 switch ( mClamping )
1645 {
1647 break;
1648
1651 {
1652 const double terrainZ = terrainHeight( x, y );
1653 if ( !std::isnan( terrainZ ) )
1654 {
1655 switch ( mClamping )
1656 {
1658 if ( std::isnan( z ) )
1659 z = terrainZ;
1660 else
1661 z += terrainZ;
1662 break;
1663
1665 z = terrainZ;
1666 break;
1667
1669 break;
1670 }
1671 }
1672 break;
1673 }
1674 }
1675
1676 return ( std::isnan( z ) ? 0 : z ) * mScale + offset;
1677}
1678
1679void QgsVectorLayerProfileGenerator::clampAltitudes( QgsLineString *lineString, const QgsPoint &centroid, double offset )
1680{
1681 for ( int i = 0; i < lineString->nCoordinates(); ++i )
1682 {
1683 if ( mFeedback->isCanceled() )
1684 break;
1685
1686 double terrainZ = 0;
1687 switch ( mClamping )
1688 {
1691 {
1692 QgsPointXY pt;
1693 switch ( mBinding )
1694 {
1696 pt.setX( lineString->xAt( i ) );
1697 pt.setY( lineString->yAt( i ) );
1698 break;
1699
1701 pt.set( centroid.x(), centroid.y() );
1702 break;
1703 }
1704
1705 terrainZ = terrainHeight( pt.x(), pt.y() );
1706 break;
1707 }
1708
1710 break;
1711 }
1712
1713 double geomZ = 0;
1714
1715 switch ( mClamping )
1716 {
1719 geomZ = lineString->zAt( i );
1720 break;
1721
1723 break;
1724 }
1725
1726 const double z = ( terrainZ + ( std::isnan( geomZ ) ? 0 : geomZ ) ) * mScale + offset;
1727 lineString->setZAt( i, z );
1728 }
1729}
1730
1731bool QgsVectorLayerProfileGenerator::clampAltitudes( QgsPolygon *polygon, double offset )
1732{
1733 if ( !polygon->is3D() )
1734 polygon->addZValue( 0 );
1735
1736 QgsPoint centroid;
1737 switch ( mBinding )
1738 {
1740 break;
1741
1743 centroid = polygon->centroid();
1744 break;
1745 }
1746
1747 QgsCurve *curve = const_cast<QgsCurve *>( polygon->exteriorRing() );
1748 QgsLineString *lineString = qgsgeometry_cast<QgsLineString *>( curve );
1749 if ( !lineString )
1750 return false;
1751
1752 clampAltitudes( lineString, centroid, offset );
1753
1754 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
1755 {
1756 if ( mFeedback->isCanceled() )
1757 break;
1758
1759 QgsCurve *curve = const_cast<QgsCurve *>( polygon->interiorRing( i ) );
1760 QgsLineString *lineString = qgsgeometry_cast<QgsLineString *>( curve );
1761 if ( !lineString )
1762 return false;
1763
1764 clampAltitudes( lineString, centroid, offset );
1765 }
1766 return true;
1767}
The Qgis class provides global constants for use throughout the application.
Definition qgis.h:54
@ Relative
Elevation is relative to terrain height (final elevation = terrain elevation + feature elevation)
@ Terrain
Elevation is clamped to terrain (final elevation = terrain elevation)
@ Absolute
Elevation is taken directly from feature and is independent of terrain height (final elevation = feat...
VertexType
Types of vertex.
Definition qgis.h:2835
@ Polygon
Polygons.
@ Unknown
Unknown types.
@ Null
No geometry.
@ Round
Use rounded joins.
@ Centroid
Clamp just centroid of feature.
@ Vertex
Clamp every vertex of feature.
@ Flat
Flat cap (in line with start/end of line)
@ ContinuousSurface
The features should be treated as representing values on a continuous surface (eg contour lines)
@ IndividualFeatures
Treat each feature as an individual object (eg buildings)
ProfileExportType
Types of export for elevation profiles.
Definition qgis.h:3929
@ Profile2D
Export profiles as 2D profile lines, with elevation stored in exported geometry Y dimension and dista...
@ Features3D
Export profiles as 3D features, with elevation values stored in exported geometry Z values.
@ DistanceVsElevationTable
Export profiles as a table of sampled distance vs elevation values.
@ Marker
Marker symbol.
@ Reverse
Reverse/inverse transform (from destination to source)
Abstract base class for all geometries.
bool is3D() const
Returns true if the geometry is 3D and contains a z-value.
virtual QgsPoint centroid() const
Returns the centroid of the geometry.
Abstract base class for objects which generate elevation profiles.
Abstract base class for storage of elevation profiles.
Abstract base class for objects which generate elevation profiles which represent a continuous surfac...
void renderResults(QgsProfileRenderContext &context) override
Renders the results to the specified context.
void copyPropertiesFromGenerator(const QgsAbstractProfileGenerator *generator) override
Copies properties from specified generator to the results object.
QVector< QgsAbstractProfileResults::Feature > asFeatures(Qgis::ProfileExportType type, QgsFeedback *feedback=nullptr) const override
Returns a list of features representing the calculated elevation results.
QVector< QgsProfileIdentifyResults > identify(const QgsProfilePoint &point, const QgsProfileIdentifyContext &context) override
Identify results visible at the specified profile point.
QgsProfileSnapResult snapPoint(const QgsProfilePoint &point, const QgsProfileSnapContext &context) override
Snaps a point to the generated elevation profile.
double valueAsDouble(int key, const QgsExpressionContext &context, double defaultValue=0.0, bool *ok=nullptr) const
Calculates the current value of the property with the specified key and interprets it as a double.
Class for doing transforms between two map coordinate systems.
void transformInPlace(double &x, double &y, double &z, Qgis::TransformDirection direction=Qgis::TransformDirection::Forward) const
Transforms an array of x, y and z double coordinates in place, from the source CRS to the destination...
Custom exception class for Coordinate Reference System related exceptions.
Curve polygon geometry type.
int numInteriorRings() const
Returns the number of interior rings contained with the curve polygon.
bool addZValue(double zValue=0) override
Adds a z-dimension to the geometry, initialized to a preset value.
const QgsCurve * exteriorRing() const
Returns the curve polygon's exterior ring.
const QgsCurve * interiorRing(int i) const
Retrieves an interior ring from the curve polygon.
Abstract base class for curved geometry type.
Definition qgscurve.h:35
virtual int numPoints() const =0
Returns the number of points in the curve.
virtual bool pointAt(int node, QgsPoint &point, Qgis::VertexType &type) const =0
Returns the point and vertex id of a point within the curve.
QgsRange which stores a range of double values.
Definition qgsrange.h:231
RAII class to pop scope from an expression context on destruction.
void setFeature(const QgsFeature &feature)
Convenience function for setting a feature for the context.
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).
QgsFeatureRequest & setCoordinateTransform(const QgsCoordinateTransform &transform)
Sets the coordinate transform which will be used to transform the feature's geometries.
QgsFeatureRequest & setFilterFids(const QgsFeatureIds &fids)
Sets the feature IDs that should be fetched.
QgsFeatureRequest & setSubsetOfAttributes(const QgsAttributeList &attrs)
Set a subset of attributes that will be fetched.
QgsFeatureRequest & setDestinationCrs(const QgsCoordinateReferenceSystem &crs, const QgsCoordinateTransformContext &context)
Sets the destination crs for feature's geometries.
void setFeedback(QgsFeedback *feedback)
Attach a feedback object that can be queried regularly by the iterator to check if it should be cance...
QgsFeatureRequest & setFilterRect(const QgsRectangle &rectangle)
Sets the rectangle from which features will be taken.
QgsFeatureRequest & setDistanceWithin(const QgsGeometry &geometry, double distance)
Sets a reference geometry and a maximum distance from this geometry to retrieve features within.
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition qgsfeature.h:58
QgsFeatureId id
Definition qgsfeature.h:66
QgsGeometry geometry
Definition qgsfeature.h:69
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
A fill symbol type, for rendering Polygon and MultiPolygon geometries.
static double pointFractionAlongLine(double x1, double y1, double x2, double y2, double px, double py)
Given the line (x1, y1) to (x2, y2) and a point (px, py) returns the fraction of the line length at w...
static bool pointsAreCollinear(double x1, double y1, double x2, double y2, double x3, double y3, double epsilon)
Given the points (x1, y1), (x2, y2) and (x3, y3) returns true if these points can be considered colli...
static bool segmentIntersection(const QgsPoint &p1, const QgsPoint &p2, const QgsPoint &q1, const QgsPoint &q2, QgsPoint &intersectionPoint, bool &isIntersection, double tolerance=1e-8, bool acceptImproperIntersection=false)
Compute the intersection between two segments.
A geometry is the spatial representation of a feature.
static QgsGeometry fromRect(const QgsRectangle &rect)
Creates a new geometry from a QgsRectangle.
QgsAbstractGeometry::const_part_iterator const_parts_begin() const
Returns STL-style const iterator pointing to the first part of the geometry.
bool boundingBoxIntersects(const QgsRectangle &rectangle) const
Returns true if the bounding box of this geometry intersects with a rectangle.
static QgsGeometry collectGeometry(const QVector< QgsGeometry > &geometries)
Creates a new multipart geometry from a list of QgsGeometry objects.
Qgis::GeometryOperationResult transform(const QgsCoordinateTransform &ct, Qgis::TransformDirection direction=Qgis::TransformDirection::Forward, bool transformZ=false)
Transforms this geometry as described by the coordinate transform ct.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
static QgsGeometry fromPointXY(const QgsPointXY &point)
Creates a new geometry from a QgsPointXY object.
Qgis::GeometryType type
QgsGeometry mergeLines() const
Merges any connected lines in a LineString/MultiLineString geometry and converts them to single line ...
bool isEmpty() const
Returns true if the geometry is empty (eg a linestring with no vertices, or a collection with no geom...
QgsAbstractGeometry::const_part_iterator const_parts_end() const
Returns STL-style iterator pointing to the imaginary part after the last part of the geometry.
static QgsGeometry unaryUnion(const QVector< QgsGeometry > &geometries, const QgsGeometryParameters &parameters=QgsGeometryParameters())
Compute the unary union on a list of geometries.
Qgis::GeometryOperationResult translate(double dx, double dy, double dz=0.0, double dm=0.0)
Translates this geometry by dx, dy, dz and dm.
Does vector analysis using the geos library and handles import, export, exception handling*.
Definition qgsgeos.h:137
Line string geometry type, with support for z-dimension and m-values.
const double * yData() const
Returns a const pointer to the y vertex data.
const double * xData() const
Returns a const pointer to the x vertex data.
const double * zData() const
Returns a const pointer to the z vertex data, or nullptr if the linestring does not have z values.
int numPoints() const override
Returns the number of points in the curve.
QgsPoint pointN(int i) const
Returns the specified point from inside the line string.
int nCoordinates() const override
Returns the number of nodes contained in the geometry.
double yAt(int index) const override
Returns the y-coordinate of the specified node in the line string.
void setZAt(int index, double z)
Sets the z-coordinate of the specified node in the line string.
bool deleteVertex(QgsVertexId position) override
Deletes a vertex within the geometry.
QgsLineString * clone() const override
Clones the geometry by performing a deep copy.
double zAt(int index) const override
Returns the z-coordinate of the specified node in the line string.
double xAt(int index) const override
Returns the x-coordinate of the specified node in the line string.
A line symbol type, for rendering LineString and MultiLineString geometries.
A marker symbol type, for rendering Point and MultiPoint geometries.
Multi line string geometry collection.
Multi point geometry collection.
Multi polygon geometry collection.
A class to represent a 2D point.
Definition qgspointxy.h:60
void setY(double y)
Sets the y value of the point.
Definition qgspointxy.h:129
void set(double x, double y)
Sets the x and y value of the point.
Definition qgspointxy.h:136
double y
Definition qgspointxy.h:64
double x
Definition qgspointxy.h:63
void setX(double x)
Sets the x value of the point.
Definition qgspointxy.h:119
Point geometry type, with support for z-dimension and m-values.
Definition qgspoint.h:49
QgsPoint * clone() const override
Clones the geometry by performing a deep copy.
Definition qgspoint.cpp:104
double z
Definition qgspoint.h:54
double x
Definition qgspoint.h:52
bool isEmpty() const override
Returns true if the geometry is empty.
Definition qgspoint.cpp:737
double y
Definition qgspoint.h:53
Polygon geometry type.
Definition qgspolygon.h:33
Polyhedral surface geometry type.
Encapsulates the context in which an elevation profile is to be generated.
Encapsulates the context of identifying profile results.
double maximumPointElevationDelta
Maximum allowed snapping delta for the elevation values when identifying a point.
double maximumPointDistanceDelta
Maximum allowed snapping delta for the distance values when identifying a point.
double maximumSurfaceElevationDelta
Maximum allowed snapping delta for the elevation values when identifying a continuous elevation surfa...
Stores identify results generated by a QgsAbstractProfileResults object.
Encapsulates a point on a distance-elevation profile.
double elevation() const
Returns the elevation of the point.
double distance() const
Returns the distance of the point.
Abstract base class for storage of elevation profiles.
const QTransform & worldTransform() const
Returns the transform from world coordinates to painter coordinates.
QgsDoubleRange elevationRange() const
Returns the range of elevations to include in the render.
QgsDoubleRange distanceRange() const
Returns the range of distances to include in the render.
QgsRenderContext & renderContext()
Returns a reference to the component QgsRenderContext.
Encapsulates properties and constraints relating to fetching elevation profiles from different source...
Encapsulates the context of snapping a profile point.
double maximumPointDistanceDelta
Maximum allowed snapping delta for the distance values when snapping to a point.
double maximumSurfaceElevationDelta
Maximum allowed snapping delta for the elevation values when snapping to a continuous elevation surfa...
double maximumPointElevationDelta
Maximum allowed snapping delta for the elevation values when snapping to a point.
Encapsulates results of snapping a profile point.
QgsProfilePoint snappedPoint
Snapped point.
bool prepare(const QgsExpressionContext &context=QgsExpressionContext()) const final
Prepares the collection against a specified expression context.
QSet< QString > referencedFields(const QgsExpressionContext &context=QgsExpressionContext(), bool ignoreContext=false) const final
Returns the set of any fields referenced by the active properties from the collection.
T lower() const
Returns the lower bound of the range.
Definition qgsrange.h:78
T upper() const
Returns the upper bound of the range.
Definition qgsrange.h:85
A rectangle specified with double values.
double xMinimum() const
Returns the x minimum value (left side of rectangle).
bool intersects(const QgsRectangle &rect) const
Returns true when rectangle intersects with other rectangle.
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
double xMaximum() const
Returns the x maximum value (right side of rectangle).
QPainter * painter()
Returns the destination QPainter for the render operation.
QgsExpressionContext & expressionContext()
Gets the expression context.
Scoped object for saving and restoring a QPainter object's state.
Abstract base class for all rendered symbols.
Definition qgssymbol.h:231
void setColor(const QColor &color) const
Sets the color for the symbol.
qreal opacity() const
Returns the opacity for the symbol.
Definition qgssymbol.h:632
QColor color() const
Returns the symbol's color.
Qgis::SymbolType type() const
Returns the symbol's type.
Definition qgssymbol.h:293
Class that takes care of tessellation of polygons into triangles.
Vector layer specific subclass of QgsMapLayerElevationProperties.
Partial snapshot of vector layer's state (only the members necessary for access to features)
Implementation of QgsAbstractProfileGenerator for vector layers.
QgsAbstractProfileResults * takeResults() override
Takes results from the generator.
bool generateProfile(const QgsProfileGenerationContext &context=QgsProfileGenerationContext()) override
Generate the profile (based on data stored in the class).
QString sourceId() const override
Returns a unique identifier representing the source of the profile.
~QgsVectorLayerProfileGenerator() override
QgsFeedback * feedback() const override
Access to feedback object of the generator (may be nullptr)
QgsVectorLayerProfileGenerator(QgsVectorLayer *layer, const QgsProfileRequest &request)
Constructor for QgsVectorLayerProfileGenerator.
QVector< QgsProfileIdentifyResults > identify(const QgsProfilePoint &point, const QgsProfileIdentifyContext &context) override
Identify results visible at the specified profile point.
std::unique_ptr< QgsMarkerSymbol > mMarkerSymbol
QVector< QgsGeometry > asGeometries() const override
Returns a list of geometries representing the calculated elevation results.
void renderResults(QgsProfileRenderContext &context) override
Renders the results to the specified context.
QVector< QgsAbstractProfileResults::Feature > asFeatures(Qgis::ProfileExportType type, QgsFeedback *feedback=nullptr) const override
Returns a list of features representing the calculated elevation results.
QgsProfileSnapResult snapPoint(const QgsProfilePoint &point, const QgsProfileSnapContext &context) override
Snaps a point to the generated elevation profile.
QString type() const override
Returns the unique string identifier for the results type.
QHash< QgsFeatureId, QVector< Feature > > features
void copyPropertiesFromGenerator(const QgsAbstractProfileGenerator *generator) override
Copies properties from specified generator to the results object.
Represents a vector layer which manages a vector based data sets.
QgsMapLayerElevationProperties * elevationProperties() override
Returns the layer's elevation properties.
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...
#define BUILTIN_UNREACHABLE
Definition qgis.h:6571
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition qgis.h:5917
QSet< QgsFeatureId > QgsFeatureIds
qint64 QgsFeatureId
64 bit feature ids negative numbers are used for uncommitted/newly added features
#define QgsDebugMsgLevel(str, level)
Definition qgslogger.h:39
#define QgsDebugError(str)
Definition qgslogger.h:38
const QgsCoordinateReferenceSystem & crs
Encapsulates information about a feature exported from the profile results.
QString layerIdentifier
Identifier for grouping output features.
QVariantMap attributes
Exported attributes.
QgsGeometry crossSectionGeometry
Cross section distance vs height geometry for feature.
QgsGeometry geometry
Feature's geometry with any terrain height adjustment and extrusion applied.
Utility class for identifying a unique vertex within a geometry.
Definition qgsvertexid.h:30