diff --git a/python/core/auto_generated/raster/qgsrasterlayer.sip.in b/python/core/auto_generated/raster/qgsrasterlayer.sip.in index 3bf9b65c26f..87cac5176a6 100644 --- a/python/core/auto_generated/raster/qgsrasterlayer.sip.in +++ b/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, diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 459eba5361d..b7e7811b3c0 100644 --- a/src/core/CMakeLists.txt +++ b/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 diff --git a/src/core/raster/qgsrasterlayer.cpp b/src/core/raster/qgsrasterlayer.cpp index 80496cb4270..00fe247c745 100644 --- a/src/core/raster/qgsrasterlayer.cpp +++ b/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 #include @@ -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 diff --git a/src/core/raster/qgsrasterlayer.h b/src/core/raster/qgsrasterlayer.h index 9480f17e96d..eae74939c9c 100644 --- a/src/core/raster/qgsrasterlayer.h +++ b/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 { diff --git a/src/core/raster/qgsrasterlayerprofilegenerator.cpp b/src/core/raster/qgsrasterlayerprofilegenerator.cpp new file mode 100644 index 00000000000..83e6f104a94 --- /dev/null +++ b/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; +} diff --git a/src/core/raster/qgsrasterlayerprofilegenerator.h b/src/core/raster/qgsrasterlayerprofilegenerator.h new file mode 100644 index 00000000000..572d6c124e7 --- /dev/null +++ b/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 + +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 diff --git a/tests/src/python/CMakeLists.txt b/tests/src/python/CMakeLists.txt index a42b830be18..a5e14382ca0 100644 --- a/tests/src/python/CMakeLists.txt +++ b/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) diff --git a/tests/src/python/test_qgsrasterlayerprofilegenerator.py b/tests/src/python/test_qgsrasterlayerprofilegenerator.py new file mode 100644 index 00000000000..c25dac61dc5 --- /dev/null +++ b/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()