QGIS API Documentation 3.30.0-'s-Hertogenbosch (f186b8efe0)
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 "qgsprofilerequest.h"
19#include "qgscurve.h"
20#include "qgsvectorlayer.h"
23#include "qgsgeos.h"
25#include "qgsterrainprovider.h"
26#include "qgspolygon.h"
27#include "qgstessellator.h"
28#include "qgsmultipolygon.h"
29#include "qgsmeshlayerutils.h"
30#include "qgsmultipoint.h"
31#include "qgsmultilinestring.h"
32#include "qgslinesymbol.h"
33#include "qgsfillsymbol.h"
34#include "qgsmarkersymbol.h"
35#include "qgsprofilepoint.h"
36#include "qgsprofilesnapping.h"
38#include <QPolygonF>
39
40//
41// QgsVectorLayerProfileResults
42//
43
45{
46 return QStringLiteral( "vector" );
47}
48
50{
51 QVector<QgsGeometry> res;
52 res.reserve( features.size() );
53 for ( auto it = features.constBegin(); it != features.constEnd(); ++it )
54 {
55 for ( const Feature &feature : it.value() )
56 {
57 res.append( feature.geometry );
58 }
59 }
60 return res;
61}
62
64{
65 switch ( profileType )
66 {
68 return snapPointToIndividualFeatures( point, context );
70 return QgsAbstractProfileSurfaceResults::snapPoint( point, context );
71 }
73}
74
75
76QVector<QgsProfileIdentifyResults> QgsVectorLayerProfileResults::identify( const QgsDoubleRange &distanceRange, const QgsDoubleRange &elevationRange, const QgsProfileIdentifyContext & )
77{
78 QgsFeatureIds ids;
79 auto visitFeature = [&ids]( QgsFeatureId featureId )
80 {
81 ids << featureId;
82 };
83
84 visitFeaturesInRange( distanceRange, elevationRange, visitFeature );
85 if ( ids.empty() )
86 return {};
87
88 QVector< QVariantMap> idsList;
89 for ( auto it = ids.constBegin(); it != ids.constEnd(); ++it )
90 idsList.append( QVariantMap( {{QStringLiteral( "id" ), *it}} ) );
91
92 return { QgsProfileIdentifyResults( mLayer, idsList ) };
93}
94
95QVector<QgsProfileIdentifyResults> QgsVectorLayerProfileResults::identify( const QgsProfilePoint &point, const QgsProfileIdentifyContext &context )
96{
97 QHash< QgsFeatureId, QVariantMap > features;
98 auto visitFeature = [&features]( QgsFeatureId featureId, double delta, double distance, double elevation )
99 {
100 auto it = features.find( featureId );
101 if ( it == features.end() )
102 {
103 features[ featureId ] = QVariantMap( {{QStringLiteral( "id" ), featureId },
104 {QStringLiteral( "delta" ), delta },
105 {QStringLiteral( "distance" ), distance },
106 {QStringLiteral( "elevation" ), elevation }
107 } );
108 }
109 else
110 {
111 const double currentDelta = it.value().value( QStringLiteral( "delta" ) ).toDouble();
112 if ( delta < currentDelta )
113 {
114 *it = QVariantMap( {{QStringLiteral( "id" ), featureId },
115 {QStringLiteral( "delta" ), delta },
116 {QStringLiteral( "distance" ), distance },
117 {QStringLiteral( "elevation" ), elevation }
118 } );
119 }
120 }
121 };
122
123 visitFeaturesAtPoint( point, context.maximumPointDistanceDelta, context.maximumPointElevationDelta, context.maximumSurfaceElevationDelta, visitFeature, true );
124
125 QVector< QVariantMap> attributes;
126 for ( auto it = features.constBegin(); it != features.constEnd(); ++it )
127 attributes.append( *it );
128
129 QVector<QgsProfileIdentifyResults> res;
130
131 if ( !attributes.empty() )
132 res.append( QgsProfileIdentifyResults( mLayer, attributes ) );
133
135 {
136 const QVector<QgsProfileIdentifyResults> surfaceResults = QgsAbstractProfileSurfaceResults::identify( point, context );
137 res.reserve( surfaceResults.size() );
138 for ( const QgsProfileIdentifyResults &surfaceResult : surfaceResults )
139 {
140 res.append( QgsProfileIdentifyResults( mLayer, surfaceResult.results() ) );
141 }
142 }
143
144 return res;
145}
146
147QgsProfileSnapResult QgsVectorLayerProfileResults::snapPointToIndividualFeatures( const QgsProfilePoint &point, const QgsProfileSnapContext &context )
148{
150 double bestSnapDistance = std::numeric_limits< double >::max();
151
152 auto visitFeature = [&bestSnapDistance, &res]( QgsFeatureId, double delta, double distance, double elevation )
153 {
154 if ( distance < bestSnapDistance )
155 {
156 bestSnapDistance = delta;
157 res.snappedPoint = QgsProfilePoint( distance, elevation );
158 }
159 };
160
161 visitFeaturesAtPoint( point, context.maximumPointDistanceDelta, context.maximumPointElevationDelta, context.maximumSurfaceElevationDelta, visitFeature, false );
162
163 return res;
164}
165
166void 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 )
167{
168 // TODO -- add spatial index if performance is an issue
169
170 const QgsPoint targetPoint( point.distance(), point.elevation() );
171
172 for ( auto it = features.constBegin(); it != features.constEnd(); ++it )
173 {
174 for ( const Feature &feature : it.value() )
175 {
176 const QgsRectangle featureBounds = feature.crossSectionGeometry.boundingBox();
177 if ( ( featureBounds.xMinimum() - maximumPointDistanceDelta <= point.distance() ) && ( featureBounds.xMaximum() + maximumPointDistanceDelta >= point.distance() ) )
178 {
179 switch ( feature.crossSectionGeometry.type() )
180 {
181 case Qgis::GeometryType::Point:
182 {
183 for ( auto partIt = feature.crossSectionGeometry.const_parts_begin(); partIt != feature.crossSectionGeometry.const_parts_end(); ++partIt )
184 {
185 if ( const QgsPoint *candidatePoint = qgsgeometry_cast< const QgsPoint * >( *partIt ) )
186 {
187 const double snapDistanceDelta = std::fabs( point.distance() - candidatePoint->x() );
188 if ( snapDistanceDelta > maximumPointDistanceDelta )
189 continue;
190
191 const double snapHeightDelta = std::fabs( point.elevation() - candidatePoint->y() );
192 if ( snapHeightDelta > maximumPointElevationDelta )
193 continue;
194
195 const double snapDistance = candidatePoint->distance( targetPoint );
196 visitor( feature.featureId, snapDistance, candidatePoint->x(), candidatePoint->y() );
197 }
198 }
199 break;
200 }
201
202 case Qgis::GeometryType::Line:
203 {
204 for ( auto partIt = feature.crossSectionGeometry.const_parts_begin(); partIt != feature.crossSectionGeometry.const_parts_end(); ++partIt )
205 {
206 if ( const QgsCurve *line = qgsgeometry_cast< const QgsCurve * >( *partIt ) )
207 {
208 // might be a vertical line
209 if ( const QgsLineString *lineString = qgsgeometry_cast< const QgsLineString * >( line ) )
210 {
211 if ( lineString->numPoints() == 2 && qgsDoubleNear( lineString->pointN( 0 ).x(), lineString->pointN( 1 ).x() ) )
212 {
213 const double snapDistanceDelta = std::fabs( point.distance() - lineString->pointN( 0 ).x() );
214 if ( snapDistanceDelta > maximumPointDistanceDelta )
215 continue;
216
217 const double snapHeightDelta = std::fabs( point.elevation() - lineString->pointN( 0 ).y() );
218 if ( snapHeightDelta <= maximumPointElevationDelta )
219 {
220 const double snapDistanceP1 = lineString->pointN( 0 ).distance( targetPoint );
221 visitor( feature.featureId, snapDistanceP1, lineString->pointN( 0 ).x(), lineString->pointN( 0 ).y() );
222 }
223
224 const double snapHeightDelta2 = std::fabs( point.elevation() - lineString->pointN( 1 ).y() );
225 if ( snapHeightDelta2 <= maximumPointElevationDelta )
226 {
227 const double snapDistanceP2 = lineString->pointN( 1 ).distance( targetPoint );
228 visitor( feature.featureId, snapDistanceP2, lineString->pointN( 1 ).x(), lineString->pointN( 1 ).y() );
229 }
230
231 if ( visitWithin )
232 {
233 double elevation1 = lineString->pointN( 0 ).y();
234 double elevation2 = lineString->pointN( 1 ).y();
235 if ( elevation1 > elevation2 )
236 std::swap( elevation1, elevation2 );
237
238 if ( point.elevation() > elevation1 && point.elevation() < elevation2 )
239 {
240 const double snapDistance = std::fabs( lineString->pointN( 0 ).x() - point.distance() );
241 visitor( feature.featureId, snapDistance, lineString->pointN( 0 ).x(), point.elevation() );
242 }
243 }
244 continue;
245 }
246 }
247
248 const QgsRectangle partBounds = ( *partIt )->boundingBox();
249 if ( point.distance() < partBounds.xMinimum() - maximumPointDistanceDelta || point.distance() > partBounds.xMaximum() + maximumPointDistanceDelta )
250 continue;
251
252 const double snappedDistance = point.distance() < partBounds.xMinimum() ? partBounds.xMinimum()
253 : point.distance() > partBounds.xMaximum() ? partBounds.xMaximum() : point.distance();
254
255 const QgsGeometry cutLine( new QgsLineString( QgsPoint( snappedDistance, minZ ), QgsPoint( snappedDistance, maxZ ) ) );
256 QgsGeos cutLineGeos( cutLine.constGet() );
257
258 const QgsGeometry points( cutLineGeos.intersection( line ) );
259
260 for ( auto vertexIt = points.vertices_begin(); vertexIt != points.vertices_end(); ++vertexIt )
261 {
262 const double snapHeightDelta = std::fabs( point.elevation() - ( *vertexIt ).y() );
263 if ( snapHeightDelta > maximumSurfaceElevationDelta )
264 continue;
265
266 const double snapDistance = ( *vertexIt ).distance( targetPoint );
267 visitor( feature.featureId, snapDistance, ( *vertexIt ).x(), ( *vertexIt ).y() );
268 }
269 }
270 }
271 break;
272 }
273
274 case Qgis::GeometryType::Polygon:
275 {
276 if ( visitWithin )
277 {
278 if ( feature.crossSectionGeometry.intersects( QgsGeometry::fromPointXY( QgsPointXY( point.distance(), point.elevation() ) ) ) )
279 {
280 visitor( feature.featureId, 0, point.distance(), point.elevation() );
281 break;
282 }
283 }
284 for ( auto partIt = feature.crossSectionGeometry.const_parts_begin(); partIt != feature.crossSectionGeometry.const_parts_end(); ++partIt )
285 {
286 if ( const QgsCurve *exterior = qgsgeometry_cast< const QgsPolygon * >( *partIt )->exteriorRing() )
287 {
288 const QgsRectangle partBounds = ( *partIt )->boundingBox();
289 if ( point.distance() < partBounds.xMinimum() - maximumPointDistanceDelta || point.distance() > partBounds.xMaximum() + maximumPointDistanceDelta )
290 continue;
291
292 const double snappedDistance = point.distance() < partBounds.xMinimum() ? partBounds.xMinimum()
293 : point.distance() > partBounds.xMaximum() ? partBounds.xMaximum() : point.distance();
294
295 const QgsGeometry cutLine( new QgsLineString( QgsPoint( snappedDistance, minZ ), QgsPoint( snappedDistance, maxZ ) ) );
296 QgsGeos cutLineGeos( cutLine.constGet() );
297
298 const QgsGeometry points( cutLineGeos.intersection( exterior ) );
299 for ( auto vertexIt = points.vertices_begin(); vertexIt != points.vertices_end(); ++vertexIt )
300 {
301 const double snapHeightDelta = std::fabs( point.elevation() - ( *vertexIt ).y() );
302 if ( snapHeightDelta > maximumSurfaceElevationDelta )
303 continue;
304
305 const double snapDistance = ( *vertexIt ).distance( targetPoint );
306 visitor( feature.featureId, snapDistance, ( *vertexIt ).x(), ( *vertexIt ).y() );
307 }
308 }
309 }
310 break;
311 }
312 case Qgis::GeometryType::Unknown:
313 case Qgis::GeometryType::Null:
314 break;
315 }
316 }
317 }
318 }
319}
320
321void QgsVectorLayerProfileResults::visitFeaturesInRange( const QgsDoubleRange &distanceRange, const QgsDoubleRange &elevationRange, const std::function<void ( QgsFeatureId )> &visitor )
322{
323 // TODO -- add spatial index if performance is an issue
324 const QgsRectangle profileRange( distanceRange.lower(), elevationRange.lower(), distanceRange.upper(), elevationRange.upper() );
325 const QgsGeometry profileRangeGeometry = QgsGeometry::fromRect( profileRange );
326 QgsGeos profileRangeGeos( profileRangeGeometry.constGet() );
327 profileRangeGeos.prepareGeometry();
328
329 for ( auto it = features.constBegin(); it != features.constEnd(); ++it )
330 {
331 for ( const Feature &feature : it.value() )
332 {
333 if ( feature.crossSectionGeometry.boundingBoxIntersects( profileRange ) )
334 {
335 switch ( feature.crossSectionGeometry.type() )
336 {
337 case Qgis::GeometryType::Point:
338 {
339 for ( auto partIt = feature.crossSectionGeometry.const_parts_begin(); partIt != feature.crossSectionGeometry.const_parts_end(); ++partIt )
340 {
341 if ( const QgsPoint *candidatePoint = qgsgeometry_cast< const QgsPoint * >( *partIt ) )
342 {
343 if ( profileRange.contains( candidatePoint->x(), candidatePoint->y() ) )
344 {
345 visitor( feature.featureId );
346 }
347 }
348 }
349 break;
350 }
351
352 case Qgis::GeometryType::Line:
353 case Qgis::GeometryType::Polygon:
354 {
355 if ( profileRangeGeos.intersects( feature.crossSectionGeometry.constGet() ) )
356 {
357 visitor( feature.featureId );
358 }
359 break;
360 }
361
362 case Qgis::GeometryType::Unknown:
363 case Qgis::GeometryType::Null:
364 break;
365 }
366 }
367 }
368 }
369}
370
372{
373 const QgsExpressionContextScopePopper scopePopper( context.renderContext().expressionContext(), mLayer ? mLayer->createExpressionContextScope() : nullptr );
374 switch ( profileType )
375 {
377 renderResultsAsIndividualFeatures( context );
378 break;
382 renderMarkersOverContinuousSurfacePlot( context );
383 break;
384 }
385}
386
387void QgsVectorLayerProfileResults::renderResultsAsIndividualFeatures( QgsProfileRenderContext &context )
388{
389 QPainter *painter = context.renderContext().painter();
390 if ( !painter )
391 return;
392
393 const QgsScopedQPainterState painterState( painter );
394
395 painter->setBrush( Qt::NoBrush );
396 painter->setPen( Qt::NoPen );
397
398 const double minDistance = context.distanceRange().lower();
399 const double maxDistance = context.distanceRange().upper();
400 const double minZ = context.elevationRange().lower();
401 const double maxZ = context.elevationRange().upper();
402
403 const QRectF visibleRegion( minDistance, minZ, maxDistance - minDistance, maxZ - minZ );
404 QPainterPath clipPath;
405 clipPath.addPolygon( context.worldTransform().map( visibleRegion ) );
406 painter->setClipPath( clipPath, Qt::ClipOperation::IntersectClip );
407
408 const QgsRectangle clipPathRect( clipPath.boundingRect() );
409
410 auto renderResult = [&context, &clipPathRect]( const Feature & profileFeature, QgsMarkerSymbol * markerSymbol, QgsLineSymbol * lineSymbol, QgsFillSymbol * fillSymbol )
411 {
412 if ( profileFeature.crossSectionGeometry.isEmpty() )
413 return;
414
415 QgsGeometry transformed = profileFeature.crossSectionGeometry;
416 transformed.transform( context.worldTransform() );
417
418 if ( !transformed.boundingBoxIntersects( clipPathRect ) )
419 return;
420
421 // we can take some shortcuts here, because we know that the geometry will already be segmentized and can't be a curved type
422 switch ( transformed.type() )
423 {
424 case Qgis::GeometryType::Point:
425 {
426 if ( const QgsPoint *point = qgsgeometry_cast< const QgsPoint * >( transformed.constGet() ) )
427 {
428 markerSymbol->renderPoint( QPointF( point->x(), point->y() ), nullptr, context.renderContext() );
429 }
430 else if ( const QgsMultiPoint *multipoint = qgsgeometry_cast< const QgsMultiPoint * >( transformed.constGet() ) )
431 {
432 const int numGeometries = multipoint->numGeometries();
433 for ( int i = 0; i < numGeometries; ++i )
434 {
435 markerSymbol->renderPoint( QPointF( multipoint->pointN( i )->x(), multipoint->pointN( i )->y() ), nullptr, context.renderContext() );
436 }
437 }
438 break;
439 }
440
441 case Qgis::GeometryType::Line:
442 {
443 if ( const QgsLineString *line = qgsgeometry_cast< const QgsLineString * >( transformed.constGet() ) )
444 {
445 lineSymbol->renderPolyline( line->asQPolygonF(), nullptr, context.renderContext() );
446 }
447 else if ( const QgsMultiLineString *multiLinestring = qgsgeometry_cast< const QgsMultiLineString * >( transformed.constGet() ) )
448 {
449 const int numGeometries = multiLinestring->numGeometries();
450 for ( int i = 0; i < numGeometries; ++i )
451 {
452 lineSymbol->renderPolyline( multiLinestring->lineStringN( i )->asQPolygonF(), nullptr, context.renderContext() );
453 }
454 }
455 break;
456 }
457
458 case Qgis::GeometryType::Polygon:
459 {
460 if ( const QgsPolygon *polygon = qgsgeometry_cast< const QgsPolygon * >( transformed.constGet() ) )
461 {
462 if ( const QgsCurve *exterior = polygon->exteriorRing() )
463 fillSymbol->renderPolygon( exterior->asQPolygonF(), nullptr, nullptr, context.renderContext() );
464 }
465 else if ( const QgsMultiPolygon *multiPolygon = qgsgeometry_cast< const QgsMultiPolygon * >( transformed.constGet() ) )
466 {
467 const int numGeometries = multiPolygon->numGeometries();
468 for ( int i = 0; i < numGeometries; ++i )
469 {
470 fillSymbol->renderPolygon( multiPolygon->polygonN( i )->exteriorRing()->asQPolygonF(), nullptr, nullptr, context.renderContext() );
471 }
472 }
473 break;
474 }
475
476 case Qgis::GeometryType::Unknown:
477 case Qgis::GeometryType::Null:
478 return;
479 }
480 };
481
483 req.setFilterFids( qgis::listToSet( features.keys() ) );
484
485 if ( respectLayerSymbology && mLayer && mLayer->renderer() )
486 {
487 std::unique_ptr< QgsFeatureRenderer > renderer( mLayer->renderer()->clone() );
488 renderer->startRender( context.renderContext(), mLayer->fields() );
489
490 // if we are respecting the layer's symbology then we'll fire off a feature request and iterate through
491 // features from the profile, rendering each in turn
492 QSet<QString> attributes = renderer->usedAttributes( context.renderContext() );
493
494 std::unique_ptr< QgsMarkerSymbol > marker( mMarkerSymbol->clone() );
495 std::unique_ptr< QgsLineSymbol > line( mLineSymbol->clone() );
496 std::unique_ptr< QgsFillSymbol > fill( mFillSymbol->clone() );
497 attributes.unite( marker->usedAttributes( context.renderContext() ) );
498 attributes.unite( line->usedAttributes( context.renderContext() ) );
499 attributes.unite( fill->usedAttributes( context.renderContext() ) );
500
501 req.setSubsetOfAttributes( attributes, mLayer->fields() );
502
503 QgsFeature feature;
504 QgsFeatureIterator it = mLayer->getFeatures( req );
505 while ( it.nextFeature( feature ) )
506 {
507 context.renderContext().expressionContext().setFeature( feature );
508 QgsSymbol *rendererSymbol = renderer->symbolForFeature( feature, context.renderContext() );
509 if ( !rendererSymbol )
510 continue;
511
512 marker->setColor( rendererSymbol->color() );
513 marker->setOpacity( rendererSymbol->opacity() );
514 line->setColor( rendererSymbol->color() );
515 line->setOpacity( rendererSymbol->opacity() );
516 fill->setColor( rendererSymbol->color() );
517 fill->setOpacity( rendererSymbol->opacity() );
518
519 marker->startRender( context.renderContext() );
520 line->startRender( context.renderContext() );
521 fill->startRender( context.renderContext() );
522
523 const QVector< Feature > profileFeatures = features.value( feature.id() );
524 for ( const Feature &profileFeature : profileFeatures )
525 {
526 renderResult( profileFeature,
527 rendererSymbol->type() == Qgis::SymbolType::Marker ? qgis::down_cast< QgsMarkerSymbol * >( rendererSymbol ) : marker.get(),
528 rendererSymbol->type() == Qgis::SymbolType::Line ? qgis::down_cast< QgsLineSymbol * >( rendererSymbol ) : line.get(),
529 rendererSymbol->type() == Qgis::SymbolType::Fill ? qgis::down_cast< QgsFillSymbol * >( rendererSymbol ) : fill.get() );
530 }
531
532 marker->stopRender( context.renderContext() );
533 line->stopRender( context.renderContext() );
534 fill->stopRender( context.renderContext() );
535 }
536
537 renderer->stopRender( context.renderContext() );
538 }
539 else
540 {
541 QSet<QString> attributes;
542 attributes.unite( mMarkerSymbol->usedAttributes( context.renderContext() ) );
543 attributes.unite( mFillSymbol->usedAttributes( context.renderContext() ) );
544 attributes.unite( mLineSymbol->usedAttributes( context.renderContext() ) );
545
546 mMarkerSymbol->startRender( context.renderContext() );
547 mFillSymbol->startRender( context.renderContext() );
548 mLineSymbol->startRender( context.renderContext() );
549 req.setSubsetOfAttributes( attributes, mLayer->fields() );
550
551 QgsFeature feature;
552 QgsFeatureIterator it = mLayer->getFeatures( req );
553 while ( it.nextFeature( feature ) )
554 {
555 context.renderContext().expressionContext().setFeature( feature );
556 const QVector< Feature > profileFeatures = features.value( feature.id() );
557 for ( const Feature &profileFeature : profileFeatures )
558 {
559 renderResult( profileFeature, mMarkerSymbol.get(), mLineSymbol.get(), mFillSymbol.get() );
560 }
561 }
562 mMarkerSymbol->stopRender( context.renderContext() );
563 mFillSymbol->stopRender( context.renderContext() );
564 mLineSymbol->stopRender( context.renderContext() );
565 }
566}
567
568void QgsVectorLayerProfileResults::renderMarkersOverContinuousSurfacePlot( QgsProfileRenderContext &context )
569{
570 QPainter *painter = context.renderContext().painter();
571 if ( !painter )
572 return;
573
574 const QgsScopedQPainterState painterState( painter );
575
576 painter->setBrush( Qt::NoBrush );
577 painter->setPen( Qt::NoPen );
578
579 const double minDistance = context.distanceRange().lower();
580 const double maxDistance = context.distanceRange().upper();
581 const double minZ = context.elevationRange().lower();
582 const double maxZ = context.elevationRange().upper();
583
584 const QRectF visibleRegion( minDistance, minZ, maxDistance - minDistance, maxZ - minZ );
585 QPainterPath clipPath;
586 clipPath.addPolygon( context.worldTransform().map( visibleRegion ) );
587 painter->setClipPath( clipPath, Qt::ClipOperation::IntersectClip );
588
589 mMarkerSymbol->startRender( context.renderContext() );
590
591 for ( auto pointIt = mDistanceToHeightMap.constBegin(); pointIt != mDistanceToHeightMap.constEnd(); ++pointIt )
592 {
593 if ( std::isnan( pointIt.value() ) )
594 continue;
595
596 mMarkerSymbol->renderPoint( context.worldTransform().map( QPointF( pointIt.key(), pointIt.value() ) ), nullptr, context.renderContext() );
597 }
598 mMarkerSymbol->stopRender( context.renderContext() );
599}
600
602{
604 const QgsVectorLayerProfileGenerator *vlGenerator = qgis::down_cast< const QgsVectorLayerProfileGenerator * >( generator );
605
606 profileType = vlGenerator->mType;
607 respectLayerSymbology = vlGenerator->mRespectLayerSymbology;
608 mMarkerSymbol.reset( vlGenerator->mProfileMarkerSymbol->clone() );
609 mShowMarkerSymbolInSurfacePlots = vlGenerator->mShowMarkerSymbolInSurfacePlots;
610}
611
612//
613// QgsVectorLayerProfileGenerator
614//
615
617 : mId( layer->id() )
618 , mFeedback( std::make_unique< QgsFeedback >() )
619 , mProfileCurve( request.profileCurve() ? request.profileCurve()->clone() : nullptr )
620 , mTerrainProvider( request.terrainProvider() ? request.terrainProvider()->clone() : nullptr )
621 , mTolerance( request.tolerance() )
622 , mSourceCrs( layer->crs() )
623 , mTargetCrs( request.crs() )
624 , mTransformContext( request.transformContext() )
625 , mExtent( layer->extent() )
626 , mSource( std::make_unique< QgsVectorLayerFeatureSource >( layer ) )
627 , mOffset( layer->elevationProperties()->zOffset() )
628 , mScale( layer->elevationProperties()->zScale() )
629 , mType( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->type() )
630 , mClamping( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->clamping() )
631 , mBinding( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->binding() )
632 , mExtrusionEnabled( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->extrusionEnabled() )
633 , mExtrusionHeight( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->extrusionHeight() )
634 , mExpressionContext( request.expressionContext() )
635 , mFields( layer->fields() )
636 , mDataDefinedProperties( layer->elevationProperties()->dataDefinedProperties() )
637 , mWkbType( layer->wkbType() )
638 , mRespectLayerSymbology( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->respectLayerSymbology() )
639 , mProfileMarkerSymbol( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->profileMarkerSymbol()->clone() )
640 , mShowMarkerSymbolInSurfacePlots( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->showMarkerSymbolInSurfacePlots() )
641 , mLayer( layer )
642{
643 if ( mTerrainProvider )
644 mTerrainProvider->prepare(); // must be done on main thread
645
646 mSymbology = qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->profileSymbology();
647 mLineSymbol.reset( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->profileLineSymbol()->clone() );
648 mFillSymbol.reset( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->profileFillSymbol()->clone() );
649}
650
652{
653 return mId;
654}
655
657
659{
660 if ( !mProfileCurve || mFeedback->isCanceled() )
661 return false;
662
663 // we need to transform the profile curve to the vector's CRS
664 mTransformedCurve.reset( mProfileCurve->clone() );
665 mLayerToTargetTransform = QgsCoordinateTransform( mSourceCrs, mTargetCrs, mTransformContext );
666 if ( mTerrainProvider )
667 mTargetToTerrainProviderTransform = QgsCoordinateTransform( mTargetCrs, mTerrainProvider->crs(), mTransformContext );
668
669 try
670 {
671 mTransformedCurve->transform( mLayerToTargetTransform, Qgis::TransformDirection::Reverse );
672 }
673 catch ( QgsCsException & )
674 {
675 QgsDebugMsg( QStringLiteral( "Error transforming profile line to vector CRS" ) );
676 return false;
677 }
678
679 const QgsRectangle profileCurveBoundingBox = mTransformedCurve->boundingBox();
680 if ( !profileCurveBoundingBox.intersects( mExtent ) )
681 return false;
682
683 if ( mFeedback->isCanceled() )
684 return false;
685
686 mResults = std::make_unique< QgsVectorLayerProfileResults >();
687 mResults->mLayer = mLayer;
688 mResults->copyPropertiesFromGenerator( this );
689
690 mProfileCurveEngine.reset( new QgsGeos( mProfileCurve.get() ) );
691 mProfileCurveEngine->prepareGeometry();
692
693 mDataDefinedProperties.prepare( mExpressionContext );
694
695 if ( mFeedback->isCanceled() )
696 return false;
697
698 switch ( QgsWkbTypes::geometryType( mWkbType ) )
699 {
700 case Qgis::GeometryType::Point:
701 if ( !generateProfileForPoints() )
702 return false;
703 break;
704
705 case Qgis::GeometryType::Line:
706 if ( !generateProfileForLines() )
707 return false;
708 break;
709
710 case Qgis::GeometryType::Polygon:
711 if ( !generateProfileForPolygons() )
712 return false;
713 break;
714
715 case Qgis::GeometryType::Unknown:
716 case Qgis::GeometryType::Null:
717 return false;
718 }
719
720 return true;
721}
722
724{
725 return mResults.release();
726}
727
729{
730 return mFeedback.get();
731}
732
733bool QgsVectorLayerProfileGenerator::generateProfileForPoints()
734{
735 // get features from layer
736 QgsFeatureRequest request;
737 request.setDestinationCrs( mTargetCrs, mTransformContext );
738 request.setDistanceWithin( QgsGeometry( mProfileCurve->clone() ), mTolerance );
739 request.setSubsetOfAttributes( mDataDefinedProperties.referencedFields( mExpressionContext ), mFields );
740 request.setFeedback( mFeedback.get() );
741
742 // our feature request is using the optimised distance within check (allowing use of spatial index)
743 // BUT this will also include points which are within the tolerance distance before/after the end of line.
744 // So we also need to double check that they fall within the flat buffered curve too.
745 std::unique_ptr< QgsAbstractGeometry > bufferedCurve( mProfileCurveEngine->buffer( mTolerance, 8, Qgis::EndCapStyle::Flat, Qgis::JoinStyle::Round, 2 ) );
746 QgsGeos bufferedCurveEngine( bufferedCurve.get() );
747 bufferedCurveEngine.prepareGeometry();
748
749 auto processPoint = [this, &bufferedCurveEngine]( const QgsFeature & feature, const QgsPoint * point )
750 {
751 if ( !bufferedCurveEngine.intersects( point ) )
752 return;
753
754 const double offset = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::ZOffset, mExpressionContext, mOffset );
755
756 const double height = featureZToHeight( point->x(), point->y(), point->z(), offset );
757 mResults->mRawPoints.append( QgsPoint( point->x(), point->y(), height ) );
758 mResults->minZ = std::min( mResults->minZ, height );
759 mResults->maxZ = std::max( mResults->maxZ, height );
760
761 QString lastError;
762 const double distance = mProfileCurveEngine->lineLocatePoint( *point, &lastError );
763 mResults->mDistanceToHeightMap.insert( distance, height );
764
766 resultFeature.featureId = feature.id();
767 if ( mExtrusionEnabled )
768 {
769 const double extrusion = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::ExtrusionHeight, mExpressionContext, mExtrusionHeight );
770
771 resultFeature.geometry = QgsGeometry( new QgsLineString( QgsPoint( point->x(), point->y(), height ),
772 QgsPoint( point->x(), point->y(), height + extrusion ) ) );
773 resultFeature.crossSectionGeometry = QgsGeometry( new QgsLineString( QgsPoint( distance, height ),
774 QgsPoint( distance, height + extrusion ) ) );
775 mResults->minZ = std::min( mResults->minZ, height + extrusion );
776 mResults->maxZ = std::max( mResults->maxZ, height + extrusion );
777 }
778 else
779 {
780 resultFeature.geometry = QgsGeometry( new QgsPoint( point->x(), point->y(), height ) );
781 resultFeature.crossSectionGeometry = QgsGeometry( new QgsPoint( distance, height ) );
782 }
783 mResults->features[resultFeature.featureId].append( resultFeature );
784 };
785
786 QgsFeature feature;
787 QgsFeatureIterator it = mSource->getFeatures( request );
788 while ( it.nextFeature( feature ) )
789 {
790 if ( mFeedback->isCanceled() )
791 return false;
792
793 mExpressionContext.setFeature( feature );
794
795 const QgsGeometry g = feature.geometry();
796 if ( g.isMultipart() )
797 {
798 for ( auto it = g.const_parts_begin(); it != g.const_parts_end(); ++it )
799 {
800 processPoint( feature, qgsgeometry_cast< const QgsPoint * >( *it ) );
801 }
802 }
803 else
804 {
805 processPoint( feature, qgsgeometry_cast< const QgsPoint * >( g.constGet() ) );
806 }
807 }
808 return true;
809}
810
811bool QgsVectorLayerProfileGenerator::generateProfileForLines()
812{
813 // get features from layer
814 QgsFeatureRequest request;
815 request.setDestinationCrs( mTargetCrs, mTransformContext );
816 request.setFilterRect( mProfileCurve->boundingBox() );
817 request.setSubsetOfAttributes( mDataDefinedProperties.referencedFields( mExpressionContext ), mFields );
818 request.setFeedback( mFeedback.get() );
819
820 auto processCurve = [this]( const QgsFeature & feature, const QgsCurve * curve )
821 {
822 QString error;
823 std::unique_ptr< QgsAbstractGeometry > intersection( mProfileCurveEngine->intersection( curve, &error ) );
824 if ( !intersection )
825 return;
826
827 if ( mFeedback->isCanceled() )
828 return;
829
830 QgsGeos curveGeos( curve );
831 curveGeos.prepareGeometry();
832
833 if ( mFeedback->isCanceled() )
834 return;
835
836 for ( auto it = intersection->const_parts_begin(); it != intersection->const_parts_end(); ++it )
837 {
838 if ( mFeedback->isCanceled() )
839 return;
840
841 if ( const QgsPoint *intersectionPoint = qgsgeometry_cast< const QgsPoint * >( *it ) )
842 {
843 // unfortunately we need to do some work to interpolate the z value for the line -- GEOS doesn't give us this
844 const double distance = curveGeos.lineLocatePoint( *intersectionPoint, &error );
845 std::unique_ptr< QgsPoint > interpolatedPoint( curve->interpolatePoint( distance ) );
846
847 const double offset = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::ZOffset, mExpressionContext, mOffset );
848
849 const double height = featureZToHeight( interpolatedPoint->x(), interpolatedPoint->y(), interpolatedPoint->z(), offset );
850 mResults->mRawPoints.append( QgsPoint( interpolatedPoint->x(), interpolatedPoint->y(), height ) );
851 mResults->minZ = std::min( mResults->minZ, height );
852 mResults->maxZ = std::max( mResults->maxZ, height );
853
854 const double distanceAlongProfileCurve = mProfileCurveEngine->lineLocatePoint( *interpolatedPoint, &error );
855 mResults->mDistanceToHeightMap.insert( distanceAlongProfileCurve, height );
856
858 resultFeature.featureId = feature.id();
859 if ( mExtrusionEnabled )
860 {
861 const double extrusion = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::ExtrusionHeight, mExpressionContext, mExtrusionHeight );
862
863 resultFeature.geometry = QgsGeometry( new QgsLineString( QgsPoint( interpolatedPoint->x(), interpolatedPoint->y(), height ),
864 QgsPoint( interpolatedPoint->x(), interpolatedPoint->y(), height + extrusion ) ) );
865 resultFeature.crossSectionGeometry = QgsGeometry( new QgsLineString( QgsPoint( distanceAlongProfileCurve, height ),
866 QgsPoint( distanceAlongProfileCurve, height + extrusion ) ) );
867 mResults->minZ = std::min( mResults->minZ, height + extrusion );
868 mResults->maxZ = std::max( mResults->maxZ, height + extrusion );
869 }
870 else
871 {
872 resultFeature.geometry = QgsGeometry( new QgsPoint( interpolatedPoint->x(), interpolatedPoint->y(), height ) );
873 resultFeature.crossSectionGeometry = QgsGeometry( new QgsPoint( distanceAlongProfileCurve, height ) );
874 }
875 mResults->features[resultFeature.featureId].append( resultFeature );
876 }
877 }
878 };
879
880 QgsFeature feature;
881 QgsFeatureIterator it = mSource->getFeatures( request );
882 while ( it.nextFeature( feature ) )
883 {
884 if ( mFeedback->isCanceled() )
885 return false;
886
887 if ( !mProfileCurveEngine->intersects( feature.geometry().constGet() ) )
888 continue;
889
890 mExpressionContext.setFeature( feature );
891
892 const QgsGeometry g = feature.geometry();
893 if ( g.isMultipart() )
894 {
895 for ( auto it = g.const_parts_begin(); it != g.const_parts_end(); ++it )
896 {
897 if ( !mProfileCurveEngine->intersects( *it ) )
898 continue;
899
900 processCurve( feature, qgsgeometry_cast< const QgsCurve * >( *it ) );
901 }
902 }
903 else
904 {
905 processCurve( feature, qgsgeometry_cast< const QgsCurve * >( g.constGet() ) );
906 }
907 }
908 return true;
909}
910
911bool QgsVectorLayerProfileGenerator::generateProfileForPolygons()
912{
913 // get features from layer
914 QgsFeatureRequest request;
915 request.setDestinationCrs( mTargetCrs, mTransformContext );
916 request.setFilterRect( mProfileCurve->boundingBox() );
917 request.setSubsetOfAttributes( mDataDefinedProperties.referencedFields( mExpressionContext ), mFields );
918 request.setFeedback( mFeedback.get() );
919
920 auto interpolatePointOnTriangle = []( const QgsPolygon * triangle, double x, double y ) -> QgsPoint
921 {
922 QgsPoint p1, p2, p3;
924 triangle->exteriorRing()->pointAt( 0, p1, vt );
925 triangle->exteriorRing()->pointAt( 1, p2, vt );
926 triangle->exteriorRing()->pointAt( 2, p3, vt );
927 const double z = QgsMeshLayerUtils::interpolateFromVerticesData( p1, p2, p3, p1.z(), p2.z(), p3.z(), QgsPointXY( x, y ) );
928 return QgsPoint( x, y, z );
929 };
930
931 std::function< void( const QgsPolygon *triangle, const QgsAbstractGeometry *intersect, QVector< QgsGeometry > &, QVector< QgsGeometry > & ) > processTriangleLineIntersect;
932 processTriangleLineIntersect = [this, &interpolatePointOnTriangle, &processTriangleLineIntersect]( const QgsPolygon * triangle, const QgsAbstractGeometry * intersect, QVector< QgsGeometry > &transformedParts, QVector< QgsGeometry > &crossSectionParts )
933 {
934 // intersect may be a (multi)point or (multi)linestring
935 switch ( QgsWkbTypes::geometryType( intersect->wkbType() ) )
936 {
937 case Qgis::GeometryType::Point:
938 if ( const QgsMultiPoint *mp = qgsgeometry_cast< const QgsMultiPoint * >( intersect ) )
939 {
940 const int numPoint = mp->numGeometries();
941 for ( int i = 0; i < numPoint; ++i )
942 {
943 processTriangleLineIntersect( triangle, mp->geometryN( i ), transformedParts, crossSectionParts );
944 }
945 }
946 else if ( const QgsPoint *p = qgsgeometry_cast< const QgsPoint * >( intersect ) )
947 {
948 const QgsPoint interpolatedPoint = interpolatePointOnTriangle( triangle, p->x(), p->y() );
949 mResults->mRawPoints.append( interpolatedPoint );
950 mResults->minZ = std::min( mResults->minZ, interpolatedPoint.z() );
951 mResults->maxZ = std::max( mResults->maxZ, interpolatedPoint.z() );
952
953 QString lastError;
954 const double distance = mProfileCurveEngine->lineLocatePoint( *p, &lastError );
955 mResults->mDistanceToHeightMap.insert( distance, interpolatedPoint.z() );
956
957 if ( mExtrusionEnabled )
958 {
959 const double extrusion = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::ExtrusionHeight, mExpressionContext, mExtrusionHeight );
960
961 transformedParts.append( QgsGeometry( new QgsLineString( interpolatedPoint,
962 QgsPoint( interpolatedPoint.x(), interpolatedPoint.y(), interpolatedPoint.z() + extrusion ) ) ) );
963 crossSectionParts.append( QgsGeometry( new QgsLineString( QgsPoint( distance, interpolatedPoint.z() ),
964 QgsPoint( distance, interpolatedPoint.z() + extrusion ) ) ) );
965 mResults->minZ = std::min( mResults->minZ, interpolatedPoint.z() + extrusion );
966 mResults->maxZ = std::max( mResults->maxZ, interpolatedPoint.z() + extrusion );
967 }
968 else
969 {
970 transformedParts.append( QgsGeometry( new QgsPoint( interpolatedPoint ) ) );
971 crossSectionParts.append( QgsGeometry( new QgsPoint( distance, interpolatedPoint.z() ) ) );
972 }
973 }
974 break;
975 case Qgis::GeometryType::Line:
976 if ( const QgsMultiLineString *ml = qgsgeometry_cast< const QgsMultiLineString * >( intersect ) )
977 {
978 const int numLines = ml->numGeometries();
979 for ( int i = 0; i < numLines; ++i )
980 {
981 processTriangleLineIntersect( triangle, ml->geometryN( i ), transformedParts, crossSectionParts );
982 }
983 }
984 else if ( const QgsLineString *ls = qgsgeometry_cast< const QgsLineString * >( intersect ) )
985 {
986 const int numPoints = ls->numPoints();
987 QVector< double > newX;
988 newX.resize( numPoints );
989 QVector< double > newY;
990 newY.resize( numPoints );
991 QVector< double > newZ;
992 newZ.resize( numPoints );
993 QVector< double > newDistance;
994 newDistance.resize( numPoints );
995
996 const double *inX = ls->xData();
997 const double *inY = ls->yData();
998 double *outX = newX.data();
999 double *outY = newY.data();
1000 double *outZ = newZ.data();
1001 double *outDistance = newDistance.data();
1002
1003 QVector< double > extrudedZ;
1004 double *extZOut = nullptr;
1005 double extrusion = 0;
1006 if ( mExtrusionEnabled )
1007 {
1008 extrudedZ.resize( numPoints );
1009 extZOut = extrudedZ.data();
1010
1011 extrusion = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::ExtrusionHeight, mExpressionContext, mExtrusionHeight );
1012 }
1013
1014 QString lastError;
1015 for ( int i = 0 ; i < numPoints; ++i )
1016 {
1017 double x = *inX++;
1018 double y = *inY++;
1019
1020 QgsPoint interpolatedPoint = interpolatePointOnTriangle( triangle, x, y );
1021 *outX++ = x;
1022 *outY++ = y;
1023 *outZ++ = interpolatedPoint.z();
1024 if ( extZOut )
1025 *extZOut++ = interpolatedPoint.z() + extrusion;
1026
1027 mResults->mRawPoints.append( interpolatedPoint );
1028 mResults->minZ = std::min( mResults->minZ, interpolatedPoint.z() );
1029 mResults->maxZ = std::max( mResults->maxZ, interpolatedPoint.z() );
1030 if ( mExtrusionEnabled )
1031 {
1032 mResults->minZ = std::min( mResults->minZ, interpolatedPoint.z() + extrusion );
1033 mResults->maxZ = std::max( mResults->maxZ, interpolatedPoint.z() + extrusion );
1034 }
1035
1036 const double distance = mProfileCurveEngine->lineLocatePoint( interpolatedPoint, &lastError );
1037 *outDistance++ = distance;
1038
1039 mResults->mDistanceToHeightMap.insert( distance, interpolatedPoint.z() );
1040 }
1041
1042 if ( mExtrusionEnabled )
1043 {
1044 std::unique_ptr< QgsLineString > ring = std::make_unique< QgsLineString >( newX, newY, newZ );
1045 std::unique_ptr< QgsLineString > extrudedRing = std::make_unique< QgsLineString >( newX, newY, extrudedZ );
1046 std::unique_ptr< QgsLineString > reversedExtrusion( extrudedRing->reversed() );
1047 ring->append( reversedExtrusion.get() );
1048 ring->close();
1049 transformedParts.append( QgsGeometry( new QgsPolygon( ring.release() ) ) );
1050
1051
1052 std::unique_ptr< QgsLineString > distanceVHeightRing = std::make_unique< QgsLineString >( newDistance, newZ );
1053 std::unique_ptr< QgsLineString > extrudedDistanceVHeightRing = std::make_unique< QgsLineString >( newDistance, extrudedZ );
1054 std::unique_ptr< QgsLineString > reversedDistanceVHeightExtrusion( extrudedDistanceVHeightRing->reversed() );
1055 distanceVHeightRing->append( reversedDistanceVHeightExtrusion.get() );
1056 distanceVHeightRing->close();
1057 crossSectionParts.append( QgsGeometry( new QgsPolygon( distanceVHeightRing.release() ) ) );
1058 }
1059 else
1060 {
1061 transformedParts.append( QgsGeometry( new QgsLineString( newX, newY, newZ ) ) );
1062 crossSectionParts.append( QgsGeometry( new QgsLineString( newDistance, newZ ) ) );
1063 }
1064 }
1065 break;
1066
1067 case Qgis::GeometryType::Polygon:
1068 case Qgis::GeometryType::Unknown:
1069 case Qgis::GeometryType::Null:
1070 return;
1071 }
1072 };
1073
1074 auto processPolygon = [this, &processTriangleLineIntersect]( const QgsCurvePolygon * polygon, QVector< QgsGeometry > &transformedParts, QVector< QgsGeometry > &crossSectionParts, double offset )
1075 {
1076 std::unique_ptr< QgsPolygon > clampedPolygon;
1077 if ( const QgsPolygon *p = qgsgeometry_cast< const QgsPolygon * >( polygon ) )
1078 {
1079 clampedPolygon.reset( p->clone() );
1080 }
1081 else
1082 {
1083 clampedPolygon.reset( qgsgeometry_cast< QgsPolygon * >( polygon->segmentize() ) );
1084 }
1085 clampAltitudes( clampedPolygon.get(), offset );
1086
1087 if ( mFeedback->isCanceled() )
1088 return;
1089
1090 const QgsRectangle bounds = clampedPolygon->boundingBox();
1091 QgsTessellator t( bounds, false, false, false, false );
1092 t.addPolygon( *clampedPolygon, 0 );
1093
1094 QgsGeometry tessellation( t.asMultiPolygon() );
1095 if ( mFeedback->isCanceled() )
1096 return;
1097
1098 tessellation.translate( bounds.xMinimum(), bounds.yMinimum() );
1099
1100 // iterate through the tessellation, finding triangles which intersect the line
1101 const int numTriangles = qgsgeometry_cast< const QgsMultiPolygon * >( tessellation.constGet() )->numGeometries();
1102 for ( int i = 0; i < numTriangles; ++i )
1103 {
1104 if ( mFeedback->isCanceled() )
1105 return;
1106
1107 const QgsPolygon *triangle = qgsgeometry_cast< const QgsPolygon * >( qgsgeometry_cast< const QgsMultiPolygon * >( tessellation.constGet() )->geometryN( i ) );
1108 if ( !mProfileCurveEngine->intersects( triangle ) )
1109 continue;
1110
1111 QString error;
1112 std::unique_ptr< QgsAbstractGeometry > intersection( mProfileCurveEngine->intersection( triangle, &error ) );
1113 if ( !intersection )
1114 continue;
1115
1116 processTriangleLineIntersect( triangle, intersection.get(), transformedParts, crossSectionParts );
1117 }
1118 };
1119
1120 QgsFeature feature;
1121 QgsFeatureIterator it = mSource->getFeatures( request );
1122 while ( it.nextFeature( feature ) )
1123 {
1124 if ( mFeedback->isCanceled() )
1125 return false;
1126
1127 if ( !mProfileCurveEngine->intersects( feature.geometry().constGet() ) )
1128 continue;
1129
1130 mExpressionContext.setFeature( feature );
1131
1132 const double offset = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::ZOffset, mExpressionContext, mOffset );
1133
1134 const QgsGeometry g = feature.geometry();
1135 QVector< QgsGeometry > transformedParts;
1136 QVector< QgsGeometry > crossSectionParts;
1137 if ( g.isMultipart() )
1138 {
1139 for ( auto it = g.const_parts_begin(); it != g.const_parts_end(); ++it )
1140 {
1141 if ( mFeedback->isCanceled() )
1142 break;
1143
1144 if ( !mProfileCurveEngine->intersects( *it ) )
1145 continue;
1146
1147 processPolygon( qgsgeometry_cast< const QgsCurvePolygon * >( *it ), transformedParts, crossSectionParts, offset );
1148 }
1149 }
1150 else
1151 {
1152 processPolygon( qgsgeometry_cast< const QgsCurvePolygon * >( g.constGet() ), transformedParts, crossSectionParts, offset );
1153 }
1154
1155 if ( mFeedback->isCanceled() )
1156 return false;
1157
1159 resultFeature.featureId = feature.id();
1160 resultFeature.geometry = transformedParts.size() > 1 ? QgsGeometry::collectGeometry( transformedParts ) : transformedParts.value( 0 );
1161 if ( !crossSectionParts.empty() )
1162 {
1163 QgsGeometry unioned = QgsGeometry::unaryUnion( crossSectionParts );
1164 if ( unioned.type() == Qgis::GeometryType::Line )
1165 unioned = unioned.mergeLines();
1166 resultFeature.crossSectionGeometry = unioned;
1167 }
1168 mResults->features[resultFeature.featureId].append( resultFeature );
1169 }
1170 return true;
1171}
1172
1173double QgsVectorLayerProfileGenerator::terrainHeight( double x, double y )
1174{
1175 if ( !mTerrainProvider )
1176 return std::numeric_limits<double>::quiet_NaN();
1177
1178 // transform feature point to terrain provider crs
1179 try
1180 {
1181 double dummyZ = 0;
1182 mTargetToTerrainProviderTransform.transformInPlace( x, y, dummyZ );
1183 }
1184 catch ( QgsCsException & )
1185 {
1186 return std::numeric_limits<double>::quiet_NaN();
1187 }
1188
1189 return mTerrainProvider->heightAt( x, y );
1190}
1191
1192double QgsVectorLayerProfileGenerator::featureZToHeight( double x, double y, double z, double offset )
1193{
1194 switch ( mClamping )
1195 {
1197 break;
1198
1201 {
1202 const double terrainZ = terrainHeight( x, y );
1203 if ( !std::isnan( terrainZ ) )
1204 {
1205 switch ( mClamping )
1206 {
1208 if ( std::isnan( z ) )
1209 z = terrainZ;
1210 else
1211 z += terrainZ;
1212 break;
1213
1215 z = terrainZ;
1216 break;
1217
1219 break;
1220 }
1221 }
1222 break;
1223 }
1224 }
1225
1226 return ( std::isnan( z ) ? 0 : z ) * mScale + offset;
1227}
1228
1229void QgsVectorLayerProfileGenerator::clampAltitudes( QgsLineString *lineString, const QgsPoint &centroid, double offset )
1230{
1231 for ( int i = 0; i < lineString->nCoordinates(); ++i )
1232 {
1233 if ( mFeedback->isCanceled() )
1234 break;
1235
1236 double terrainZ = 0;
1237 switch ( mClamping )
1238 {
1241 {
1242 QgsPointXY pt;
1243 switch ( mBinding )
1244 {
1246 pt.setX( lineString->xAt( i ) );
1247 pt.setY( lineString->yAt( i ) );
1248 break;
1249
1251 pt.set( centroid.x(), centroid.y() );
1252 break;
1253 }
1254
1255 terrainZ = terrainHeight( pt.x(), pt.y() );
1256 break;
1257 }
1258
1260 break;
1261 }
1262
1263 double geomZ = 0;
1264
1265 switch ( mClamping )
1266 {
1269 geomZ = lineString->zAt( i );
1270 break;
1271
1273 break;
1274 }
1275
1276 const double z = ( terrainZ + ( std::isnan( geomZ ) ? 0 : geomZ ) ) * mScale + offset;
1277 lineString->setZAt( i, z );
1278 }
1279}
1280
1281bool QgsVectorLayerProfileGenerator::clampAltitudes( QgsPolygon *polygon, double offset )
1282{
1283 if ( !polygon->is3D() )
1284 polygon->addZValue( 0 );
1285
1287 switch ( mBinding )
1288 {
1290 break;
1291
1293 centroid = polygon->centroid();
1294 break;
1295 }
1296
1297 QgsCurve *curve = const_cast<QgsCurve *>( polygon->exteriorRing() );
1298 QgsLineString *lineString = qgsgeometry_cast<QgsLineString *>( curve );
1299 if ( !lineString )
1300 return false;
1301
1302 clampAltitudes( lineString, centroid, offset );
1303
1304 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
1305 {
1306 if ( mFeedback->isCanceled() )
1307 break;
1308
1309 QgsCurve *curve = const_cast<QgsCurve *>( polygon->interiorRing( i ) );
1310 QgsLineString *lineString = qgsgeometry_cast<QgsLineString *>( curve );
1311 if ( !lineString )
1312 return false;
1313
1314 clampAltitudes( lineString, centroid, offset );
1315 }
1316 return true;
1317}
1318
@ 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:1895
@ Centroid
Clamp just centroid of feature.
@ Vertex
Clamp every vertex of feature.
@ 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)
@ Marker
Marker symbol.
@ Line
Line symbol.
@ Fill
Fill symbol.
Abstract base class for all geometries.
bool is3D() const SIP_HOLDGIL
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.
std::unique_ptr< QgsFillSymbol > mFillSymbol
std::unique_ptr< QgsLineSymbol > mLineSymbol
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< 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 SIP_THROW(QgsCsException)
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.
Definition: qgsexception.h:66
Curve polygon geometry type.
bool addZValue(double zValue=0) override
Adds a z-dimension to the geometry, initialized to a preset value.
const QgsCurve * interiorRing(int i) const SIP_HOLDGIL
Retrieves an interior ring from the curve polygon.
const QgsCurve * exteriorRing() const SIP_HOLDGIL
Returns the curve polygon's exterior ring.
QgsAbstractGeometry * segmentize(double tolerance=M_PI_2/90, SegmentationToleranceType toleranceType=MaximumAngle) const override
Returns a geometry without curves.
int numInteriorRings() const SIP_HOLDGIL
Returns the number of interior rings contained with the curve polygon.
Abstract base class for curved geometry type.
Definition: qgscurve.h:36
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:203
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)
This class wraps a request for features to a vector layer (or directly its vector data provider).
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:56
QgsGeometry geometry
Definition: qgsfeature.h:67
Q_GADGET QgsFeatureId id
Definition: qgsfeature.h:64
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition: qgsfeedback.h:45
A fill symbol type, for rendering Polygon and MultiPolygon geometries.
Definition: qgsfillsymbol.h:30
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:164
const QgsAbstractGeometry * constGet() const SIP_HOLDGIL
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
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) SIP_THROW(QgsCsException)
Transforms this geometry as described by the coordinate transform ct.
static QgsGeometry fromRect(const QgsRectangle &rect) SIP_HOLDGIL
Creates a new geometry from a QgsRectangle.
bool isMultipart() const SIP_HOLDGIL
Returns true if WKB of the geometry is of WKBMulti* type.
Qgis::GeometryType type
Definition: qgsgeometry.h:167
static QgsGeometry fromPointXY(const QgsPointXY &point) SIP_HOLDGIL
Creates a new geometry from a QgsPointXY object.
QgsGeometry mergeLines() const
Merges any connected lines in a LineString/MultiLineString geometry and converts them to single line ...
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.
Does vector analysis using the geos library and handles import, export, exception handling*.
Definition: qgsgeos.h:99
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:45
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.
int nCoordinates() const override SIP_HOLDGIL
Returns the number of nodes contained in the geometry.
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.
Definition: qgslinesymbol.h:30
A marker symbol type, for rendering Point and MultiPoint geometries.
Multi line string geometry collection.
Multi point geometry collection.
Definition: qgsmultipoint.h:30
Multi polygon geometry collection.
A class to represent a 2D point.
Definition: qgspointxy.h:59
void set(double x, double y) SIP_HOLDGIL
Sets the x and y value of the point.
Definition: qgspointxy.h:139
void setX(double x) SIP_HOLDGIL
Sets the x value of the point.
Definition: qgspointxy.h:122
double y
Definition: qgspointxy.h:63
Q_GADGET double x
Definition: qgspointxy.h:62
void setY(double y) SIP_HOLDGIL
Sets the y value of the point.
Definition: qgspointxy.h:132
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 y
Definition: qgspoint.h:53
Polygon geometry type.
Definition: qgspolygon.h:34
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 SIP_HOLDGIL
Returns the elevation of the point.
double distance() const SIP_HOLDGIL
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.
QSet< QString > referencedFields(const QgsExpressionContext &context=QgsExpressionContext(), bool ignoreContext=false) const override
Returns the set of any fields referenced by the active properties from the collection.
bool prepare(const QgsExpressionContext &context=QgsExpressionContext()) const override
Prepares the collection against a specified expression context.
T lower() const
Returns the lower bound of the range.
Definition: qgsrange.h:66
T upper() const
Returns the upper bound of the range.
Definition: qgsrange.h:73
A rectangle specified with double values.
Definition: qgsrectangle.h:42
bool intersects(const QgsRectangle &rect) const SIP_HOLDGIL
Returns true when rectangle intersects with other rectangle.
Definition: qgsrectangle.h:349
double xMaximum() const SIP_HOLDGIL
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:183
double xMinimum() const SIP_HOLDGIL
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:188
double yMinimum() const SIP_HOLDGIL
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:198
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:93
void setColor(const QColor &color) const
Sets the color for the symbol.
Definition: qgssymbol.cpp:898
qreal opacity() const
Returns the opacity for the symbol.
Definition: qgssymbol.h:495
QColor color() const
Returns the symbol's color.
Definition: qgssymbol.cpp:908
Qgis::SymbolType type() const
Returns the symbol's type.
Definition: qgssymbol.h:152
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.
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) SIP_HOLDGIL
Returns the geometry type for a WKB type, e.g., both MultiPolygon and CurvePolygon would have a Polyg...
Definition: qgswkbtypes.h:865
CORE_EXPORT QgsMeshVertex centroid(const QgsMeshFace &face, const QVector< QgsMeshVertex > &vertices)
Returns the centroid of the face.
#define BUILTIN_UNREACHABLE
Definition: qgis.h:4180
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:3509
QSet< QgsFeatureId > QgsFeatureIds
Definition: qgsfeatureid.h:37
qint64 QgsFeatureId
64 bit feature ids negative numbers are used for uncommitted/newly added features
Definition: qgsfeatureid.h:28
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
const QgsCoordinateReferenceSystem & crs
QgsGeometry crossSectionGeometry
Cross section distance vs height geometry for feature.
QgsGeometry geometry
Feature's geometry with any terrain height adjustment and extrusion applied.