QGIS API Documentation 3.40.0-Bratislava (b56115d8743)
Loading...
Searching...
No Matches
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 "qgsprofilepoint.h"
27#include "qgsfillsymbol.h"
28#include "qgsthreadingutils.h"
29
30#include <QPolygonF>
31#include <QThread>
32
33//
34// QgsRasterLayerProfileResults
35//
36
38{
39 return QStringLiteral( "raster" );
40}
41
42QVector<QgsProfileIdentifyResults> QgsRasterLayerProfileResults::identify( const QgsProfilePoint &point, const QgsProfileIdentifyContext &context )
43{
44 const QVector<QgsProfileIdentifyResults> noLayerResults = QgsAbstractProfileSurfaceResults::identify( point, context );
45
46 // we have to make a new list, with the correct layer reference set
47 QVector<QgsProfileIdentifyResults> res;
48 res.reserve( noLayerResults.size() );
49 for ( const QgsProfileIdentifyResults &result : noLayerResults )
50 {
51 res.append( QgsProfileIdentifyResults( mLayer, result.results() ) );
52 }
53 return res;
54}
55
56
57
58//
59// QgsRasterLayerProfileGenerator
60//
61
64 , mId( layer->id() )
65 , mFeedback( std::make_unique< QgsRasterBlockFeedback >() )
66 , mProfileCurve( request.profileCurve() ? request.profileCurve()->clone() : nullptr )
67 , mSourceCrs( layer->crs() )
68 , mTargetCrs( request.crs() )
69 , mTransformContext( request.transformContext() )
70 , mOffset( layer->elevationProperties()->zOffset() )
71 , mScale( layer->elevationProperties()->zScale() )
72 , mLayer( layer )
73 , mBand( qgis::down_cast< QgsRasterLayerElevationProperties * >( layer->elevationProperties() )->bandNumber() )
74 , mRasterUnitsPerPixelX( layer->rasterUnitsPerPixelX() )
75 , mRasterUnitsPerPixelY( layer->rasterUnitsPerPixelY() )
76 , mStepDistance( request.stepDistance() )
77{
78 mRasterProvider.reset( layer->dataProvider()->clone() );
79 mRasterProvider->moveToThread( nullptr );
80
81 mSymbology = qgis::down_cast< QgsRasterLayerElevationProperties * >( layer->elevationProperties() )->profileSymbology();
82 mElevationLimit = qgis::down_cast< QgsRasterLayerElevationProperties * >( layer->elevationProperties() )->elevationLimit();
83 mLineSymbol.reset( qgis::down_cast< QgsRasterLayerElevationProperties * >( layer->elevationProperties() )->profileLineSymbol()->clone() );
84 mFillSymbol.reset( qgis::down_cast< QgsRasterLayerElevationProperties * >( layer->elevationProperties() )->profileFillSymbol()->clone() );
85}
86
88{
89 return mId;
90}
91
96
98
100{
101 if ( !mProfileCurve || mFeedback->isCanceled() )
102 return false;
103
104 QgsScopedAssignObjectToCurrentThread assignProviderToCurrentThread( mRasterProvider.get() );
105
106 const double startDistanceOffset = std::max( !context.distanceRange().isInfinite() ? context.distanceRange().lower() : 0, 0.0 );
107 const double endDistance = context.distanceRange().upper();
108
109 std::unique_ptr< QgsCurve > trimmedCurve;
110 QgsCurve *sourceCurve = nullptr;
111 if ( startDistanceOffset > 0 || endDistance < mProfileCurve->length() )
112 {
113 trimmedCurve.reset( mProfileCurve->curveSubstring( startDistanceOffset, endDistance ) );
114 sourceCurve = trimmedCurve.get();
115 }
116 else
117 {
118 sourceCurve = mProfileCurve.get();
119 }
120
121 // we need to transform the profile curve to the raster's CRS
122 std::unique_ptr< QgsCurve > transformedCurve( sourceCurve->clone() );
123 const QgsCoordinateTransform rasterToTargetTransform( mSourceCrs, mTargetCrs, mTransformContext );
124 try
125 {
126 transformedCurve->transform( rasterToTargetTransform, Qgis::TransformDirection::Reverse );
127 }
128 catch ( QgsCsException & )
129 {
130 QgsDebugError( QStringLiteral( "Error transforming profile line to raster CRS" ) );
131 return false;
132 }
133
134 if ( mFeedback->isCanceled() )
135 return false;
136
137 const QgsRectangle profileCurveBoundingBox = transformedCurve->boundingBox();
138 if ( !profileCurveBoundingBox.intersects( mRasterProvider->extent() ) )
139 return false;
140
141 if ( mFeedback->isCanceled() )
142 return false;
143
144 mResults = std::make_unique< QgsRasterLayerProfileResults >();
145 mResults->mLayer = mLayer;
146 mResults->mId = mId;
147 mResults->copyPropertiesFromGenerator( this );
148
149 std::unique_ptr< QgsGeometryEngine > curveEngine( QgsGeometry::createGeometryEngine( transformedCurve.get() ) );
150 curveEngine->prepareGeometry();
151
152 if ( mFeedback->isCanceled() )
153 return false;
154
155 double stepDistance = mStepDistance;
156 if ( !std::isnan( context.maximumErrorMapUnits() ) )
157 {
158 // convert the maximum error in curve units to a step distance
159 // TODO -- there's no point in this being << pixel size!
160 if ( std::isnan( stepDistance ) || context.maximumErrorMapUnits() > stepDistance )
161 {
162 stepDistance = context.maximumErrorMapUnits();
163 }
164 }
165
166 QSet< QgsPointXY > profilePoints;
167 if ( !std::isnan( stepDistance ) )
168 {
169 // if specific step distance specified, use this to generate points along the curve
170 QgsGeometry densifiedCurve( sourceCurve->clone() );
171 densifiedCurve = densifiedCurve.densifyByDistance( stepDistance );
172 densifiedCurve.transform( rasterToTargetTransform, Qgis::TransformDirection::Reverse );
173 profilePoints.reserve( densifiedCurve.constGet()->nCoordinates() );
174 for ( auto it = densifiedCurve.vertices_begin(); it != densifiedCurve.vertices_end(); ++it )
175 {
176 profilePoints.insert( *it );
177 }
178 }
179
180 if ( mFeedback->isCanceled() )
181 return false;
182
183 // calculate the portion of the raster which actually covers the curve
184 int subRegionWidth = 0;
185 int subRegionHeight = 0;
186 int subRegionLeft = 0;
187 int subRegionTop = 0;
188 QgsRectangle rasterSubRegion = mRasterProvider->xSize() > 0 && mRasterProvider->ySize() > 0 ?
190 mRasterProvider->extent(),
191 mRasterProvider->xSize(),
192 mRasterProvider->ySize(),
193 transformedCurve->boundingBox(),
194 subRegionWidth,
195 subRegionHeight,
196 subRegionLeft,
197 subRegionTop ) : transformedCurve->boundingBox();
198
199 const bool zeroXYSize = mRasterProvider->xSize() == 0 || mRasterProvider->ySize() == 0;
200 if ( zeroXYSize )
201 {
202 const double curveLengthInPixels = sourceCurve->length() / context.mapUnitsPerDistancePixel();
203 const double conversionFactor = curveLengthInPixels / transformedCurve->length();
204 subRegionWidth = rasterSubRegion.width() * conversionFactor;
205 subRegionHeight = rasterSubRegion.height() * conversionFactor;
206
207 // ensure we fetch at least 1 pixel wide/high blocks, otherwise exactly vertical/horizontal profile lines will result in zero size blocks
208 // see https://github.com/qgis/QGIS/issues/51196
209 if ( subRegionWidth == 0 )
210 {
211 subRegionWidth = 1;
212 rasterSubRegion.setXMaximum( rasterSubRegion.xMinimum() + 1 / conversionFactor );
213 }
214 if ( subRegionHeight == 0 )
215 {
216 subRegionHeight = 1;
217 rasterSubRegion.setYMaximum( rasterSubRegion.yMinimum() + 1 / conversionFactor );
218 }
219 }
220
221 // iterate over the raster blocks, throwing away any which don't intersect the profile curve
222 QgsRasterIterator it( mRasterProvider.get() );
223 // we use smaller tile sizes vs the default, as we will be skipping over tiles which don't intersect the curve at all,
224 // and we expect that to be the VAST majority of the tiles in the raster.
225 // => Smaller tile sizes = more regions we can shortcut over = less pixels to iterate over = faster runtime
226 it.setMaximumTileHeight( 64 );
227 it.setMaximumTileWidth( 64 );
228
229 it.startRasterRead( mBand, subRegionWidth, subRegionHeight, rasterSubRegion );
230
231 const double halfPixelSizeX = mRasterUnitsPerPixelX / 2.0;
232 const double halfPixelSizeY = mRasterUnitsPerPixelY / 2.0;
233 int blockColumns = 0;
234 int blockRows = 0;
235 int blockTopLeftColumn = 0;
236 int blockTopLeftRow = 0;
237 QgsRectangle blockExtent;
238
239 while ( it.next( mBand, blockColumns, blockRows, blockTopLeftColumn, blockTopLeftRow, blockExtent ) )
240 {
241 if ( mFeedback->isCanceled() )
242 return false;
243
244 const QgsGeometry blockExtentGeom = QgsGeometry::fromRect( blockExtent );
245 if ( !curveEngine->intersects( blockExtentGeom.constGet() ) )
246 continue;
247
248 std::unique_ptr< QgsRasterBlock > block( mRasterProvider->block( mBand, blockExtent, blockColumns, blockRows, mFeedback.get() ) );
249 if ( mFeedback->isCanceled() )
250 return false;
251
252 if ( !block )
253 continue;
254
255 bool isNoData = false;
256
257 // 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
258 // sample at specific points
259 if ( !std::isnan( stepDistance ) )
260 {
261 auto it = profilePoints.begin();
262 while ( it != profilePoints.end() )
263 {
264 if ( mFeedback->isCanceled() )
265 return false;
266
267 // convert point to a pixel and sample, if it's in this block
268 if ( blockExtent.contains( *it ) )
269 {
270 int row, col;
271 if ( zeroXYSize )
272 {
273 row = std::clamp( static_cast< int >( std::round( ( blockExtent.yMaximum() - it->y() ) * blockRows / blockExtent.height() ) ), 0, blockRows - 1 );
274 col = std::clamp( static_cast< int >( std::round( ( it->x() - blockExtent.xMinimum() ) * blockColumns / blockExtent.width() ) ), 0, blockColumns - 1 );
275 }
276 else
277 {
278 row = std::clamp( static_cast< int >( std::round( ( blockExtent.yMaximum() - it->y() ) / mRasterUnitsPerPixelY ) ), 0, blockRows - 1 );
279 col = std::clamp( static_cast< int >( std::round( ( it->x() - blockExtent.xMinimum() ) / mRasterUnitsPerPixelX ) ), 0, blockColumns - 1 );
280 }
281 double val = block->valueAndNoData( row, col, isNoData );
282 if ( !isNoData )
283 {
284 val = val * mScale + mOffset;
285 }
286 else
287 {
288 val = std::numeric_limits<double>::quiet_NaN();
289 }
290
291 QgsPoint pixel( it->x(), it->y(), val );
292 try
293 {
294 pixel.transform( rasterToTargetTransform );
295 }
296 catch ( QgsCsException & )
297 {
298 continue;
299 }
300 mResults->mRawPoints.append( pixel );
301
302 it = profilePoints.erase( it );
303 }
304 else
305 {
306 it++;
307 }
308 }
309
310 if ( profilePoints.isEmpty() )
311 break; // all done!
312 }
313 else
314 {
315 double currentY = blockExtent.yMaximum() - 0.5 * mRasterUnitsPerPixelY;
316 for ( int row = 0; row < blockRows; ++row )
317 {
318 if ( mFeedback->isCanceled() )
319 return false;
320
321 double currentX = blockExtent.xMinimum() + 0.5 * mRasterUnitsPerPixelX;
322 for ( int col = 0; col < blockColumns; ++col, currentX += mRasterUnitsPerPixelX )
323 {
324 const double val = block->valueAndNoData( row, col, isNoData );
325
326 // does pixel intersect curve?
327 QgsGeometry pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - halfPixelSizeX,
328 currentY - halfPixelSizeY,
329 currentX + halfPixelSizeX,
330 currentY + halfPixelSizeY ) );
331 if ( !curveEngine->intersects( pixelRectGeometry.constGet() ) )
332 continue;
333
334 QgsPoint pixel( currentX, currentY, isNoData ? std::numeric_limits<double>::quiet_NaN() : val * mScale + mOffset );
335 try
336 {
337 pixel.transform( rasterToTargetTransform );
338 }
339 catch ( QgsCsException & )
340 {
341 continue;
342 }
343 mResults->mRawPoints.append( pixel );
344 }
345 currentY -= mRasterUnitsPerPixelY;
346 }
347 }
348 }
349
350 if ( mFeedback->isCanceled() )
351 return false;
352
353 // convert x/y values back to distance/height values
354 QgsGeos originalCurveGeos( sourceCurve );
355 originalCurveGeos.prepareGeometry();
356 QString lastError;
357 for ( const QgsPoint &pixel : std::as_const( mResults->mRawPoints ) )
358 {
359 if ( mFeedback->isCanceled() )
360 return false;
361
362 const double distance = originalCurveGeos.lineLocatePoint( pixel, &lastError ) + startDistanceOffset;
363
364 if ( !std::isnan( pixel.z() ) )
365 {
366 mResults->minZ = std::min( pixel.z(), mResults->minZ );
367 mResults->maxZ = std::max( pixel.z(), mResults->maxZ );
368 }
369 mResults->mDistanceToHeightMap.insert( distance, pixel.z() );
370 }
371
372 return true;
373}
374
379
381{
382 return mFeedback.get();
383}
@ RespectsMaximumErrorMapUnit
Generated profile respects the QgsProfileGenerationContext::maximumErrorMapUnits() property.
@ RespectsDistanceRange
Generated profile respects the QgsProfileGenerationContext::distanceRange() property.
QFlags< ProfileGeneratorFlag > ProfileGeneratorFlags
Definition qgis.h:3920
@ Reverse
Reverse/inverse transform (from destination to source)
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.
Abstract base class for objects which generate elevation profiles which represent a continuous surfac...
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.
Abstract base class for curved geometry type.
Definition qgscurve.h:35
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:285
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition qgsfeedback.h:44
A geometry is the spatial representation of a feature.
static QgsGeometry fromRect(const QgsRectangle &rect)
Creates a new geometry from a QgsRectangle.
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)
Transforms this geometry as described by the coordinate transform ct.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
QgsAbstractGeometry::vertex_iterator vertices_begin() const
Returns STL-style iterator pointing to the first vertex of the geometry.
static QgsGeometryEngine * createGeometryEngine(const QgsAbstractGeometry *geometry, double precision=0.0)
Creates and returns a new geometry engine representing the specified geometry using precision on a gr...
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:137
void prepareGeometry() override
Prepares the geometry, so that subsequent calls to spatial relation methods are much faster.
Definition qgsgeos.cpp:289
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:3153
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
Transforms the geometry using a coordinate transform.
Definition qgspoint.cpp:383
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:78
T upper() const
Returns the upper bound of the range.
Definition qgsrange.h:85
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.
bool contains(const QgsRectangle &rect) const
Returns true when rectangle contains other rectangle.
double xMinimum() const
Returns the x minimum value (left side of rectangle).
bool intersects(const QgsRectangle &rect) const
Returns true when rectangle intersects with other rectangle.
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
double width() const
Returns the width of the rectangle.
double yMaximum() const
Returns the y maximum value (top side of rectangle).
void setYMaximum(double y)
Set the maximum y value.
void setXMaximum(double x)
Set the maximum x value.
double height() const
Returns the height of the rectangle.
Temporarily moves a QObject to the current thread, then resets it back to nullptr thread on destructi...
#define QgsDebugError(str)
Definition qgslogger.h:38
const QgsCoordinateReferenceSystem & crs