Implement a profile generator for QgsRasterLayers

This class implements an optimised method for generating a profile
from a raster layer in a thread-safe way:

On the main thread:

- The data provider is cloned during preparation, and all other
required properties from the raster layer are copied and stored
for later thread-safe use on a background thread

On the background thread:

- The profile line is transformed to the raster's native CRS
- We then iterate over the portion of the raster which intersects
the profile line's bounding box in small tiles. We use small tiles
here as we will shortcut by skipping straight over any tiles which
don't intersect the profile line at all, without requesting their
raster data at all. Since the profile line will only cover a very
small portion of an overall raster extent, by using small tiles
we end up shortcutting and avoiding the costly tile pixel iteration
for most of the raster's coverage.
- For any tiles which DO intersect the profile curve, we fetch the
tile data and then iterate over the pixels, keeping only those
which actually intersect the profile curve. These pixel centroids
are then transformed back to the original CRS of the profile line.
- After collecting the filtered pixels centroids and their raster
(height) values, we then convert the pixel x/y locations to a
distance/chainage along the profile line, giving us an array of
distance vs height values for all pixels which intersect the profile
line.

python/core/auto_generated/raster/qgsrasterlayer.sip.in

@@ -12,7 +12,7 @@
 
 typedef QList < QPair< QString, QColor > > QgsLegendColorList;
 
