QGIS API Documentation 3.38.0-Grenoble (exported)
Loading...
Searching...
No Matches
qgsalgorithmrasterdtmslopebasedfilter.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsalgorithmrasterdtmslopebasedfilter.cpp
3 ---------------------
4 begin : July 2023
5 copyright : (C) 2023 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 ***************************************************************************/
17
19#include "qgsrasterfilewriter.h"
20#include <algorithm>
21
23
24QString QgsRasterDtmSlopeBasedFilterAlgorithm::name() const
25{
26 return QStringLiteral( "dtmslopebasedfilter" );
27}
28
29QString QgsRasterDtmSlopeBasedFilterAlgorithm::displayName() const
30{
31 return QObject::tr( "DTM filter (slope-based)" );
32}
33
34QStringList QgsRasterDtmSlopeBasedFilterAlgorithm::tags() const
35{
36 return QObject::tr( "dem,filter,slope,dsm,dtm,terrain" ).split( ',' );
37}
38
39QString QgsRasterDtmSlopeBasedFilterAlgorithm::group() const
40{
41 return QObject::tr( "Raster terrain analysis" );
42}
43
44QString QgsRasterDtmSlopeBasedFilterAlgorithm::groupId() const
45{
46 return QStringLiteral( "rasterterrainanalysis" );
47}
48
49QString QgsRasterDtmSlopeBasedFilterAlgorithm::shortHelpString() const
50{
51 return QObject::tr( "This algorithm can be used to filter a digital elevation model in order to classify its cells into ground and object (non-ground) cells.\n\n"
52 "The tool uses concepts as described by Vosselman (2000) and is based on the assumption that a large height difference between two nearby "
53 "cells is unlikely to be caused by a steep slope in the terrain. The probability that the higher cell might be non-ground increases when "
54 "the distance between the two cells decreases. Therefore the filter defines a maximum height difference (<i>dz_max</i>) between two cells as a "
55 "function of the distance (<i>d</i>) between the cells (<i>dz_max( d ) = d</i>).\n\n"
56 "A cell is classified as terrain if there is no cell within the kernel radius to which the height difference is larger than the allowed "
57 "maximum height difference at the distance between these two cells.\n\n"
58 "The approximate terrain slope (<i>s</i>) parameter is used to modify the filter function to match the overall slope in the study "
59 "area (<i>dz_max( d ) = d * s</i>).\n\n"
60 "A 5 % confidence interval (<i>ci = 1.65 * sqrt( 2 * stddev )</i>) may be used to modify the filter function even further by either "
61 "relaxing (<i>dz_max( d ) = d * s + ci</i>) or amplifying (<i>dz_max( d ) = d * s - ci</i>) the filter criterium.\n\n"
62 "References: Vosselman, G. (2000): Slope based filtering of laser altimetry data. IAPRS, Vol. XXXIII, Part B3, Amsterdam, The Netherlands, 935-942\n\n"
63 "This algorithm is a port of the SAGA 'DTM Filter (slope-based)' tool." );
64}
65
66void QgsRasterDtmSlopeBasedFilterAlgorithm::initAlgorithm( const QVariantMap & )
67{
68 addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT" ),
69 QObject::tr( "Input layer" ) ) );
70
71 addParameter( new QgsProcessingParameterBand( QStringLiteral( "BAND" ),
72 QObject::tr( "Band number" ), 1, QStringLiteral( "INPUT" ) ) );
73
74 std::unique_ptr< QgsProcessingParameterNumber > radiusParam = std::make_unique< QgsProcessingParameterNumber >( QStringLiteral( "RADIUS" ),
75 QObject::tr( "Kernel radius (pixels)" ), Qgis::ProcessingNumberParameterType::Integer, 5, false, 1, 1000 );
76 radiusParam->setHelp( QObject::tr( "The radius of the filter kernel (in pixels). Must be large enough to reach ground cells next to non-ground objects." ) );
77 addParameter( radiusParam.release() );
78
79 std::unique_ptr< QgsProcessingParameterNumber > terrainSlopeParam = std::make_unique< QgsProcessingParameterNumber >( QStringLiteral( "TERRAIN_SLOPE" ),
80 QObject::tr( "Terrain slope (%, pixel size/vertical units)" ), Qgis::ProcessingNumberParameterType::Double, 30, false, 0, 1000 );
81 terrainSlopeParam->setHelp( QObject::tr( "The approximate terrain slope in %. The terrain slope must be adjusted to account for the ratio of height units vs raster pixel dimensions. Used to relax the filter criterium in steeper terrain." ) );
82 addParameter( terrainSlopeParam.release() );
83
84 std::unique_ptr< QgsProcessingParameterEnum > filterModificationParam = std::make_unique< QgsProcessingParameterEnum >( QStringLiteral( "FILTER_MODIFICATION" ),
85 QObject::tr( "Filter modification" ), QStringList{ QObject::tr( "None" ), QObject::tr( "Relax filter" ), QObject::tr( "Amplify" ) }, false, 0 );
86 filterModificationParam->setHelp( QObject::tr( "Choose whether to apply the filter kernel without modification or to use a confidence interval to relax or amplify the height criterium." ) );
87 addParameter( filterModificationParam.release() );
88
89 std::unique_ptr< QgsProcessingParameterNumber > stDevParam = std::make_unique< QgsProcessingParameterNumber >( QStringLiteral( "STANDARD_DEVIATION" ),
90 QObject::tr( "Standard deviation" ), Qgis::ProcessingNumberParameterType::Double, 0.1, false, 0, 1000 );
91 stDevParam->setHelp( QObject::tr( "The standard deviation used to calculate a 5% confidence interval applied to the height threshold." ) );
92 addParameter( stDevParam.release() );
93
94 std::unique_ptr< QgsProcessingParameterRasterDestination > outputLayerGroundParam = std::make_unique<QgsProcessingParameterRasterDestination>( QStringLiteral( "OUTPUT_GROUND" ),
95 QObject::tr( "Output layer (ground)" ), QVariant(), true, true );
96 outputLayerGroundParam->setHelp( QObject::tr( "The filtered DEM containing only cells classified as ground." ) );
97 addParameter( outputLayerGroundParam.release() );
98
99 std::unique_ptr< QgsProcessingParameterRasterDestination > outputLayerNonGroundParam = std::make_unique<QgsProcessingParameterRasterDestination>( QStringLiteral( "OUTPUT_NONGROUND" ),
100 QObject::tr( "Output layer (non-ground objects)" ), QVariant(), true, false );
101 outputLayerNonGroundParam->setHelp( QObject::tr( "The non-ground objects removed by the filter." ) );
102 addParameter( outputLayerNonGroundParam.release() );
103}
104
105QgsRasterDtmSlopeBasedFilterAlgorithm *QgsRasterDtmSlopeBasedFilterAlgorithm::createInstance() const
106{
107 return new QgsRasterDtmSlopeBasedFilterAlgorithm();
108}
109
110bool QgsRasterDtmSlopeBasedFilterAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
111{
112 QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT" ), context );
113 if ( !layer )
114 throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT" ) ) );
115
116 const int band = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
117
118 mBand = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
119 if ( mBand < 1 || mBand > layer->bandCount() )
120 throw QgsProcessingException( QObject::tr( "Invalid band number for BAND (%1): Valid values for input raster are 1 to %2" ).arg( mBand )
121 .arg( layer->bandCount() ) );
122
123 mInterface.reset( layer->dataProvider()->clone() );
124 mHasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( band );
125 mLayerWidth = layer->width();
126 mLayerHeight = layer->height();
127 mExtent = layer->extent();
128 mCrs = layer->crs();
129 mRasterUnitsPerPixelX = layer->rasterUnitsPerPixelX();
130 mRasterUnitsPerPixelY = layer->rasterUnitsPerPixelY();
131 mDataType = layer->dataProvider()->dataType( mBand );
132 mNoData = layer->dataProvider()->sourceNoDataValue( mBand );
133 return true;
134}
135
136QVariantMap QgsRasterDtmSlopeBasedFilterAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
137{
138 const QString groundOutputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT_GROUND" ), context );
139 std::unique_ptr< QgsRasterFileWriter > groundWriter;
140 std::unique_ptr<QgsRasterDataProvider > groundDestProvider;
141
142 if ( !groundOutputFile.isEmpty() )
143 {
144 const QFileInfo fi( groundOutputFile );
145 const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
146
147 groundWriter = std::make_unique< QgsRasterFileWriter >( groundOutputFile );
148 groundWriter->setOutputProviderKey( QStringLiteral( "gdal" ) );
149 groundWriter->setOutputFormat( outputFormat );
150
151 groundDestProvider.reset( groundWriter->createOneBandRaster( mDataType, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
152
153 if ( !groundDestProvider )
154 throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( groundOutputFile ) );
155 if ( !groundDestProvider->isValid() )
156 throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( groundOutputFile, groundDestProvider->error().message( QgsErrorMessage::Text ) ) );
157
158 groundDestProvider->setNoDataValue( 1, mNoData );
159 groundDestProvider->setEditable( true );
160 }
161
162 const QString nonGroundOutputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT_NONGROUND" ), context );
163 std::unique_ptr< QgsRasterFileWriter > nonGroundWriter;
164 std::unique_ptr<QgsRasterDataProvider > nonGroundDestProvider;
165
166 if ( !nonGroundOutputFile.isEmpty() )
167 {
168 const QFileInfo fi( groundOutputFile );
169 const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
170
171 nonGroundWriter = std::make_unique< QgsRasterFileWriter >( nonGroundOutputFile );
172 nonGroundWriter->setOutputProviderKey( QStringLiteral( "gdal" ) );
173 nonGroundWriter->setOutputFormat( outputFormat );
174
175 nonGroundDestProvider.reset( nonGroundWriter->createOneBandRaster( mDataType, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
176
177 if ( !nonGroundDestProvider )
178 throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( nonGroundOutputFile ) );
179 if ( !nonGroundDestProvider->isValid() )
180 throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( nonGroundOutputFile, nonGroundDestProvider->error().message( QgsErrorMessage::Text ) ) );
181
182 nonGroundDestProvider->setNoDataValue( 1, mNoData );
183 nonGroundDestProvider->setEditable( true );
184 }
185
188 const int numBlocksX = static_cast< int>( std::ceil( 1.0 * mLayerWidth / blockWidth ) );
189 const int numBlocksY = static_cast< int >( std::ceil( 1.0 * mLayerHeight / blockHeight ) );
190 const int numBlocks = numBlocksX * numBlocksY;
191
192 const int radius = parameterAsInt( parameters, QStringLiteral( "RADIUS" ), context );
193
194 const double terrainSlopePercent = parameterAsDouble( parameters, QStringLiteral( "TERRAIN_SLOPE" ), context ) / 100; //20.0 / 100 * 0.143;
195 const int filterModification = parameterAsEnum( parameters, QStringLiteral( "FILTER_MODIFICATION" ), context );
196 const double standardDeviation = parameterAsDouble( parameters, QStringLiteral( "STANDARD_DEVIATION" ), context );
197
198 // create kernel
199 QVector<double> kernel;
200 kernel.reserve( ( radius * 2 ) * ( radius * 2 ) );
201 int kernelSize = 0;
202 for ( int y = -radius; y <= radius; y++ )
203 {
204 for ( int x = -radius; x <= radius; x++ )
205 {
206 const double distance = std::sqrt( x * x + y * y );
207 if ( distance < radius )
208 {
209 kernelSize++;
210 kernel.push_back( x );
211 kernel.push_back( y );
212 switch ( filterModification )
213 {
214 case 0:
215 kernel.push_back( distance * terrainSlopePercent );
216 break;
217
218 case 1:
219 kernel.push_back( distance * terrainSlopePercent + 1.65 * std::sqrt( 2 * standardDeviation ) );
220 break;
221
222 case 2:
223 {
224 const double dz = distance * terrainSlopePercent - 1.65 * std::sqrt( 2 * standardDeviation );
225 kernel.push_back( dz > 0 ? dz : 0 );
226 break;
227 }
228 }
229 }
230 }
231 }
232
233 QgsRasterIterator iter( mInterface.get(), radius );
234 iter.startRasterRead( 1, mLayerWidth, mLayerHeight, mExtent );
235 int iterLeft = 0;
236 int iterTop = 0;
237 int iterCols = 0;
238 int iterRows = 0;
239 int tileLeft = 0;
240 int tileTop = 0;
241 int tileCols = 0;
242 int tileRows = 0;
243
244 QgsRectangle blockExtent;
245
246 std::unique_ptr< QgsRasterBlock > inputBlock;
247 int blockIndex = 0;
248 while ( iter.readNextRasterPart( 1, iterCols, iterRows, inputBlock, iterLeft, iterTop, &blockExtent, &tileCols, &tileRows, &tileLeft, &tileTop ) )
249 {
250 std::unique_ptr< QgsRasterBlock > outputGroundBlock;
251 if ( groundDestProvider )
252 outputGroundBlock = std::make_unique< QgsRasterBlock >( mDataType, tileCols, tileRows );
253
254 std::unique_ptr< QgsRasterBlock > outputNonGroundBlock;
255 if ( nonGroundDestProvider )
256 outputNonGroundBlock = std::make_unique< QgsRasterBlock >( mDataType, tileCols, tileRows );
257
258 double baseProgress = static_cast< double >( blockIndex ) / numBlocks;
259 feedback->setProgress( 100.0 * baseProgress );
260 blockIndex++;
261 if ( feedback->isCanceled() )
262 break;
263
264 const int tileBoundaryLeft = tileLeft - iterLeft;
265 const int tileBoundaryTop = tileTop - iterTop;
266
267 const double rowProgressStep = 1.0 / numBlocks / tileRows;
268 double rowProgress = 0;
269 for ( int row = tileBoundaryTop; row < tileBoundaryTop + tileRows; row++ )
270 {
271 if ( feedback->isCanceled() )
272 break;
273
274 feedback->setProgress( 100.0 * ( baseProgress + rowProgress ) );
275 rowProgress += rowProgressStep;
276
277 for ( int col = tileBoundaryLeft; col < tileBoundaryLeft + tileCols; col++ )
278 {
279 if ( feedback->isCanceled() )
280 break;
281
282 bool isNoData = false;
283 const double val = inputBlock->valueAndNoData( row, col, isNoData );
284 if ( isNoData )
285 {
286 if ( outputGroundBlock )
287 outputGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
288 if ( outputNonGroundBlock )
289 outputNonGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
290 }
291 else
292 {
293 bool nonGround = false;
294 const double *kernelData = kernel.constData();
295 for ( int i = 0; i < kernelSize; ++i )
296 {
297 const int dx = static_cast< int >( *kernelData++ );
298 const int dy = static_cast< int >( *kernelData++ );
299 const double distance = *kernelData++;
300 const int rCol = col + dx;
301 const int rRow = row + dy;
302 if ( rCol >= 0 && ( rCol < ( iterLeft + iterCols ) ) && rRow >= 0 && ( rRow < ( iterTop + iterRows ) ) )
303 {
304 bool otherIsNoData = false;
305 const double otherVal = inputBlock->valueAndNoData( rRow, rCol, otherIsNoData );
306 if ( !otherIsNoData )
307 {
308 const double dz = val - otherVal;
309 if ( dz > 0 && dz > distance )
310 {
311 nonGround = true;
312 break;
313 }
314 }
315 }
316 }
317 if ( nonGround )
318 {
319 if ( outputGroundBlock )
320 outputGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
321 if ( outputNonGroundBlock )
322 outputNonGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, val );
323 }
324 else
325 {
326 if ( outputGroundBlock )
327 outputGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, val );
328 if ( outputNonGroundBlock )
329 outputNonGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
330 }
331 }
332 }
333 }
334 if ( groundDestProvider )
335 groundDestProvider->writeBlock( outputGroundBlock.get(), mBand, tileLeft, tileTop );
336 if ( nonGroundDestProvider )
337 nonGroundDestProvider->writeBlock( outputNonGroundBlock.get(), mBand, tileLeft, tileTop );
338 }
339 if ( groundDestProvider )
340 groundDestProvider->setEditable( false );
341 if ( nonGroundDestProvider )
342 nonGroundDestProvider->setEditable( false );
343
344 QVariantMap outputs;
345 outputs.insert( QStringLiteral( "OUTPUT_GROUND" ), groundOutputFile );
346 outputs.insert( QStringLiteral( "OUTPUT_NONGROUND" ), nonGroundOutputFile );
347 return outputs;
348}
349
350
352
353
354
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition qgsfeedback.h:53
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition qgsfeedback.h:61
virtual QgsRectangle extent() const
Returns the extent of the layer.
QgsCoordinateReferenceSystem crs
Definition qgsmaplayer.h:82
Contains information about the context in which a processing algorithm is executed.
Custom exception class for processing related exceptions.
Base class for providing feedback from a processing algorithm.
A raster band parameter for Processing algorithms.
A raster layer parameter for processing algorithms.
QgsRasterDataProvider * clone() const override=0
Clone itself, create deep copy.
virtual bool sourceHasNoDataValue(int bandNo) const
Returns true if source band has no data value.
virtual double sourceNoDataValue(int bandNo) const
Value representing no data value.
Qgis::DataType dataType(int bandNo) const override=0
Returns data type for the band specified by number.
static QString driverForExtension(const QString &extension)
Returns the GDAL driver name for a specified file extension.
Iterator for sequentially processing raster cells.
static const int DEFAULT_MAXIMUM_TILE_WIDTH
Default maximum tile width.
static const int DEFAULT_MAXIMUM_TILE_HEIGHT
Default maximum tile height.
Represents a raster layer.
int height() const
Returns the height of the (unclipped) raster.
int bandCount() const
Returns the number of bands in this layer.
double rasterUnitsPerPixelX() const
Returns the number of raster units per each raster pixel in X axis.
QgsRasterDataProvider * dataProvider() override
Returns the source data provider.
double rasterUnitsPerPixelY() const
Returns the number of raster units per each raster pixel in Y axis.
int width() const
Returns the width of the (unclipped) raster.
A rectangle specified with double values.