QGIS API Documentation 3.31.0-Master (9f23a2c1dc)
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 if ( mLayer )
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 mElevationLimit = qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->elevationLimit();
648
649 mLineSymbol.reset( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->profileLineSymbol()->clone() );
650 mFillSymbol.reset( qgis::down_cast< QgsVectorLayerElevationProperties * >( layer->elevationProperties() )->profileFillSymbol()->clone() );
651}
652
654{
655 return mId;
656}
657
659
661{
662 if ( !mProfileCurve || mFeedback->isCanceled() )
663 return false;
664
665 // we need to transform the profile curve to the vector's CRS
666 mTransformedCurve.reset( mProfileCurve->clone() );
667 mLayerToTargetTransform = QgsCoordinateTransform( mSourceCrs, mTargetCrs, mTransformContext );
668 if ( mTerrainProvider )
669 mTargetToTerrainProviderTransform = QgsCoordinateTransform( mTargetCrs, mTerrainProvider->crs(), mTransformContext );
670
671 try
672 {
673 mTransformedCurve->transform( mLayerToTargetTransform, Qgis::TransformDirection::Reverse );
674 }
675 catch ( QgsCsException & )
676 {
677 QgsDebugError( QStringLiteral( "Error transforming profile line to vector CRS" ) );
678 return false;
679 }
680
681 const QgsRectangle profileCurveBoundingBox = mTransformedCurve->boundingBox();
682 if ( !profileCurveBoundingBox.intersects( mExtent ) )
683 return false;
684
685 if ( mFeedback->isCanceled() )
686 return false;
687
688 mResults = std::make_unique< QgsVectorLayerProfileResults >();
689 mResults->mLayer = mLayer;
690 mResults->copyPropertiesFromGenerator( this );
691
692 mProfileCurveEngine.reset( new QgsGeos( mProfileCurve.get() ) );
693 mProfileCurveEngine->prepareGeometry();
694
695 mDataDefinedProperties.prepare( mExpressionContext );
696
697 if ( mFeedback->isCanceled() )
698 return false;
699
700 switch ( QgsWkbTypes::geometryType( mWkbType ) )
701 {
702 case Qgis::GeometryType::Point:
703 if ( !generateProfileForPoints() )
704 return false;
705 break;
706
707 case Qgis::GeometryType::Line:
708 if ( !generateProfileForLines() )
709 return false;
710 break;
711
712 case Qgis::GeometryType::Polygon:
713 if ( !generateProfileForPolygons() )
714 return false;
715 break;
716
717 case Qgis::GeometryType::Unknown:
718 case Qgis::GeometryType::Null:
719 return false;
720 }
721
722 return true;
723}
724
726{
727 return mResults.release();
728}
729
731{
732 return mFeedback.get();
733}
734
735bool QgsVectorLayerProfileGenerator::generateProfileForPoints()
736{
737 // get features from layer
738 QgsFeatureRequest request;
739 request.setDestinationCrs( mTargetCrs, mTransformContext );
740 request.setDistanceWithin( QgsGeometry( mProfileCurve->clone() ), mTolerance );
741 request.setSubsetOfAttributes( mDataDefinedProperties.referencedFields( mExpressionContext ), mFields );
742 request.setFeedback( mFeedback.get() );
743
744 // our feature request is using the optimised distance within check (allowing use of spatial index)
745 // BUT this will also include points which are within the tolerance distance before/after the end of line.
746 // So we also need to double check that they fall within the flat buffered curve too.
747 std::unique_ptr< QgsAbstractGeometry > bufferedCurve( mProfileCurveEngine->buffer( mTolerance, 8, Qgis::EndCapStyle::Flat, Qgis::JoinStyle::Round, 2 ) );
748 QgsGeos bufferedCurveEngine( bufferedCurve.get() );
749 bufferedCurveEngine.prepareGeometry();
750
751 auto processPoint = [this, &bufferedCurveEngine]( const QgsFeature & feature, const QgsPoint * point )
752 {
753 if ( !bufferedCurveEngine.intersects( point ) )
754 return;
755
756 const double offset = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::ZOffset, mExpressionContext, mOffset );
757
758 const double height = featureZToHeight( point->x(), point->y(), point->z(), offset );
759 mResults->mRawPoints.append( QgsPoint( point->x(), point->y(), height ) );
760 mResults->minZ = std::min( mResults->minZ, height );
761 mResults->maxZ = std::max( mResults->maxZ, height );
762
763 QString lastError;
764 const double distance = mProfileCurveEngine->lineLocatePoint( *point, &lastError );
765 mResults->mDistanceToHeightMap.insert( distance, height );
766
768 resultFeature.featureId = feature.id();
769 if ( mExtrusionEnabled )
770 {
771 const double extrusion = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::ExtrusionHeight, mExpressionContext, mExtrusionHeight );
772
773 resultFeature.geometry = QgsGeometry( new QgsLineString( QgsPoint( point->x(), point->y(), height ),
774 QgsPoint( point->x(), point->y(), height + extrusion ) ) );
775 resultFeature.crossSectionGeometry = QgsGeometry( new QgsLineString( QgsPoint( distance, height ),
776 QgsPoint( distance, height + extrusion ) ) );
777 mResults->minZ = std::min( mResults->minZ, height + extrusion );
778 mResults->maxZ = std::max( mResults->maxZ, height + extrusion );
779 }
780 else
781 {
782 resultFeature.geometry = QgsGeometry( new QgsPoint( point->x(), point->y(), height ) );
783 resultFeature.crossSectionGeometry = QgsGeometry( new QgsPoint( distance, height ) );
784 }
785 mResults->features[resultFeature.featureId].append( resultFeature );
786 };
787
788 QgsFeature feature;
789 QgsFeatureIterator it = mSource->getFeatures( request );
790 while ( it.nextFeature( feature ) )
791 {
792 if ( mFeedback->isCanceled() )
793 return false;
794
795 mExpressionContext.setFeature( feature );
796
797 const QgsGeometry g = feature.geometry();
798 if ( g.isMultipart() )
799 {
800 for ( auto it = g.const_parts_begin(); it != g.const_parts_end(); ++it )
801 {
802 processPoint( feature, qgsgeometry_cast< const QgsPoint * >( *it ) );
803 }
804 }
805 else
806 {
807 processPoint( feature, qgsgeometry_cast< const QgsPoint * >( g.constGet() ) );
808 }
809 }
810 return true;
811}
812
813bool QgsVectorLayerProfileGenerator::generateProfileForLines()
814{
815 // get features from layer
816 QgsFeatureRequest request;
817 request.setDestinationCrs( mTargetCrs, mTransformContext );
818 request.setFilterRect( mProfileCurve->boundingBox() );
819 request.setSubsetOfAttributes( mDataDefinedProperties.referencedFields( mExpressionContext ), mFields );
820 request.setFeedback( mFeedback.get() );
821
822 auto processCurve = [this]( const QgsFeature & feature, const QgsCurve * curve )
823 {
824 QString error;
825 std::unique_ptr< QgsAbstractGeometry > intersection( mProfileCurveEngine->intersection( curve, &error ) );
826 if ( !intersection )
827 return;
828
829 if ( mFeedback->isCanceled() )
830 return;
831
832 QgsGeos curveGeos( curve );
833 curveGeos.prepareGeometry();
834
835 if ( mFeedback->isCanceled() )
836 return;
837
838 for ( auto it = intersection->const_parts_begin(); it != intersection->const_parts_end(); ++it )
839 {
840 if ( mFeedback->isCanceled() )
841 return;
842
843 if ( const QgsPoint *intersectionPoint = qgsgeometry_cast< const QgsPoint * >( *it ) )
844 {
845 // unfortunately we need to do some work to interpolate the z value for the line -- GEOS doesn't give us this
846 const double distance = curveGeos.lineLocatePoint( *intersectionPoint, &error );
847 std::unique_ptr< QgsPoint > interpolatedPoint( curve->interpolatePoint( distance ) );
848
849 const double offset = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::ZOffset, mExpressionContext, mOffset );
850
851 const double height = featureZToHeight( interpolatedPoint->x(), interpolatedPoint->y(), interpolatedPoint->z(), offset );
852 mResults->mRawPoints.append( QgsPoint( interpolatedPoint->x(), interpolatedPoint->y(), height ) );
853 mResults->minZ = std::min( mResults->minZ, height );
854 mResults->maxZ = std::max( mResults->maxZ, height );
855
856 const double distanceAlongProfileCurve = mProfileCurveEngine->lineLocatePoint( *interpolatedPoint, &error );
857 mResults->mDistanceToHeightMap.insert( distanceAlongProfileCurve, height );
858
860 resultFeature.featureId = feature.id();
861 if ( mExtrusionEnabled )
862 {
863 const double extrusion = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::ExtrusionHeight, mExpressionContext, mExtrusionHeight );
864
865 resultFeature.geometry = QgsGeometry( new QgsLineString( QgsPoint( interpolatedPoint->x(), interpolatedPoint->y(), height ),
866 QgsPoint( interpolatedPoint->x(), interpolatedPoint->y(), height + extrusion ) ) );
867 resultFeature.crossSectionGeometry = QgsGeometry( new QgsLineString( QgsPoint( distanceAlongProfileCurve, height ),
868 QgsPoint( distanceAlongProfileCurve, height + extrusion ) ) );
869 mResults->minZ = std::min( mResults->minZ, height + extrusion );
870 mResults->maxZ = std::max( mResults->maxZ, height + extrusion );
871 }
872 else
873 {
874 resultFeature.geometry = QgsGeometry( new QgsPoint( interpolatedPoint->x(), interpolatedPoint->y(), height ) );
875 resultFeature.crossSectionGeometry = QgsGeometry( new QgsPoint( distanceAlongProfileCurve, height ) );
876 }
877 mResults->features[resultFeature.featureId].append( resultFeature );
878 }
879 }
880 };
881
882 QgsFeature feature;
883 QgsFeatureIterator it = mSource->getFeatures( request );
884 while ( it.nextFeature( feature ) )
885 {
886 if ( mFeedback->isCanceled() )
887 return false;
888
889 if ( !mProfileCurveEngine->intersects( feature.geometry().constGet() ) )
890 continue;
891
892 mExpressionContext.setFeature( feature );
893
894 const QgsGeometry g = feature.geometry();
895 if ( g.isMultipart() )
896 {
897 for ( auto it = g.const_parts_begin(); it != g.const_parts_end(); ++it )
898 {
899 if ( !mProfileCurveEngine->intersects( *it ) )
900 continue;
901
902 processCurve( feature, qgsgeometry_cast< const QgsCurve * >( *it ) );
903 }
904 }
905 else
906 {
907 processCurve( feature, qgsgeometry_cast< const QgsCurve * >( g.constGet() ) );
908 }
909 }
910 return true;
911}
912
913bool QgsVectorLayerProfileGenerator::generateProfileForPolygons()
914{
915 // get features from layer
916 QgsFeatureRequest request;
917 request.setDestinationCrs( mTargetCrs, mTransformContext );
918 request.setFilterRect( mProfileCurve->boundingBox() );
919 request.setSubsetOfAttributes( mDataDefinedProperties.referencedFields( mExpressionContext ), mFields );
920 request.setFeedback( mFeedback.get() );
921
922 auto interpolatePointOnTriangle = []( const QgsPolygon * triangle, double x, double y ) -> QgsPoint
923 {
924 QgsPoint p1, p2, p3;
926 triangle->exteriorRing()->pointAt( 0, p1, vt );
927 triangle->exteriorRing()->pointAt( 1, p2, vt );
928 triangle->exteriorRing()->pointAt( 2, p3, vt );
929 const double z = QgsMeshLayerUtils::interpolateFromVerticesData( p1, p2, p3, p1.z(), p2.z(), p3.z(), QgsPointXY( x, y ) );
930 return QgsPoint( x, y, z );
931 };
932
933 std::function< void( const QgsPolygon *triangle, const QgsAbstractGeometry *intersect, QVector< QgsGeometry > &, QVector< QgsGeometry > & ) > processTriangleLineIntersect;
934 processTriangleLineIntersect = [this, &interpolatePointOnTriangle, &processTriangleLineIntersect]( const QgsPolygon * triangle, const QgsAbstractGeometry * intersect, QVector< QgsGeometry > &transformedParts, QVector< QgsGeometry > &crossSectionParts )
935 {
936 // intersect may be a (multi)point or (multi)linestring
937 switch ( QgsWkbTypes::geometryType( intersect->wkbType() ) )
938 {
939 case Qgis::GeometryType::Point:
940 if ( const QgsMultiPoint *mp = qgsgeometry_cast< const QgsMultiPoint * >( intersect ) )
941 {
942 const int numPoint = mp->numGeometries();
943 for ( int i = 0; i < numPoint; ++i )
944 {
945 processTriangleLineIntersect( triangle, mp->geometryN( i ), transformedParts, crossSectionParts );
946 }
947 }
948 else if ( const QgsPoint *p = qgsgeometry_cast< const QgsPoint * >( intersect ) )
949 {
950 const QgsPoint interpolatedPoint = interpolatePointOnTriangle( triangle, p->x(), p->y() );
951 mResults->mRawPoints.append( interpolatedPoint );
952 mResults->minZ = std::min( mResults->minZ, interpolatedPoint.z() );
953 mResults->maxZ = std::max( mResults->maxZ, interpolatedPoint.z() );
954
955 QString lastError;
956 const double distance = mProfileCurveEngine->lineLocatePoint( *p, &lastError );
957 mResults->mDistanceToHeightMap.insert( distance, interpolatedPoint.z() );
958
959 if ( mExtrusionEnabled )
960 {
961 const double extrusion = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::ExtrusionHeight, mExpressionContext, mExtrusionHeight );
962
963 transformedParts.append( QgsGeometry( new QgsLineString( interpolatedPoint,
964 QgsPoint( interpolatedPoint.x(), interpolatedPoint.y(), interpolatedPoint.z() + extrusion ) ) ) );
965 crossSectionParts.append( QgsGeometry( new QgsLineString( QgsPoint( distance, interpolatedPoint.z() ),
966 QgsPoint( distance, interpolatedPoint.z() + extrusion ) ) ) );
967 mResults->minZ = std::min( mResults->minZ, interpolatedPoint.z() + extrusion );
968 mResults->maxZ = std::max( mResults->maxZ, interpolatedPoint.z() + extrusion );
969 }
970 else
971 {
972 transformedParts.append( QgsGeometry( new QgsPoint( interpolatedPoint ) ) );
973 crossSectionParts.append( QgsGeometry( new QgsPoint( distance, interpolatedPoint.z() ) ) );
974 }
975 }
976 break;
977 case Qgis::GeometryType::Line:
978 if ( const QgsMultiLineString *ml = qgsgeometry_cast< const QgsMultiLineString * >( intersect ) )
979 {
980 const int numLines = ml->numGeometries();
981 for ( int i = 0; i < numLines; ++i )
982 {
983 processTriangleLineIntersect( triangle, ml->geometryN( i ), transformedParts, crossSectionParts );
984 }
985 }
986 else if ( const QgsLineString *ls = qgsgeometry_cast< const QgsLineString * >( intersect ) )
987 {
988 const int numPoints = ls->numPoints();
989 QVector< double > newX;
990 newX.resize( numPoints );
991 QVector< double > newY;
992 newY.resize( numPoints );
993 QVector< double > newZ;
994 newZ.resize( numPoints );
995 QVector< double > newDistance;
996 newDistance.resize( numPoints );
997
998 const double *inX = ls->xData();
999 const double *inY = ls->yData();
1000 double *outX = newX.data();
1001 double *outY = newY.data();
1002 double *outZ = newZ.data();
1003 double *outDistance = newDistance.data();
1004
1005 QVector< double > extrudedZ;
1006 double *extZOut = nullptr;
1007 double extrusion = 0;
1008 if ( mExtrusionEnabled )
1009 {
1010 extrudedZ.resize( numPoints );
1011 extZOut = extrudedZ.data();
1012
1013 extrusion = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::ExtrusionHeight, mExpressionContext, mExtrusionHeight );
1014 }
1015
1016 QString lastError;
1017 for ( int i = 0 ; i < numPoints; ++i )
1018 {
1019 double x = *inX++;
1020 double y = *inY++;
1021
1022 QgsPoint interpolatedPoint = interpolatePointOnTriangle( triangle, x, y );
1023 *outX++ = x;
1024 *outY++ = y;
1025 *outZ++ = interpolatedPoint.z();
1026 if ( extZOut )
1027 *extZOut++ = interpolatedPoint.z() + extrusion;
1028
1029 mResults->mRawPoints.append( interpolatedPoint );
1030 mResults->minZ = std::min( mResults->minZ, interpolatedPoint.z() );
1031 mResults->maxZ = std::max( mResults->maxZ, interpolatedPoint.z() );
1032 if ( mExtrusionEnabled )
1033 {
1034 mResults->minZ = std::min( mResults->minZ, interpolatedPoint.z() + extrusion );
1035 mResults->maxZ = std::max( mResults->maxZ, interpolatedPoint.z() + extrusion );
1036 }
1037
1038 const double distance = mProfileCurveEngine->lineLocatePoint( interpolatedPoint, &lastError );
1039 *outDistance++ = distance;
1040
1041 mResults->mDistanceToHeightMap.insert( distance, interpolatedPoint.z() );
1042 }
1043
1044 if ( mExtrusionEnabled )
1045 {
1046 std::unique_ptr< QgsLineString > ring = std::make_unique< QgsLineString >( newX, newY, newZ );
1047 std::unique_ptr< QgsLineString > extrudedRing = std::make_unique< QgsLineString >( newX, newY, extrudedZ );
1048 std::unique_ptr< QgsLineString > reversedExtrusion( extrudedRing->reversed() );
1049 ring->append( reversedExtrusion.get() );
1050 ring->close();
1051 transformedParts.append( QgsGeometry( new QgsPolygon( ring.release() ) ) );
1052
1053
1054 std::unique_ptr< QgsLineString > distanceVHeightRing = std::make_unique< QgsLineString >( newDistance, newZ );
1055 std::unique_ptr< QgsLineString > extrudedDistanceVHeightRing = std::make_unique< QgsLineString >( newDistance, extrudedZ );
1056 std::unique_ptr< QgsLineString > reversedDistanceVHeightExtrusion( extrudedDistanceVHeightRing->reversed() );
1057 distanceVHeightRing->append( reversedDistanceVHeightExtrusion.get() );
1058 distanceVHeightRing->close();
1059 crossSectionParts.append( QgsGeometry( new QgsPolygon( distanceVHeightRing.release() ) ) );
1060 }
1061 else
1062 {
1063 transformedParts.append( QgsGeometry( new QgsLineString( newX, newY, newZ ) ) );
1064 crossSectionParts.append( QgsGeometry( new QgsLineString( newDistance, newZ ) ) );
1065 }
1066 }
1067 break;
1068
1069 case Qgis::GeometryType::Polygon:
1070 case Qgis::GeometryType::Unknown:
1071 case Qgis::GeometryType::Null:
1072 return;
1073 }
1074 };
1075
1076 auto processPolygon = [this, &processTriangleLineIntersect]( const QgsCurvePolygon * polygon, QVector< QgsGeometry > &transformedParts, QVector< QgsGeometry > &crossSectionParts, double offset )
1077 {
1078 std::unique_ptr< QgsPolygon > clampedPolygon;
1079 if ( const QgsPolygon *p = qgsgeometry_cast< const QgsPolygon * >( polygon ) )
1080 {
1081 clampedPolygon.reset( p->clone() );
1082 }
1083 else
1084 {
1085 clampedPolygon.reset( qgsgeometry_cast< QgsPolygon * >( polygon->segmentize() ) );
1086 }
1087 clampAltitudes( clampedPolygon.get(), offset );
1088
1089 if ( mFeedback->isCanceled() )
1090 return;
1091
1092 const QgsRectangle bounds = clampedPolygon->boundingBox();
1093 QgsTessellator t( bounds, false, false, false, false );
1094 t.addPolygon( *clampedPolygon, 0 );
1095
1096 QgsGeometry tessellation( t.asMultiPolygon() );
1097 if ( mFeedback->isCanceled() )
1098 return;
1099
1100 tessellation.translate( bounds.xMinimum(), bounds.yMinimum() );
1101
1102 // iterate through the tessellation, finding triangles which intersect the line
1103 const int numTriangles = qgsgeometry_cast< const QgsMultiPolygon * >( tessellation.constGet() )->numGeometries();
1104 for ( int i = 0; i < numTriangles; ++i )
1105 {
1106 if ( mFeedback->isCanceled() )
1107 return;
1108
1109 const QgsPolygon *triangle = qgsgeometry_cast< const QgsPolygon * >( qgsgeometry_cast< const QgsMultiPolygon * >( tessellation.constGet() )->geometryN( i ) );
1110 if ( !mProfileCurveEngine->intersects( triangle ) )
1111 continue;
1112
1113 QString error;
1114 std::unique_ptr< QgsAbstractGeometry > intersection( mProfileCurveEngine->intersection( triangle, &error ) );
1115 if ( !intersection )
1116 continue;
1117
1118 processTriangleLineIntersect( triangle, intersection.get(), transformedParts, crossSectionParts );
1119 }
1120 };
1121
1122 QgsFeature feature;
1123 QgsFeatureIterator it = mSource->getFeatures( request );
1124 while ( it.nextFeature( feature ) )
1125 {
1126 if ( mFeedback->isCanceled() )
1127 return false;
1128
1129 if ( !mProfileCurveEngine->intersects( feature.geometry().constGet() ) )
1130 continue;
1131
1132 mExpressionContext.setFeature( feature );
1133
1134 const double offset = mDataDefinedProperties.valueAsDouble( QgsMapLayerElevationProperties::ZOffset, mExpressionContext, mOffset );
1135
1136 const QgsGeometry g = feature.geometry();
1137 QVector< QgsGeometry > transformedParts;
1138 QVector< QgsGeometry > crossSectionParts;
1139 if ( g.isMultipart() )
1140 {
1141 for ( auto it = g.const_parts_begin(); it != g.const_parts_end(); ++it )
1142 {
1143 if ( mFeedback->isCanceled() )
1144 break;
1145
1146 if ( !mProfileCurveEngine->intersects( *it ) )
1147 continue;
1148
1149 processPolygon( qgsgeometry_cast< const QgsCurvePolygon * >( *it ), transformedParts, crossSectionParts, offset );
1150 }
1151 }
1152 else
1153 {
1154 processPolygon( qgsgeometry_cast< const QgsCurvePolygon * >( g.constGet() ), transformedParts, crossSectionParts, offset );
1155 }
1156
1157 if ( mFeedback->isCanceled() )
1158 return false;
1159
1161 resultFeature.featureId = feature.id();
1162 resultFeature.geometry = transformedParts.size() > 1 ? QgsGeometry::collectGeometry( transformedParts ) : transformedParts.value( 0 );
1163 if ( !crossSectionParts.empty() )
1164 {
1165 QgsGeometry unioned = QgsGeometry::unaryUnion( crossSectionParts );
1166 if ( unioned.type() == Qgis::GeometryType::Line )
1167 unioned = unioned.mergeLines();
1168 resultFeature.crossSectionGeometry = unioned;
1169 }
1170 mResults->features[resultFeature.featureId].append( resultFeature );
1171 }
1172 return true;
1173}
1174
1175double QgsVectorLayerProfileGenerator::terrainHeight( double x, double y )
1176{
1177 if ( !mTerrainProvider )
1178 return std::numeric_limits<double>::quiet_NaN();
1179
1180 // transform feature point to terrain provider crs
1181 try
1182 {
1183 double dummyZ = 0;
1184 mTargetToTerrainProviderTransform.transformInPlace( x, y, dummyZ );
1185 }
1186 catch ( QgsCsException & )
1187 {
1188 return std::numeric_limits<double>::quiet_NaN();
1189 }
1190
1191 return mTerrainProvider->heightAt( x, y );
1192}
1193
1194double QgsVectorLayerProfileGenerator::featureZToHeight( double x, double y, double z, double offset )
1195{
1196 switch ( mClamping )
1197 {
1199 break;
1200
1203 {
1204 const double terrainZ = terrainHeight( x, y );
1205 if ( !std::isnan( terrainZ ) )
1206 {
1207 switch ( mClamping )
1208 {
1210 if ( std::isnan( z ) )
1211 z = terrainZ;
1212 else
1213 z += terrainZ;
1214 break;
1215
1217 z = terrainZ;
1218 break;
1219
1221 break;
1222 }
1223 }
1224 break;
1225 }
1226 }
1227
1228 return ( std::isnan( z ) ? 0 : z ) * mScale + offset;
1229}
1230
1231void QgsVectorLayerProfileGenerator::clampAltitudes( QgsLineString *lineString, const QgsPoint &centroid, double offset )
1232{
1233 for ( int i = 0; i < lineString->nCoordinates(); ++i )
1234 {
1235 if ( mFeedback->isCanceled() )
1236 break;
1237
1238 double terrainZ = 0;
1239 switch ( mClamping )
1240 {
1243 {
1244 QgsPointXY pt;
1245 switch ( mBinding )
1246 {
1248 pt.setX( lineString->xAt( i ) );
1249 pt.setY( lineString->yAt( i ) );
1250 break;
1251
1253 pt.set( centroid.x(), centroid.y() );
1254 break;
1255 }
1256
1257 terrainZ = terrainHeight( pt.x(), pt.y() );
1258 break;
1259 }
1260
1262 break;
1263 }
1264
1265 double geomZ = 0;
1266
1267 switch ( mClamping )
1268 {
1271 geomZ = lineString->zAt( i );
1272 break;
1273
1275 break;
1276 }
1277
1278 const double z = ( terrainZ + ( std::isnan( geomZ ) ? 0 : geomZ ) ) * mScale + offset;
1279 lineString->setZAt( i, z );
1280 }
1281}
1282
1283bool QgsVectorLayerProfileGenerator::clampAltitudes( QgsPolygon *polygon, double offset )
1284{
1285 if ( !polygon->is3D() )
1286 polygon->addZValue( 0 );
1287
1289 switch ( mBinding )
1290 {
1292 break;
1293
1295 centroid = polygon->centroid();
1296 break;
1297 }
1298
1299 QgsCurve *curve = const_cast<QgsCurve *>( polygon->exteriorRing() );
1300 QgsLineString *lineString = qgsgeometry_cast<QgsLineString *>( curve );
1301 if ( !lineString )
1302 return false;
1303
1304 clampAltitudes( lineString, centroid, offset );
1305
1306 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
1307 {
1308 if ( mFeedback->isCanceled() )
1309 break;
1310
1311 QgsCurve *curve = const_cast<QgsCurve *>( polygon->interiorRing( i ) );
1312 QgsLineString *lineString = qgsgeometry_cast<QgsLineString *>( curve );
1313 if ( !lineString )
1314 return false;
1315
1316 clampAltitudes( lineString, centroid, offset );
1317 }
1318 return true;
1319}
1320
@ 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:2010
@ 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:67
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:901
qreal opacity() const
Returns the opacity for the symbol.
Definition: qgssymbol.h:495
QColor color() const
Returns the symbol's color.
Definition: qgssymbol.cpp:911
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:4574
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:3903
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 QgsDebugError(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.