-class QgsRasterLayer : QgsMapLayer
+class QgsRasterLayer : QgsMapLayer, QgsAbstractProfileSource
 {
 %Docstring(signature="appended")
 
@@ -107,6 +107,9 @@ created for the same data source and renderer is cloned too.
 .. versionadded:: 3.0
 %End
 
+    virtual QgsAbstractProfileGenerator *createProfileGenerator( const QgsProfileRequest &request ) /Factory/;
+
+
     enum ColorShadingAlgorithm
     {
       UndefinedShader,

src/core/CMakeLists.txt

@@ -665,6 +665,7 @@ set(QGIS_CORE_SRCS
   raster/qgsrasteriterator.cpp
   raster/qgsrasterlayer.cpp
   raster/qgsrasterlayerelevationproperties.cpp
+  raster/qgsrasterlayerprofilegenerator.cpp
   raster/qgsrasterlayerrenderer.cpp
   raster/qgsrasterlayertemporalproperties.cpp
   raster/qgsrasterminmaxorigin.cpp
@@ -1676,6 +1677,7 @@ set(QGIS_CORE_HDRS
   raster/qgsrasteriterator.h
   raster/qgsrasterlayer.h
   raster/qgsrasterlayerelevationproperties.h
+  raster/qgsrasterlayerprofilegenerator.h
   raster/qgsrasterlayerrenderer.h
   raster/qgsrasterlayertemporalproperties.h
   raster/qgsrasterminmaxorigin.h

src/core/raster/qgsrasterlayer.cpp

@@ -60,6 +60,7 @@ email                : tim at linfiniti.com
 #include "qgsmaplayerfactory.h"
 #include "qgsrasterpipe.h"
 #include "qgsrasterlayerelevationproperties.h"
+#include "qgsrasterlayerprofilegenerator.h"
 
 #include <cmath>
 #include <cstdio>
@@ -176,6 +177,14 @@ QgsRasterLayer *QgsRasterLayer::clone() const
   return layer;
 }
 
+QgsAbstractProfileGenerator *QgsRasterLayer::createProfileGenerator( const QgsProfileRequest &request )
+{
+  if ( !mElevationProperties->isEnabled() )
+    return nullptr;
+
+  return new QgsRasterLayerProfileGenerator( this, request );
+}
+
 //////////////////////////////////////////////////////////
 //
 // Static Methods and members

src/core/raster/qgsrasterlayer.h

@@ -38,6 +38,7 @@
 #include "qgsrasterviewport.h"
 #include "qgsrasterminmaxorigin.h"
 #include "qgscontrastenhancement.h"
+#include "qgsabstractprofilesource.h"
 
 class QgsMapToPixel;
 class QgsRasterRenderer;
@@ -72,7 +73,7 @@ typedef QList < QPair< QString, QColor > > QgsLegendColorList;
  *     my_raster_layer = QgsRasterLayer("/path/to/file.tif", "my layer")
  * \endcode
  */
-class CORE_EXPORT QgsRasterLayer : public QgsMapLayer
+class CORE_EXPORT QgsRasterLayer : public QgsMapLayer, public QgsAbstractProfileSource
 {
     Q_OBJECT
   public:
@@ -178,6 +179,8 @@ class CORE_EXPORT QgsRasterLayer : public QgsMapLayer
      */
     QgsRasterLayer *clone() const override SIP_FACTORY;
 
+    QgsAbstractProfileGenerator *createProfileGenerator( const QgsProfileRequest &request ) override SIP_FACTORY;
+
     //! \brief This enumerator describes the types of shading that can be used
     enum ColorShadingAlgorithm
     {

src/core/raster/qgsrasterlayerprofilegenerator.cpp

@@ -0,0 +1,161 @@
+/***************************************************************************
+                         qgsrasterlayerprofilegenerator.cpp
+                         ---------------
+    begin                : March 2022
+    copyright            : (C) 2022 by Nyall Dawson
+    email                : nyall dot dawson at gmail dot com
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+#include "qgsrasterlayerprofilegenerator.h"
+#include "qgsprofilerequest.h"
+#include "qgscurve.h"
+#include "qgsrasterlayer.h"
+#include "qgsrasterlayerelevationproperties.h"
+#include "qgsrasteriterator.h"
+#include "qgsgeometryengine.h"
+#include "qgsgeos.h"
+
+QgsRasterLayerProfileGenerator::QgsRasterLayerProfileGenerator( QgsRasterLayer *layer, const QgsProfileRequest &request )
+  : mProfileCurve( request.profileCurve() ? request.profileCurve()->clone() : nullptr )
+  , mSourceCrs( layer->crs() )
+  , mTargetCrs( request.crs() )
+  , mTransformContext( request.transformContext() )
+  , mOffset( layer->elevationProperties()->zOffset() )
+  , mScale( layer->elevationProperties()->zScale() )
+  , mRasterUnitsPerPixelX( layer->rasterUnitsPerPixelX() )
+  , mRasterUnitsPerPixelY( layer->rasterUnitsPerPixelY() )
+{
+  mRasterProvider.reset( layer->dataProvider()->clone() );
+}
+
+QgsRasterLayerProfileGenerator::~QgsRasterLayerProfileGenerator() = default;
+
+bool QgsRasterLayerProfileGenerator::generateProfile()
+{
+  if ( !mProfileCurve )
+    return false;
+
+  // we need to transform the profile curve to the raster's CRS
+  std::unique_ptr< QgsCurve > transformedCurve( mProfileCurve->clone() );
+  const QgsCoordinateTransform rasterToTargetTransform( mSourceCrs, mTargetCrs, mTransformContext );
+  try
+  {
+    transformedCurve->transform( rasterToTargetTransform, Qgis::TransformDirection::Reverse );
+  }
+  catch ( QgsCsException & )
+  {
+    QgsDebugMsg( QStringLiteral( "Error transforming profile line to raster CRS" ) );
+    return false;
+  }
+
+  const QgsRectangle profileCurveBoundingBox = transformedCurve->boundingBox();
+  if ( !profileCurveBoundingBox.intersects( mRasterProvider->extent() ) )
+    return false;
+
+  std::unique_ptr< QgsGeometryEngine > curveEngine( QgsGeometry::createGeometryEngine( transformedCurve.get() ) );
+  curveEngine->prepareGeometry();
+
+  // calculate the portion of the raster which actually covers the curve
+  int subRegionWidth = 0;
+  int subRegionHeight = 0;
+  int subRegionLeft = 0;
+  int subRegionTop = 0;
+  const QgsRectangle rasterSubRegion = QgsRasterIterator::subRegion(
+                                         mRasterProvider->extent(),
+                                         mRasterProvider->xSize(),
+                                         mRasterProvider->ySize(),
+                                         transformedCurve->boundingBox(),
+                                         subRegionWidth,
+                                         subRegionHeight,
+                                         subRegionLeft,
+                                         subRegionTop );
+
+  // iterate over the raster blocks, throwing away any which don't intersect the profile curve
+  QgsRasterIterator it( mRasterProvider.get() );
+  // we use smaller tile sizes vs the default, as we will be skipping over tiles which don't intersect the curve at all,
+  // and we expect that to be the VAST majority of the tiles in the raster.
+  // => Smaller tile sizes = more regions we can shortcut over = less pixels to iterate over = faster runtime
+  it.setMaximumTileHeight( 64 );
+  it.setMaximumTileWidth( 64 );
+
+  it.startRasterRead( mBand, subRegionWidth, subRegionHeight, rasterSubRegion );
+
+  const double halfPixelSizeX = mRasterUnitsPerPixelX / 2.0;
+  const double halfPixelSizeY = mRasterUnitsPerPixelY / 2.0;
+  int blockColumns = 0;
+  int blockRows = 0;
+  int blockTopLeftColumn = 0;
+  int blockTopLeftRow = 0;
+  QgsRectangle blockExtent;
+
+  while ( it.next( mBand, blockColumns, blockRows, blockTopLeftColumn, blockTopLeftRow, blockExtent ) )
+  {
+    const QgsGeometry blockExtentGeom = QgsGeometry::fromRect( blockExtent );
+    if ( !curveEngine->intersects( blockExtentGeom.constGet() ) )
+      continue;
+
+    // TODO -- add feedback argument!
+    std::unique_ptr< QgsRasterBlock > block( mRasterProvider->block( mBand, blockExtent, blockColumns, blockRows ) );
+    if ( !block )
+      continue;
+
+    bool isNoData = false;
+
+    double currentY = blockExtent.yMaximum() - 0.5 * mRasterUnitsPerPixelY;
+    for ( int row = 0; row < blockRows; ++row )
+    {
+      double currentX = blockExtent.xMinimum() + 0.5 * mRasterUnitsPerPixelX;
+      for ( int col = 0; col < blockColumns; ++col, currentX += mRasterUnitsPerPixelX )
+      {
+        const double val = block->valueAndNoData( row, col, isNoData );
+        if ( isNoData )
+          continue;
+
+        // does pixel intersect curve?
+        QgsGeometry pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - halfPixelSizeX,
+                                        currentY - halfPixelSizeY,
+                                        currentX + halfPixelSizeX,
+                                        currentY + halfPixelSizeY ) );
+        if ( !curveEngine->intersects( pixelRectGeometry.constGet() ) )
+          continue;
+
+        QgsPoint pixel( currentX, currentY, val * mScale + mOffset );
+        try
+        {
+          pixel.transform( rasterToTargetTransform );
+        }
+        catch ( QgsCsException & )
+        {
+          continue;
+        }
+        mRawPoints.append( pixel );
+      }
+      currentY -= mRasterUnitsPerPixelY;
+    }
+  }
+
+  // convert x/y values back to distance/height values
+  QgsGeos originalCurveGeos( mProfileCurve.get() );
+  originalCurveGeos.prepareGeometry();
+  mResults.reserve( mRawPoints.size() );
+  QString lastError;
+  for ( const QgsPoint &pixel : std::as_const( mRawPoints ) )
+  {
+    const double distance = originalCurveGeos.lineLocatePoint( pixel, &lastError );
+
+    Result res;
+    res.distance = distance;
+    res.height = pixel.z();
+    mResults.push_back( res );
+  }
+
+  return true;
+}

src/core/raster/qgsrasterlayerprofilegenerator.h

@@ -0,0 +1,77 @@
+/***************************************************************************
+                         qgsrasterlayerprofilegenerator.h
+                         ---------------
+    begin                : March 2022
+    copyright            : (C) 2022 by Nyall Dawson
+    email                : nyall dot dawson at gmail dot com
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+#ifndef QGSRASTERLAYERPROFILEGENERATOR_H
+#define QGSRASTERLAYERPROFILEGENERATOR_H
+
+#include "qgis_core.h"
+#include "qgis_sip.h"
+#include "qgsabstractprofilegenerator.h"
+#include "qgscoordinatereferencesystem.h"
+#include "qgscoordinatetransformcontext.h"
+
+#include <memory>
+
+class QgsProfileRequest;
+class QgsCurve;
+class QgsRasterLayer;
+class QgsRasterDataProvider;
+
+
+#define SIP_NO_FILE
+
+/**
+ * \brief Implementation of QgsAbstractProfileGenerator for raster layers.
+ *
+ * \note Not available in Python bindings
+ * \ingroup core
+ * \since QGIS 3.26
+ */
+class CORE_EXPORT QgsRasterLayerProfileGenerator : public QgsAbstractProfileGenerator
+{
+
+  public:
+
+    /**
+     * Constructor for QgsRasterLayerProfileGenerator.
+     */
+    QgsRasterLayerProfileGenerator( QgsRasterLayer *layer, const QgsProfileRequest &request );
+
+    ~QgsRasterLayerProfileGenerator() override;
+
+    bool generateProfile() override;
+
+  private:
+
+    std::unique_ptr< QgsCurve > mProfileCurve;
+
+    QgsCoordinateReferenceSystem mSourceCrs;
+    QgsCoordinateReferenceSystem mTargetCrs;
+    QgsCoordinateTransformContext mTransformContext;
+
+    double mOffset = 0;
+    double mScale = 1;
+
+    std::unique_ptr< QgsRasterDataProvider > mRasterProvider;
+
+    int mBand = 1;
+    double mRasterUnitsPerPixelX = 1;
+    double mRasterUnitsPerPixelY = 1;
+
+
+};
+
+#endif // QGSRASTERLAYERPROFILEGENERATOR_H

tests/src/python/CMakeLists.txt

@@ -277,6 +277,7 @@ ADD_PYTHON_TEST(PyQgsRasterFileWriter test_qgsrasterfilewriter.py)
 ADD_PYTHON_TEST(PyQgsRasterFileWriterTask test_qgsrasterfilewritertask.py)
 ADD_PYTHON_TEST(PyQgsRasterLayer test_qgsrasterlayer.py)
 ADD_PYTHON_TEST(PyQgsRasterLayerElevationProperties test_qgsrasterlayerelevationproperties.py)
+ADD_PYTHON_TEST(PyQgsRasterLayerProfileGenerator test_qgsrasterlayerprofilegenerator.py)
 ADD_PYTHON_TEST(PyQgsRasterLayerProperties test_qgsrasterlayerproperties.py)
 ADD_PYTHON_TEST(PyQgsRasterLayerRenderer test_qgsrasterlayerrenderer.py)
 ADD_PYTHON_TEST(PyQgsRasterColorRampShader test_qgsrastercolorrampshader.py)

tests/src/python/test_qgsrasterlayerprofilegenerator.py

@@ -0,0 +1,71 @@
+# -*- coding: utf-8 -*-
+"""QGIS Unit tests for QgsRasterLayer profile generation
+
+.. note:: This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+"""
+__author__ = 'Nyall Dawson'
+__date__ = '18/03/2022'
+__copyright__ = 'Copyright 2022, The QGIS Project'
+
+import os
+
+import qgis  # NOQA
+
+from qgis.PyQt.QtCore import QTemporaryDir
+
+from qgis.core import (
+    QgsRasterLayer,
+    QgsLineString,
+    QgsProfileRequest,
+    QgsCoordinateReferenceSystem,
+    QgsCoordinateTransformContext,
+    QgsFlatTerrainProvider,
+    QgsMeshTerrainProvider
+)
+
+from qgis.PyQt.QtXml import QDomDocument
+
+from qgis.testing import start_app, unittest
+from utilities import unitTestDataPath
+
+start_app()
+
+
+class TestQgsRasterLayerProfileGenerator(unittest.TestCase):
+
+    def testGeneration(self):
+        rl = QgsRasterLayer(os.path.join(unitTestDataPath(), '3d', 'dtm.tif'), 'DTM')
+        self.assertTrue(rl.isValid())
+
+        curve = QgsLineString()
+        curve.fromWkt('LineString (-348095.18706532847136259 6633687.0235139261931181, -347271.57799367723055184 6633093.13086318597197533, -346140.60267287614988163 6632697.89590711053460836, -345777.013075890194159 6631575.50219972990453243)')
+        req = QgsProfileRequest(curve)
+
+        generator = rl.createProfileGenerator(req)
+        # not enabled for elevation
+        self.assertIsNone(generator)
+
+        rl.elevationProperties().setEnabled(True)
+
+        generator = rl.createProfileGenerator(req)
+        self.assertIsNotNone(generator)
+        # the request did not have the crs of the linestring set, so the whole linestring falls outside the raster extent
+        self.assertFalse(generator.generateProfile())
+
+        # set correct crs for linestring and re-try
+        req.setCrs(QgsCoordinateReferenceSystem('EPSG:3857'))
+        generator = rl.createProfileGenerator(req)
+        self.assertTrue(generator.generateProfile())
+
+        results = {r.distance: r.height for r in generator.results()}
+        first_point = min(results.keys())
+        last_point = max(results.keys())
+        self.assertEqual(results[first_point], 154)
+        self.assertEqual(results[last_point], 99)
+
+
+if __name__ == '__main__':
+    unittest.main()