101 if ( !mProfileCurve || mFeedback->isCanceled() )
109 std::unique_ptr< QgsCurve > trimmedCurve;
111 if ( startDistanceOffset > 0 || endDistance < mProfileCurve->length() )
113 trimmedCurve.reset( mProfileCurve->curveSubstring( startDistanceOffset, endDistance ) );
114 sourceCurve = trimmedCurve.get();
118 sourceCurve = mProfileCurve.get();
122 std::unique_ptr< QgsCurve > transformedCurve( sourceCurve->
clone() );
130 QgsDebugError( QStringLiteral(
"Error transforming profile line to raster CRS" ) );
134 if ( mFeedback->isCanceled() )
137 const QgsRectangle profileCurveBoundingBox = transformedCurve->boundingBox();
138 if ( !profileCurveBoundingBox.
intersects( mRasterProvider->extent() ) )
141 if ( mFeedback->isCanceled() )
144 mResults = std::make_unique< QgsRasterLayerProfileResults >();
145 mResults->mLayer = mLayer;
147 mResults->copyPropertiesFromGenerator(
this );
150 curveEngine->prepareGeometry();
152 if ( mFeedback->isCanceled() )
155 double stepDistance = mStepDistance;
166 QSet< QgsPointXY > profilePoints;
167 if ( !std::isnan( stepDistance ) )
176 profilePoints.insert( *it );
180 if ( mFeedback->isCanceled() )
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(),
197 subRegionTop ) : transformedCurve->boundingBox();
199 const bool zeroXYSize = mRasterProvider->xSize() == 0 || mRasterProvider->ySize() == 0;
203 const double conversionFactor = curveLengthInPixels / transformedCurve->length();
204 subRegionWidth = rasterSubRegion.
width() * conversionFactor;
205 subRegionHeight = rasterSubRegion.
height() * conversionFactor;
209 if ( subRegionWidth == 0 )
214 if ( subRegionHeight == 0 )
229 it.
startRasterRead( mBand, subRegionWidth, subRegionHeight, rasterSubRegion );
231 const double halfPixelSizeX = mRasterUnitsPerPixelX / 2.0;
232 const double halfPixelSizeY = mRasterUnitsPerPixelY / 2.0;
233 int blockColumns = 0;
235 int blockTopLeftColumn = 0;
236 int blockTopLeftRow = 0;
239 while ( it.
next( mBand, blockColumns, blockRows, blockTopLeftColumn, blockTopLeftRow, blockExtent ) )
241 if ( mFeedback->isCanceled() )
245 if ( !curveEngine->intersects( blockExtentGeom.
constGet() ) )
248 std::unique_ptr< QgsRasterBlock > block( mRasterProvider->block( mBand, blockExtent, blockColumns, blockRows, mFeedback.get() ) );
249 if ( mFeedback->isCanceled() )
255 bool isNoData =
false;
259 if ( !std::isnan( stepDistance ) )
261 auto it = profilePoints.begin();
262 while ( it != profilePoints.end() )
264 if ( mFeedback->isCanceled() )
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 );
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 );
281 double val = block->valueAndNoData( row, col, isNoData );
284 val = val * mScale + mOffset;
288 val = std::numeric_limits<double>::quiet_NaN();
291 QgsPoint pixel( it->x(), it->y(), val );
294 pixel.
transform( rasterToTargetTransform );
300 mResults->mRawPoints.append( pixel );
302 it = profilePoints.erase( it );
310 if ( profilePoints.isEmpty() )
315 double currentY = blockExtent.
yMaximum() - 0.5 * mRasterUnitsPerPixelY;
316 for (
int row = 0; row < blockRows; ++row )
318 if ( mFeedback->isCanceled() )
321 double currentX = blockExtent.
xMinimum() + 0.5 * mRasterUnitsPerPixelX;
322 for (
int col = 0; col < blockColumns; ++col, currentX += mRasterUnitsPerPixelX )
324 const double val = block->valueAndNoData( row, col, isNoData );
328 currentY - halfPixelSizeY,
329 currentX + halfPixelSizeX,
330 currentY + halfPixelSizeY ) );
331 if ( !curveEngine->intersects( pixelRectGeometry.
constGet() ) )
334 QgsPoint pixel( currentX, currentY, isNoData ? std::numeric_limits<double>::quiet_NaN() : val * mScale + mOffset );
337 pixel.
transform( rasterToTargetTransform );
343 mResults->mRawPoints.append( pixel );
345 currentY -= mRasterUnitsPerPixelY;
350 if ( mFeedback->isCanceled() )
354 QgsGeos originalCurveGeos( sourceCurve );
357 for (
const QgsPoint &pixel : std::as_const( mResults->mRawPoints ) )
359 if ( mFeedback->isCanceled() )
362 const double distance = originalCurveGeos.
lineLocatePoint( pixel, &lastError ) + startDistanceOffset;
364 if ( !std::isnan( pixel.z() ) )
366 mResults->minZ = std::min( pixel.z(), mResults->minZ );
367 mResults->maxZ = std::max( pixel.z(), mResults->maxZ );
369 mResults->mDistanceToHeightMap.insert( distance, pixel.z() );