24QString QgsRasterDtmSlopeBasedFilterAlgorithm::name()
const
26 return QStringLiteral(
"dtmslopebasedfilter" );
29QString QgsRasterDtmSlopeBasedFilterAlgorithm::displayName()
const
31 return QObject::tr(
"DTM filter (slope-based)" );
34QStringList QgsRasterDtmSlopeBasedFilterAlgorithm::tags()
const
36 return QObject::tr(
"dem,filter,slope,dsm,dtm,terrain" ).split(
',' );
39QString QgsRasterDtmSlopeBasedFilterAlgorithm::group()
const
41 return QObject::tr(
"Raster terrain analysis" );
44QString QgsRasterDtmSlopeBasedFilterAlgorithm::groupId()
const
46 return QStringLiteral(
"rasterterrainanalysis" );
49QString QgsRasterDtmSlopeBasedFilterAlgorithm::shortHelpString()
const
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." );
66QString QgsRasterDtmSlopeBasedFilterAlgorithm::shortDescription()
const
68 return QObject::tr(
"Filters a Digital Elevation Model in order to classify its cells into ground and object (non-ground) cells." );
71void QgsRasterDtmSlopeBasedFilterAlgorithm::initAlgorithm(
const QVariantMap & )
75 addParameter(
new QgsProcessingParameterBand( QStringLiteral(
"BAND" ), QObject::tr(
"Band number" ), 1, QStringLiteral(
"INPUT" ) ) );
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() );
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() );
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() );
90 stDevParam->setHelp( QObject::tr(
"The standard deviation used to calculate a 5% confidence interval applied to the height threshold." ) );
91 addParameter( stDevParam.release() );
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" ) } } ) } } ) );
98 addParameter( createOptsParam.release() );
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" ) } } ) } } ) );
103 addParameter( creationOptsParam.release() );
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() );
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() );
114QgsRasterDtmSlopeBasedFilterAlgorithm *QgsRasterDtmSlopeBasedFilterAlgorithm::createInstance()
const
116 return new QgsRasterDtmSlopeBasedFilterAlgorithm();
121 QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral(
"INPUT" ), context );
125 const int band = parameterAsInt( parameters, QStringLiteral(
"BAND" ), context );
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() ) );
133 mLayerWidth = layer->
width();
134 mLayerHeight = layer->
height();
135 mExtent = layer->
extent();
146 QString creationOptions = parameterAsString( parameters, QStringLiteral(
"CREATION_OPTIONS" ), context ).trimmed();
148 const QString optionsString = parameterAsString( parameters, QStringLiteral(
"CREATE_OPTIONS" ), context );
149 if ( !optionsString.isEmpty() )
150 creationOptions = optionsString;
152 const QString groundOutputFile = parameterAsOutputLayer( parameters, QStringLiteral(
"OUTPUT_GROUND" ), context );
153 std::unique_ptr<QgsRasterFileWriter> groundWriter;
154 std::unique_ptr<QgsRasterDataProvider> groundDestProvider;
156 if ( !groundOutputFile.isEmpty() )
158 const QFileInfo fi( groundOutputFile );
161 groundWriter = std::make_unique<QgsRasterFileWriter>( groundOutputFile );
162 groundWriter->setOutputProviderKey( QStringLiteral(
"gdal" ) );
163 if ( !creationOptions.isEmpty() )
165 groundWriter->setCreationOptions( creationOptions.split(
'|' ) );
167 groundWriter->setOutputFormat( outputFormat );
169 groundDestProvider.reset( groundWriter->createOneBandRaster( mDataType, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
171 if ( !groundDestProvider )
172 throw QgsProcessingException( QObject::tr(
"Could not create raster output: %1" ).arg( groundOutputFile ) );
173 if ( !groundDestProvider->isValid() )
176 groundDestProvider->setNoDataValue( 1, mNoData );
177 groundDestProvider->setEditable(
true );
180 const QString nonGroundOutputFile = parameterAsOutputLayer( parameters, QStringLiteral(
"OUTPUT_NONGROUND" ), context );
181 std::unique_ptr<QgsRasterFileWriter> nonGroundWriter;
182 std::unique_ptr<QgsRasterDataProvider> nonGroundDestProvider;
184 if ( !nonGroundOutputFile.isEmpty() )
186 const QFileInfo fi( groundOutputFile );
189 nonGroundWriter = std::make_unique<QgsRasterFileWriter>( nonGroundOutputFile );
190 nonGroundWriter->setOutputProviderKey( QStringLiteral(
"gdal" ) );
191 if ( !creationOptions.isEmpty() )
193 nonGroundWriter->setCreationOptions( creationOptions.split(
'|' ) );
195 nonGroundWriter->setOutputFormat( outputFormat );
197 nonGroundDestProvider.reset( nonGroundWriter->createOneBandRaster( mDataType, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
199 if ( !nonGroundDestProvider )
200 throw QgsProcessingException( QObject::tr(
"Could not create raster output: %1" ).arg( nonGroundOutputFile ) );
201 if ( !nonGroundDestProvider->isValid() )
204 nonGroundDestProvider->setNoDataValue( 1, mNoData );
205 nonGroundDestProvider->setEditable(
true );
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;
214 const int radius = parameterAsInt( parameters, QStringLiteral(
"RADIUS" ), context );
216 const double terrainSlopePercent = parameterAsDouble( parameters, QStringLiteral(
"TERRAIN_SLOPE" ), context ) / 100;
217 const int filterModification = parameterAsEnum( parameters, QStringLiteral(
"FILTER_MODIFICATION" ), context );
218 const double standardDeviation = parameterAsDouble( parameters, QStringLiteral(
"STANDARD_DEVIATION" ), context );
221 QVector<double> kernel;
222 kernel.reserve( ( radius * 2 ) * ( radius * 2 ) );
224 for (
int y = -radius; y <= radius; y++ )
226 for (
int x = -radius; x <= radius; x++ )
228 const double distance = std::sqrt( x * x + y * y );
229 if ( distance < radius )
232 kernel.push_back( x );
233 kernel.push_back( y );
234 switch ( filterModification )
237 kernel.push_back( distance * terrainSlopePercent );
241 kernel.push_back( distance * terrainSlopePercent + 1.65 * std::sqrt( 2 * standardDeviation ) );
246 const double dz = distance * terrainSlopePercent - 1.65 * std::sqrt( 2 * standardDeviation );
247 kernel.push_back( dz > 0 ? dz : 0 );
256 iter.startRasterRead( 1, mLayerWidth, mLayerHeight, mExtent );
268 std::unique_ptr<QgsRasterBlock> inputBlock;
270 while ( iter.readNextRasterPart( 1, iterCols, iterRows, inputBlock, iterLeft, iterTop, &blockExtent, &tileCols, &tileRows, &tileLeft, &tileTop ) )
272 std::unique_ptr<QgsRasterBlock> outputGroundBlock;
273 if ( groundDestProvider )
274 outputGroundBlock = std::make_unique<QgsRasterBlock>( mDataType, tileCols, tileRows );
276 std::unique_ptr<QgsRasterBlock> outputNonGroundBlock;
277 if ( nonGroundDestProvider )
278 outputNonGroundBlock = std::make_unique<QgsRasterBlock>( mDataType, tileCols, tileRows );
280 double baseProgress =
static_cast<double>( blockIndex ) / numBlocks;
286 const int tileBoundaryLeft = tileLeft - iterLeft;
287 const int tileBoundaryTop = tileTop - iterTop;
289 const double rowProgressStep = 1.0 / numBlocks / tileRows;
290 double rowProgress = 0;
291 for (
int row = tileBoundaryTop; row < tileBoundaryTop + tileRows; row++ )
296 feedback->
setProgress( 100.0 * ( baseProgress + rowProgress ) );
297 rowProgress += rowProgressStep;
299 for (
int col = tileBoundaryLeft; col < tileBoundaryLeft + tileCols; col++ )
304 bool isNoData =
false;
305 const double val = inputBlock->valueAndNoData( row, col, isNoData );
308 if ( outputGroundBlock )
309 outputGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
310 if ( outputNonGroundBlock )
311 outputNonGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
315 bool nonGround =
false;
316 const double *kernelData = kernel.constData();
317 for (
int i = 0; i < kernelSize; ++i )
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 ) ) )
326 bool otherIsNoData =
false;
327 const double otherVal = inputBlock->valueAndNoData( rRow, rCol, otherIsNoData );
328 if ( !otherIsNoData )
330 const double dz = val - otherVal;
331 if ( dz > 0 && dz > distance )
341 if ( outputGroundBlock )
342 outputGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
343 if ( outputNonGroundBlock )
344 outputNonGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, val );
348 if ( outputGroundBlock )
349 outputGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, val );
350 if ( outputNonGroundBlock )
351 outputNonGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
356 if ( groundDestProvider )
358 if ( !groundDestProvider->writeBlock( outputGroundBlock.get(), mBand, tileLeft, tileTop ) )
360 throw QgsProcessingException( QObject::tr(
"Could not write raster block: %1" ).arg( groundDestProvider->error().summary() ) );
363 if ( nonGroundDestProvider )
365 if ( !nonGroundDestProvider->writeBlock( outputNonGroundBlock.get(), mBand, tileLeft, tileTop ) )
367 throw QgsProcessingException( QObject::tr(
"Could not write raster block: %1" ).arg( nonGroundDestProvider->error().summary() ) );
371 if ( groundDestProvider )
372 groundDestProvider->setEditable(
false );
373 if ( nonGroundDestProvider )
374 nonGroundDestProvider->setEditable(
false );
377 outputs.insert( QStringLiteral(
"OUTPUT_GROUND" ), groundOutputFile );
378 outputs.insert( QStringLiteral(
"OUTPUT_NONGROUND" ), nonGroundOutputFile );
@ 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.
@ Double
Double/float values.
bool isCanceled() const
Tells whether the operation has been canceled already.
void setProgress(double progress)
Sets the current progress for the feedback object.
virtual QgsRectangle extent() const
Returns the extent of the layer.
QgsCoordinateReferenceSystem crs
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.