QGIS API Documentation 3.29.0-Master (8c80f25a4f)
qgsrasterlayerprofilegenerator.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsrasterlayerprofilegenerator.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 "qgsrasterlayer.h"
22#include "qgsrasteriterator.h"
23#include "qgsgeometryengine.h"
24#include "qgsgeos.h"
25#include "qgslinesymbol.h"
26#include "qgsgeometryutils.h"
27#include "qgsprofilepoint.h"
28#include "qgsfillsymbol.h"
29
30#include <QPolygonF>
31
32//
33// QgsRasterLayerProfileResults
34//
35
37{
38 return QStringLiteral( "raster" );
39}
40
41QVector<QgsProfileIdentifyResults> QgsRasterLayerProfileResults::identify( const QgsProfilePoint &point, const QgsProfileIdentifyContext &context )
42{
43 const QVector<QgsProfileIdentifyResults> noLayerResults = QgsAbstractProfileSurfaceResults::identify( point, context );
44
45 // we have to make a new list, with the correct layer reference set
46 QVector<QgsProfileIdentifyResults> res;
47 res.reserve( noLayerResults.size() );
48 for ( const QgsProfileIdentifyResults &result : noLayerResults )
49 {
50 res.append( QgsProfileIdentifyResults( mLayer, result.results() ) );
51 }
52 return res;
53}
54
55
56
57//
58// QgsRasterLayerProfileGenerator
59//
60
62 : mId( layer->id() )
63 , mFeedback( std::make_unique< QgsRasterBlockFeedback >() )
64 , mProfileCurve( request.profileCurve() ? request.profileCurve()->clone() : nullptr )
65 , mSourceCrs( layer->crs() )
66 , mTargetCrs( request.crs() )
67 , mTransformContext( request.transformContext() )
68 , mOffset( layer->elevationProperties()->zOffset() )
69 , mScale( layer->elevationProperties()->zScale() )
70 , mLayer( layer )
71 , mBand( qgis::down_cast< QgsRasterLayerElevationProperties * >( layer->elevationProperties() )->bandNumber() )
72 , mRasterUnitsPerPixelX( layer->rasterUnitsPerPixelX() )
73 , mRasterUnitsPerPixelY( layer->rasterUnitsPerPixelY() )
74 , mStepDistance( request.stepDistance() )
75{
76 mRasterProvider.reset( layer->dataProvider()->clone() );
77
78 mSymbology = qgis::down_cast< QgsRasterLayerElevationProperties * >( layer->elevationProperties() )->profileSymbology();
79 mLineSymbol.reset( qgis::down_cast< QgsRasterLayerElevationProperties * >( layer->elevationProperties() )->profileLineSymbol()->clone() );
80 mFillSymbol.reset( qgis::down_cast< QgsRasterLayerElevationProperties * >( layer->elevationProperties() )->profileFillSymbol()->clone() );
81}
82
84{
85 return mId;
86}
87
88Qgis::ProfileGeneratorFlags QgsRasterLayerProfileGenerator::flags() const
89{
91}
92
94
96{
97 if ( !mProfileCurve || mFeedback->isCanceled() )
98 return false;
99
100 const double startDistanceOffset = std::max( !context.distanceRange().isInfinite() ? context.distanceRange().lower() : 0, 0.0 );
101 const double endDistance = context.distanceRange().upper();
102
103 std::unique_ptr< QgsCurve > trimmedCurve;
104 QgsCurve *sourceCurve = nullptr;
105 if ( startDistanceOffset > 0 || endDistance < mProfileCurve->length() )
106 {
107 trimmedCurve.reset( mProfileCurve->curveSubstring( startDistanceOffset, endDistance ) );
108 sourceCurve = trimmedCurve.get();
109 }
110 else
111 {
112 sourceCurve = mProfileCurve.get();
113 }
114
115 // we need to transform the profile curve to the raster's CRS
116 std::unique_ptr< QgsCurve > transformedCurve( sourceCurve->clone() );
117 const QgsCoordinateTransform rasterToTargetTransform( mSourceCrs, mTargetCrs, mTransformContext );
118 try
119 {
120 transformedCurve->transform( rasterToTargetTransform, Qgis::TransformDirection::Reverse );
121 }
122 catch ( QgsCsException & )
123 {
124 QgsDebugMsg( QStringLiteral( "Error transforming profile line to raster CRS" ) );
125 return false;
126 }
127
128 if ( mFeedback->isCanceled() )
129 return false;
130
131 const QgsRectangle profileCurveBoundingBox = transformedCurve->boundingBox();
132 if ( !profileCurveBoundingBox.intersects( mRasterProvider->extent() ) )
133 return false;
134
135 if ( mFeedback->isCanceled() )
136 return false;
137
138 mResults = std::make_unique< QgsRasterLayerProfileResults >();
139 mResults->mLayer = mLayer;
140 mResults->copyPropertiesFromGenerator( this );
141
142 std::unique_ptr< QgsGeometryEngine > curveEngine( QgsGeometry::createGeometryEngine( transformedCurve.get() ) );
143 curveEngine->prepareGeometry();
144
145 if ( mFeedback->isCanceled() )
146 return false;
147
148 double stepDistance = mStepDistance;
149 if ( !std::isnan( context.maximumErrorMapUnits() ) )
150 {
151 // convert the maximum error in curve units to a step distance
152 // TODO -- there's no point in this being << pixel size!
153 if ( std::isnan( stepDistance ) || context.maximumErrorMapUnits() > stepDistance )
154 {
155 stepDistance = context.maximumErrorMapUnits();
156 }
157 }
158
159 QSet< QgsPointXY > profilePoints;
160 if ( !std::isnan( stepDistance ) )
161 {
162 // if specific step distance specified, use this to generate points along the curve
163 QgsGeometry densifiedCurve( sourceCurve->clone() );
164 densifiedCurve = densifiedCurve.densifyByDistance( stepDistance );
165 densifiedCurve.transform( rasterToTargetTransform, Qgis::TransformDirection::Reverse );
166 profilePoints.reserve( densifiedCurve.constGet()->nCoordinates() );
167 for ( auto it = densifiedCurve.vertices_begin(); it != densifiedCurve.vertices_end(); ++it )
168 {
169 profilePoints.insert( *it );
170 }
171 }
172
173 if ( mFeedback->isCanceled() )
174 return false;
175
176 // calculate the portion of the raster which actually covers the curve
177 int subRegionWidth = 0;
178 int subRegionHeight = 0;
179 int subRegionLeft = 0;
180 int subRegionTop = 0;
181 const QgsRectangle rasterSubRegion = mRasterProvider->xSize() > 0 && mRasterProvider->ySize() > 0 ?
183 mRasterProvider->extent(),
184 mRasterProvider->xSize(),
185 mRasterProvider->ySize(),
186 transformedCurve->boundingBox(),
187 subRegionWidth,
188 subRegionHeight,
189 subRegionLeft,
190 subRegionTop ) : transformedCurve->boundingBox();
191
192 const bool zeroXYSize = mRasterProvider->xSize() == 0 || mRasterProvider->ySize() == 0;
193 if ( zeroXYSize )
194 {
195 const double curveLengthInPixels = sourceCurve->length() / context.mapUnitsPerDistancePixel();
196 const double conversionFactor = curveLengthInPixels / transformedCurve->length();
197 subRegionWidth = rasterSubRegion.width() * conversionFactor;
198 subRegionHeight = rasterSubRegion.height() * conversionFactor;
199 }
200
201 // iterate over the raster blocks, throwing away any which don't intersect the profile curve
202 QgsRasterIterator it( mRasterProvider.get() );
203 // we use smaller tile sizes vs the default, as we will be skipping over tiles which don't intersect the curve at all,
204 // and we expect that to be the VAST majority of the tiles in the raster.
205 // => Smaller tile sizes = more regions we can shortcut over = less pixels to iterate over = faster runtime
206 it.setMaximumTileHeight( 64 );
207 it.setMaximumTileWidth( 64 );
208
209 it.startRasterRead( mBand, subRegionWidth, subRegionHeight, rasterSubRegion );
210
211 const double halfPixelSizeX = mRasterUnitsPerPixelX / 2.0;
212 const double halfPixelSizeY = mRasterUnitsPerPixelY / 2.0;
213 int blockColumns = 0;
214 int blockRows = 0;
215 int blockTopLeftColumn = 0;
216 int blockTopLeftRow = 0;
217 QgsRectangle blockExtent;
218
219 while ( it.next( mBand, blockColumns, blockRows, blockTopLeftColumn, blockTopLeftRow, blockExtent ) )
220 {
221 if ( mFeedback->isCanceled() )
222 return false;
223
224 const QgsGeometry blockExtentGeom = QgsGeometry::fromRect( blockExtent );
225 if ( !curveEngine->intersects( blockExtentGeom.constGet() ) )
226 continue;
227
228 std::unique_ptr< QgsRasterBlock > block( mRasterProvider->block( mBand, blockExtent, blockColumns, blockRows, mFeedback.get() ) );
229 if ( mFeedback->isCanceled() )
230 return false;
231
232 if ( !block )
233 continue;
234
235 bool isNoData = false;
236
237 // there's two potential code paths we use here, depending on if we want to sample at every pixel, or if we only want to
238 // sample at specific points
239 if ( !std::isnan( stepDistance ) )
240 {
241 auto it = profilePoints.begin();
242 while ( it != profilePoints.end() )
243 {
244 if ( mFeedback->isCanceled() )
245 return false;
246
247 // convert point to a pixel and sample, if it's in this block
248 if ( blockExtent.contains( *it ) )
249 {
250 int row, col;
251 if ( zeroXYSize )
252 {
253 row = std::clamp( static_cast< int >( std::round( ( blockExtent.yMaximum() - it->y() ) * blockRows / blockExtent.height() ) ), 0, blockRows - 1 );
254 col = std::clamp( static_cast< int >( std::round( ( it->x() - blockExtent.xMinimum() ) * blockColumns / blockExtent.width() ) ), 0, blockColumns - 1 );
255 }
256 else
257 {
258 row = std::clamp( static_cast< int >( std::round( ( blockExtent.yMaximum() - it->y() ) / mRasterUnitsPerPixelY ) ), 0, blockRows - 1 );
259 col = std::clamp( static_cast< int >( std::round( ( it->x() - blockExtent.xMinimum() ) / mRasterUnitsPerPixelX ) ), 0, blockColumns - 1 );
260 }
261 double val = block->valueAndNoData( row, col, isNoData );
262 if ( !isNoData )
263 {
264 val = val * mScale + mOffset;
265 }
266 else
267 {
268 val = std::numeric_limits<double>::quiet_NaN();
269 }
270
271 QgsPoint pixel( it->x(), it->y(), val );
272 try
273 {
274 pixel.transform( rasterToTargetTransform );
275 }
276 catch ( QgsCsException & )
277 {
278 continue;
279 }
280 mResults->mRawPoints.append( pixel );
281
282 it = profilePoints.erase( it );
283 }
284 else
285 {
286 it++;
287 }
288 }
289
290 if ( profilePoints.isEmpty() )
291 break; // all done!
292 }
293 else
294 {
295 double currentY = blockExtent.yMaximum() - 0.5 * mRasterUnitsPerPixelY;
296 for ( int row = 0; row < blockRows; ++row )
297 {
298 if ( mFeedback->isCanceled() )
299 return false;
300
301 double currentX = blockExtent.xMinimum() + 0.5 * mRasterUnitsPerPixelX;
302 for ( int col = 0; col < blockColumns; ++col, currentX += mRasterUnitsPerPixelX )
303 {
304 const double val = block->valueAndNoData( row, col, isNoData );
305
306 // does pixel intersect curve?
307 QgsGeometry pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - halfPixelSizeX,
308 currentY - halfPixelSizeY,
309 currentX + halfPixelSizeX,
310 currentY + halfPixelSizeY ) );
311 if ( !curveEngine->intersects( pixelRectGeometry.constGet() ) )
312 continue;
313
314 QgsPoint pixel( currentX, currentY, isNoData ? std::numeric_limits<double>::quiet_NaN() : val * mScale + mOffset );
315 try
316 {
317 pixel.transform( rasterToTargetTransform );
318 }
319 catch ( QgsCsException & )
320 {
321 continue;
322 }
323 mResults->mRawPoints.append( pixel );
324 }
325 currentY -= mRasterUnitsPerPixelY;
326 }
327 }
328 }
329
330 if ( mFeedback->isCanceled() )
331 return false;
332
333 // convert x/y values back to distance/height values
334 QgsGeos originalCurveGeos( sourceCurve );
335 originalCurveGeos.prepareGeometry();
336 QString lastError;
337 for ( const QgsPoint &pixel : std::as_const( mResults->mRawPoints ) )
338 {
339 if ( mFeedback->isCanceled() )
340 return false;
341
342 const double distance = originalCurveGeos.lineLocatePoint( pixel, &lastError ) + startDistanceOffset;
343
344 if ( !std::isnan( pixel.z() ) )
345 {
346 mResults->minZ = std::min( pixel.z(), mResults->minZ );
347 mResults->maxZ = std::max( pixel.z(), mResults->maxZ );
348 }
349 mResults->mDistanceToHeightMap.insert( distance, pixel.z() );
350 }
351
352 return true;
353}
354
356{
357 return mResults.release();
358}
359
361{
362 return mFeedback.get();
363}
@ RespectsMaximumErrorMapUnit
Generated profile respects the QgsProfileGenerationContext::maximumErrorMapUnits() property.
@ RespectsDistanceRange
Generated profile respects the QgsProfileGenerationContext::distanceRange() property.
virtual int nCoordinates() const
Returns the number of nodes contained in the geometry.
virtual double length() const
Returns the planar, 2-dimensional length of the geometry.
Abstract base class for storage of elevation profiles.
QVector< QgsProfileIdentifyResults > identify(const QgsProfilePoint &point, const QgsProfileIdentifyContext &context) override
Identify results visible at the specified profile point.
Class for doing transforms between two map coordinate systems.
Custom exception class for Coordinate Reference System related exceptions.
Definition: qgsexception.h:66
Abstract base class for curved geometry type.
Definition: qgscurve.h:36
QgsCurve * clone() const override=0
Clones the geometry by performing a deep copy.
bool isInfinite() const
Returns true if the range consists of all possible values.
Definition: qgsrange.h:247
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition: qgsfeedback.h:45
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.
QgsGeometry densifyByDistance(double distance) const
Densifies the geometry by adding regularly placed extra nodes inside each segment so that the maximum...
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.
static QgsGeometryEngine * createGeometryEngine(const QgsAbstractGeometry *geometry)
Creates and returns a new geometry engine representing the specified geometry.
QgsAbstractGeometry::vertex_iterator vertices_begin() const
Returns STL-style iterator pointing to the first vertex of the geometry.
QgsAbstractGeometry::vertex_iterator vertices_end() const
Returns STL-style iterator pointing to the imaginary vertex after the last vertex of the geometry.
Does vector analysis using the geos library and handles import, export, exception handling*.
Definition: qgsgeos.h:99
void prepareGeometry() override
Prepares the geometry, so that subsequent calls to spatial relation methods are much faster.
Definition: qgsgeos.cpp:252
double lineLocatePoint(const QgsPoint &point, QString *errorMsg=nullptr) const
Returns a distance representing the location along this linestring of the closest point on this lines...
Definition: qgsgeos.cpp:2810
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:49
void transform(const QgsCoordinateTransform &ct, Qgis::TransformDirection d=Qgis::TransformDirection::Forward, bool transformZ=false) override SIP_THROW(QgsCsException)
Transforms the geometry using a coordinate transform.
Definition: qgspoint.cpp:382
Encapsulates the context in which an elevation profile is to be generated.
double maximumErrorMapUnits() const
Returns the maximum allowed error in the generated result, in profile curve map units.
double mapUnitsPerDistancePixel() const
Returns the number of map units per pixel in the distance dimension.
QgsDoubleRange distanceRange() const
Returns the range of distances to include in the generation.
Encapsulates the context of identifying profile results.
Stores identify results generated by a QgsAbstractProfileResults object.
Encapsulates a point on a distance-elevation profile.
Encapsulates properties and constraints relating to fetching elevation profiles from different source...
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
Feedback object tailored for raster block reading.
QgsRasterDataProvider * clone() const override=0
Clone itself, create deep copy.
Iterator for sequentially processing raster cells.
bool next(int bandNumber, int &columns, int &rows, int &topLeftColumn, int &topLeftRow, QgsRectangle &blockExtent)
Fetches details of the next part of the raster data.
void setMaximumTileWidth(int w)
Sets the maximum tile width returned during iteration.
static QgsRectangle subRegion(const QgsRectangle &rasterExtent, int rasterWidth, int rasterHeight, const QgsRectangle &subRegion, int &subRegionWidth, int &subRegionHeight, int &subRegionLeft, int &subRegionTop)
Given an overall raster extent and width and height in pixels, calculates the sub region of the raste...
void startRasterRead(int bandNumber, qgssize nCols, qgssize nRows, const QgsRectangle &extent, QgsRasterBlockFeedback *feedback=nullptr)
Start reading of raster band.
void setMaximumTileHeight(int h)
Sets the minimum tile height returned during iteration.
Raster layer specific subclass of QgsMapLayerElevationProperties.
QgsAbstractProfileResults * takeResults() override
Takes results from the generator.
QgsRasterLayerProfileGenerator(QgsRasterLayer *layer, const QgsProfileRequest &request)
Constructor for QgsRasterLayerProfileGenerator.
QgsFeedback * feedback() const override
Access to feedback object of the generator (may be nullptr)
bool generateProfile(const QgsProfileGenerationContext &context=QgsProfileGenerationContext()) override
Generate the profile (based on data stored in the class).
Qgis::ProfileGeneratorFlags flags() const override
Returns flags which reflect how the profile generator operates.
QString sourceId() const override
Returns a unique identifier representing the source of the profile.
~QgsRasterLayerProfileGenerator() override
QVector< QgsProfileIdentifyResults > identify(const QgsProfilePoint &point, const QgsProfileIdentifyContext &context) override
Identify results visible at the specified profile point.
QString type() const override
Returns the unique string identifier for the results type.
Represents a raster layer.
QgsRasterDataProvider * dataProvider() override
Returns the source data provider.
QgsMapLayerElevationProperties * elevationProperties() override
Returns the layer's elevation properties.
A rectangle specified with double values.
Definition: qgsrectangle.h:42
double yMaximum() const SIP_HOLDGIL
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:193
bool intersects(const QgsRectangle &rect) const SIP_HOLDGIL
Returns true when rectangle intersects with other rectangle.
Definition: qgsrectangle.h:349
double xMinimum() const SIP_HOLDGIL
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:188
double height() const SIP_HOLDGIL
Returns the height of the rectangle.
Definition: qgsrectangle.h:230
double width() const SIP_HOLDGIL
Returns the width of the rectangle.
Definition: qgsrectangle.h:223
bool contains(const QgsRectangle &rect) const SIP_HOLDGIL
Returns true when rectangle contains other rectangle.
Definition: qgsrectangle.h:363
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
const QgsCoordinateReferenceSystem & crs