QGIS API Documentation 3.43.0-Master (e01d6d7c4c0)
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
66QString QgsRasterDtmSlopeBasedFilterAlgorithm::shortDescription() const
67{
68 return QObject::tr( "Filters a Digital Elevation Model in order to classify its cells into ground and object (non-ground) cells." );
69}
70
71void QgsRasterDtmSlopeBasedFilterAlgorithm::initAlgorithm( const QVariantMap & )
72{
73 addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT" ), QObject::tr( "Input layer" ) ) );
74
75 addParameter( new QgsProcessingParameterBand( QStringLiteral( "BAND" ), QObject::tr( "Band number" ), 1, QStringLiteral( "INPUT" ) ) );
76
77 auto radiusParam = std::make_unique<QgsProcessingParameterNumber>( QStringLiteral( "RADIUS" ), QObject::tr( "Kernel radius (pixels)" ), Qgis::ProcessingNumberParameterType::Integer, 5, false, 1, 1000 );
78 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." ) );
79 addParameter( radiusParam.release() );
80
81 auto terrainSlopeParam = std::make_unique<QgsProcessingParameterNumber>( QStringLiteral( "TERRAIN_SLOPE" ), QObject::tr( "Terrain slope (%, pixel size/vertical units)" ), Qgis::ProcessingNumberParameterType::Double, 30, false, 0, 1000 );
82 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." ) );
83 addParameter( terrainSlopeParam.release() );
84
85 auto filterModificationParam = std::make_unique<QgsProcessingParameterEnum>( QStringLiteral( "FILTER_MODIFICATION" ), 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 auto stDevParam = std::make_unique<QgsProcessingParameterNumber>( QStringLiteral( "STANDARD_DEVIATION" ), QObject::tr( "Standard deviation" ), Qgis::ProcessingNumberParameterType::Double, 0.1, false, 0, 1000 );
90 stDevParam->setHelp( QObject::tr( "The standard deviation used to calculate a 5% confidence interval applied to the height threshold." ) );
91 addParameter( stDevParam.release() );
92
93 // backwards compatibility parameter
94 // TODO QGIS 4: remove parameter and related logic
95 auto createOptsParam = std::make_unique<QgsProcessingParameterString>( QStringLiteral( "CREATE_OPTIONS" ), QObject::tr( "Creation options" ), QVariant(), false, true );
96 createOptsParam->setMetadata( QVariantMap( { { QStringLiteral( "widget_wrapper" ), QVariantMap( { { QStringLiteral( "widget_type" ), QStringLiteral( "rasteroptions" ) } } ) } } ) );
97 createOptsParam->setFlags( createOptsParam->flags() | Qgis::ProcessingParameterFlag::Hidden );
98 addParameter( createOptsParam.release() );
99
100 auto creationOptsParam = std::make_unique<QgsProcessingParameterString>( QStringLiteral( "CREATION_OPTIONS" ), QObject::tr( "Creation options" ), QVariant(), false, true );
101 creationOptsParam->setMetadata( QVariantMap( { { QStringLiteral( "widget_wrapper" ), QVariantMap( { { QStringLiteral( "widget_type" ), QStringLiteral( "rasteroptions" ) } } ) } } ) );
102 creationOptsParam->setFlags( creationOptsParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
103 addParameter( creationOptsParam.release() );
104
105 auto outputLayerGroundParam = std::make_unique<QgsProcessingParameterRasterDestination>( QStringLiteral( "OUTPUT_GROUND" ), QObject::tr( "Output layer (ground)" ), QVariant(), true, true );
106 outputLayerGroundParam->setHelp( QObject::tr( "The filtered DEM containing only cells classified as ground." ) );
107 addParameter( outputLayerGroundParam.release() );
108
109 auto outputLayerNonGroundParam = std::make_unique<QgsProcessingParameterRasterDestination>( QStringLiteral( "OUTPUT_NONGROUND" ), QObject::tr( "Output layer (non-ground objects)" ), QVariant(), true, false );
110 outputLayerNonGroundParam->setHelp( QObject::tr( "The non-ground objects removed by the filter." ) );
111 addParameter( outputLayerNonGroundParam.release() );
112}
113
114QgsRasterDtmSlopeBasedFilterAlgorithm *QgsRasterDtmSlopeBasedFilterAlgorithm::createInstance() const
115{
116 return new QgsRasterDtmSlopeBasedFilterAlgorithm();
117}
118
119bool QgsRasterDtmSlopeBasedFilterAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
120{
121 QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT" ), context );
122 if ( !layer )
123 throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT" ) ) );
124
125 const int band = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
126
127 mBand = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
128 if ( mBand < 1 || mBand > layer->bandCount() )
129 throw QgsProcessingException( QObject::tr( "Invalid band number for BAND (%1): Valid values for input raster are 1 to %2" ).arg( mBand ).arg( layer->bandCount() ) );
130
131 mInterface.reset( layer->dataProvider()->clone() );
132 mHasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( band );
133 mLayerWidth = layer->width();
134 mLayerHeight = layer->height();
135 mExtent = layer->extent();
136 mCrs = layer->crs();
137 mRasterUnitsPerPixelX = layer->rasterUnitsPerPixelX();
138 mRasterUnitsPerPixelY = layer->rasterUnitsPerPixelY();
139 mDataType = layer->dataProvider()->dataType( mBand );
140 mNoData = layer->dataProvider()->sourceNoDataValue( mBand );
141 return true;
142}
143
144QVariantMap QgsRasterDtmSlopeBasedFilterAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
145{
146 QString creationOptions = parameterAsString( parameters, QStringLiteral( "CREATION_OPTIONS" ), context ).trimmed();
147 // handle backwards compatibility parameter CREATE_OPTIONS
148 const QString optionsString = parameterAsString( parameters, QStringLiteral( "CREATE_OPTIONS" ), context );
149 if ( !optionsString.isEmpty() )
150 creationOptions = optionsString;
151
152 const QString groundOutputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT_GROUND" ), context );
153 std::unique_ptr<QgsRasterFileWriter> groundWriter;
154 std::unique_ptr<QgsRasterDataProvider> groundDestProvider;
155
156 if ( !groundOutputFile.isEmpty() )
157 {
158 const QFileInfo fi( groundOutputFile );
159 const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
160
161 groundWriter = std::make_unique<QgsRasterFileWriter>( groundOutputFile );
162 groundWriter->setOutputProviderKey( QStringLiteral( "gdal" ) );
163 if ( !creationOptions.isEmpty() )
164 {
165 groundWriter->setCreationOptions( creationOptions.split( '|' ) );
166 }
167 groundWriter->setOutputFormat( outputFormat );
168
169 groundDestProvider.reset( groundWriter->createOneBandRaster( mDataType, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
170
171 if ( !groundDestProvider )
172 throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( groundOutputFile ) );
173 if ( !groundDestProvider->isValid() )
174 throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( groundOutputFile, groundDestProvider->error().message( QgsErrorMessage::Text ) ) );
175
176 groundDestProvider->setNoDataValue( 1, mNoData );
177 groundDestProvider->setEditable( true );
178 }
179
180 const QString nonGroundOutputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT_NONGROUND" ), context );
181 std::unique_ptr<QgsRasterFileWriter> nonGroundWriter;
182 std::unique_ptr<QgsRasterDataProvider> nonGroundDestProvider;
183
184 if ( !nonGroundOutputFile.isEmpty() )
185 {
186 const QFileInfo fi( groundOutputFile );
187 const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
188
189 nonGroundWriter = std::make_unique<QgsRasterFileWriter>( nonGroundOutputFile );
190 nonGroundWriter->setOutputProviderKey( QStringLiteral( "gdal" ) );
191 if ( !creationOptions.isEmpty() )
192 {
193 nonGroundWriter->setCreationOptions( creationOptions.split( '|' ) );
194 }
195 nonGroundWriter->setOutputFormat( outputFormat );
196
197 nonGroundDestProvider.reset( nonGroundWriter->createOneBandRaster( mDataType, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
198
199 if ( !nonGroundDestProvider )
200 throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( nonGroundOutputFile ) );
201 if ( !nonGroundDestProvider->isValid() )
202 throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( nonGroundOutputFile, nonGroundDestProvider->error().message( QgsErrorMessage::Text ) ) );
203
204 nonGroundDestProvider->setNoDataValue( 1, mNoData );
205 nonGroundDestProvider->setEditable( true );
206 }
207
210 const int numBlocksX = static_cast<int>( std::ceil( 1.0 * mLayerWidth / blockWidth ) );
211 const int numBlocksY = static_cast<int>( std::ceil( 1.0 * mLayerHeight / blockHeight ) );
212 const int numBlocks = numBlocksX * numBlocksY;
213
214 const int radius = parameterAsInt( parameters, QStringLiteral( "RADIUS" ), context );
215
216 const double terrainSlopePercent = parameterAsDouble( parameters, QStringLiteral( "TERRAIN_SLOPE" ), context ) / 100; //20.0 / 100 * 0.143;
217 const int filterModification = parameterAsEnum( parameters, QStringLiteral( "FILTER_MODIFICATION" ), context );
218 const double standardDeviation = parameterAsDouble( parameters, QStringLiteral( "STANDARD_DEVIATION" ), context );
219
220 // create kernel
221 QVector<double> kernel;
222 kernel.reserve( ( radius * 2 ) * ( radius * 2 ) );
223 int kernelSize = 0;
224 for ( int y = -radius; y <= radius; y++ )
225 {
226 for ( int x = -radius; x <= radius; x++ )
227 {
228 const double distance = std::sqrt( x * x + y * y );
229 if ( distance < radius )
230 {
231 kernelSize++;
232 kernel.push_back( x );
233 kernel.push_back( y );
234 switch ( filterModification )
235 {
236 case 0:
237 kernel.push_back( distance * terrainSlopePercent );
238 break;
239
240 case 1:
241 kernel.push_back( distance * terrainSlopePercent + 1.65 * std::sqrt( 2 * standardDeviation ) );
242 break;
243
244 case 2:
245 {
246 const double dz = distance * terrainSlopePercent - 1.65 * std::sqrt( 2 * standardDeviation );
247 kernel.push_back( dz > 0 ? dz : 0 );
248 break;
249 }
250 }
251 }
252 }
253 }
254
255 QgsRasterIterator iter( mInterface.get(), radius );
256 iter.startRasterRead( 1, mLayerWidth, mLayerHeight, mExtent );
257 int iterLeft = 0;
258 int iterTop = 0;
259 int iterCols = 0;
260 int iterRows = 0;
261 int tileLeft = 0;
262 int tileTop = 0;
263 int tileCols = 0;
264 int tileRows = 0;
265
266 QgsRectangle blockExtent;
267
268 std::unique_ptr<QgsRasterBlock> inputBlock;
269 int blockIndex = 0;
270 while ( iter.readNextRasterPart( 1, iterCols, iterRows, inputBlock, iterLeft, iterTop, &blockExtent, &tileCols, &tileRows, &tileLeft, &tileTop ) )
271 {
272 std::unique_ptr<QgsRasterBlock> outputGroundBlock;
273 if ( groundDestProvider )
274 outputGroundBlock = std::make_unique<QgsRasterBlock>( mDataType, tileCols, tileRows );
275
276 std::unique_ptr<QgsRasterBlock> outputNonGroundBlock;
277 if ( nonGroundDestProvider )
278 outputNonGroundBlock = std::make_unique<QgsRasterBlock>( mDataType, tileCols, tileRows );
279
280 double baseProgress = static_cast<double>( blockIndex ) / numBlocks;
281 feedback->setProgress( 100.0 * baseProgress );
282 blockIndex++;
283 if ( feedback->isCanceled() )
284 break;
285
286 const int tileBoundaryLeft = tileLeft - iterLeft;
287 const int tileBoundaryTop = tileTop - iterTop;
288
289 const double rowProgressStep = 1.0 / numBlocks / tileRows;
290 double rowProgress = 0;
291 for ( int row = tileBoundaryTop; row < tileBoundaryTop + tileRows; row++ )
292 {
293 if ( feedback->isCanceled() )
294 break;
295
296 feedback->setProgress( 100.0 * ( baseProgress + rowProgress ) );
297 rowProgress += rowProgressStep;
298
299 for ( int col = tileBoundaryLeft; col < tileBoundaryLeft + tileCols; col++ )
300 {
301 if ( feedback->isCanceled() )
302 break;
303
304 bool isNoData = false;
305 const double val = inputBlock->valueAndNoData( row, col, isNoData );
306 if ( isNoData )
307 {
308 if ( outputGroundBlock )
309 outputGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
310 if ( outputNonGroundBlock )
311 outputNonGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
312 }
313 else
314 {
315 bool nonGround = false;
316 const double *kernelData = kernel.constData();
317 for ( int i = 0; i < kernelSize; ++i )
318 {
319 const int dx = static_cast<int>( *kernelData++ );
320 const int dy = static_cast<int>( *kernelData++ );
321 const double distance = *kernelData++;
322 const int rCol = col + dx;
323 const int rRow = row + dy;
324 if ( rCol >= 0 && ( rCol < ( iterLeft + iterCols ) ) && rRow >= 0 && ( rRow < ( iterTop + iterRows ) ) )
325 {
326 bool otherIsNoData = false;
327 const double otherVal = inputBlock->valueAndNoData( rRow, rCol, otherIsNoData );
328 if ( !otherIsNoData )
329 {
330 const double dz = val - otherVal;
331 if ( dz > 0 && dz > distance )
332 {
333 nonGround = true;
334 break;
335 }
336 }
337 }
338 }
339 if ( nonGround )
340 {
341 if ( outputGroundBlock )
342 outputGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
343 if ( outputNonGroundBlock )
344 outputNonGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, val );
345 }
346 else
347 {
348 if ( outputGroundBlock )
349 outputGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, val );
350 if ( outputNonGroundBlock )
351 outputNonGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
352 }
353 }
354 }
355 }
356 if ( groundDestProvider )
357 {
358 if ( !groundDestProvider->writeBlock( outputGroundBlock.get(), mBand, tileLeft, tileTop ) )
359 {
360 throw QgsProcessingException( QObject::tr( "Could not write raster block: %1" ).arg( groundDestProvider->error().summary() ) );
361 }
362 }
363 if ( nonGroundDestProvider )
364 {
365 if ( !nonGroundDestProvider->writeBlock( outputNonGroundBlock.get(), mBand, tileLeft, tileTop ) )
366 {
367 throw QgsProcessingException( QObject::tr( "Could not write raster block: %1" ).arg( nonGroundDestProvider->error().summary() ) );
368 }
369 }
370 }
371 if ( groundDestProvider )
372 groundDestProvider->setEditable( false );
373 if ( nonGroundDestProvider )
374 nonGroundDestProvider->setEditable( false );
375
376 QVariantMap outputs;
377 outputs.insert( QStringLiteral( "OUTPUT_GROUND" ), groundOutputFile );
378 outputs.insert( QStringLiteral( "OUTPUT_NONGROUND" ), nonGroundOutputFile );
379 return outputs;
380}
381
382
@ Hidden
Parameter is hidden and should not be shown to users.
@ Advanced
Parameter is an advanced parameter which should be hidden from users by default.
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:84
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.