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." );
66void QgsRasterDtmSlopeBasedFilterAlgorithm::initAlgorithm(
const QVariantMap & )
69 QObject::tr(
"Input layer" ) ) );
72 QObject::tr(
"Band number" ), 1, QStringLiteral(
"INPUT" ) ) );
74 std::unique_ptr< QgsProcessingParameterNumber > radiusParam = std::make_unique< QgsProcessingParameterNumber >( QStringLiteral(
"RADIUS" ),
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() );
79 std::unique_ptr< QgsProcessingParameterNumber > terrainSlopeParam = std::make_unique< QgsProcessingParameterNumber >( QStringLiteral(
"TERRAIN_SLOPE" ),
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() );
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() );
89 std::unique_ptr< QgsProcessingParameterNumber > stDevParam = std::make_unique< QgsProcessingParameterNumber >( QStringLiteral(
"STANDARD_DEVIATION" ),
91 stDevParam->setHelp( QObject::tr(
"The standard deviation used to calculate a 5% confidence interval applied to the height threshold." ) );
92 addParameter( stDevParam.release() );
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() );
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() );
105QgsRasterDtmSlopeBasedFilterAlgorithm *QgsRasterDtmSlopeBasedFilterAlgorithm::createInstance()
const
107 return new QgsRasterDtmSlopeBasedFilterAlgorithm();
112 QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral(
"INPUT" ), context );
116 const int band = parameterAsInt( parameters, QStringLiteral(
"BAND" ), context );
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 )
125 mLayerWidth = layer->
width();
126 mLayerHeight = layer->
height();
127 mExtent = layer->
extent();
138 const QString groundOutputFile = parameterAsOutputLayer( parameters, QStringLiteral(
"OUTPUT_GROUND" ), context );
139 std::unique_ptr< QgsRasterFileWriter > groundWriter;
140 std::unique_ptr<QgsRasterDataProvider > groundDestProvider;
142 if ( !groundOutputFile.isEmpty() )
144 const QFileInfo fi( groundOutputFile );
147 groundWriter = std::make_unique< QgsRasterFileWriter >( groundOutputFile );
148 groundWriter->setOutputProviderKey( QStringLiteral(
"gdal" ) );
149 groundWriter->setOutputFormat( outputFormat );
151 groundDestProvider.reset( groundWriter->createOneBandRaster( mDataType, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
153 if ( !groundDestProvider )
154 throw QgsProcessingException( QObject::tr(
"Could not create raster output: %1" ).arg( groundOutputFile ) );
155 if ( !groundDestProvider->isValid() )
158 groundDestProvider->setNoDataValue( 1, mNoData );
159 groundDestProvider->setEditable(
true );
162 const QString nonGroundOutputFile = parameterAsOutputLayer( parameters, QStringLiteral(
"OUTPUT_NONGROUND" ), context );
163 std::unique_ptr< QgsRasterFileWriter > nonGroundWriter;
164 std::unique_ptr<QgsRasterDataProvider > nonGroundDestProvider;
166 if ( !nonGroundOutputFile.isEmpty() )
168 const QFileInfo fi( groundOutputFile );
171 nonGroundWriter = std::make_unique< QgsRasterFileWriter >( nonGroundOutputFile );
172 nonGroundWriter->setOutputProviderKey( QStringLiteral(
"gdal" ) );
173 nonGroundWriter->setOutputFormat( outputFormat );
175 nonGroundDestProvider.reset( nonGroundWriter->createOneBandRaster( mDataType, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
177 if ( !nonGroundDestProvider )
178 throw QgsProcessingException( QObject::tr(
"Could not create raster output: %1" ).arg( nonGroundOutputFile ) );
179 if ( !nonGroundDestProvider->isValid() )
182 nonGroundDestProvider->setNoDataValue( 1, mNoData );
183 nonGroundDestProvider->setEditable(
true );
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;
192 const int radius = parameterAsInt( parameters, QStringLiteral(
"RADIUS" ), context );
194 const double terrainSlopePercent = parameterAsDouble( parameters, QStringLiteral(
"TERRAIN_SLOPE" ), context ) / 100;
195 const int filterModification = parameterAsEnum( parameters, QStringLiteral(
"FILTER_MODIFICATION" ), context );
196 const double standardDeviation = parameterAsDouble( parameters, QStringLiteral(
"STANDARD_DEVIATION" ), context );
199 QVector<double> kernel;
200 kernel.reserve( ( radius * 2 ) * ( radius * 2 ) );
202 for (
int y = -radius; y <= radius; y++ )
204 for (
int x = -radius; x <= radius; x++ )
206 const double distance = std::sqrt( x * x + y * y );
207 if ( distance < radius )
210 kernel.push_back( x );
211 kernel.push_back( y );
212 switch ( filterModification )
215 kernel.push_back( distance * terrainSlopePercent );
219 kernel.push_back( distance * terrainSlopePercent + 1.65 * std::sqrt( 2 * standardDeviation ) );
224 const double dz = distance * terrainSlopePercent - 1.65 * std::sqrt( 2 * standardDeviation );
225 kernel.push_back( dz > 0 ? dz : 0 );
234 iter.startRasterRead( 1, mLayerWidth, mLayerHeight, mExtent );
246 std::unique_ptr< QgsRasterBlock > inputBlock;
248 while ( iter.readNextRasterPart( 1, iterCols, iterRows, inputBlock, iterLeft, iterTop, &blockExtent, &tileCols, &tileRows, &tileLeft, &tileTop ) )
250 std::unique_ptr< QgsRasterBlock > outputGroundBlock;
251 if ( groundDestProvider )
252 outputGroundBlock = std::make_unique< QgsRasterBlock >( mDataType, tileCols, tileRows );
254 std::unique_ptr< QgsRasterBlock > outputNonGroundBlock;
255 if ( nonGroundDestProvider )
256 outputNonGroundBlock = std::make_unique< QgsRasterBlock >( mDataType, tileCols, tileRows );
258 double baseProgress =
static_cast< double >( blockIndex ) / numBlocks;
264 const int tileBoundaryLeft = tileLeft - iterLeft;
265 const int tileBoundaryTop = tileTop - iterTop;
267 const double rowProgressStep = 1.0 / numBlocks / tileRows;
268 double rowProgress = 0;
269 for (
int row = tileBoundaryTop; row < tileBoundaryTop + tileRows; row++ )
274 feedback->
setProgress( 100.0 * ( baseProgress + rowProgress ) );
275 rowProgress += rowProgressStep;
277 for (
int col = tileBoundaryLeft; col < tileBoundaryLeft + tileCols; col++ )
282 bool isNoData =
false;
283 const double val = inputBlock->valueAndNoData( row, col, isNoData );
286 if ( outputGroundBlock )
287 outputGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
288 if ( outputNonGroundBlock )
289 outputNonGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
293 bool nonGround =
false;
294 const double *kernelData = kernel.constData();
295 for (
int i = 0; i < kernelSize; ++i )
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 ) ) )
304 bool otherIsNoData =
false;
305 const double otherVal = inputBlock->valueAndNoData( rRow, rCol, otherIsNoData );
306 if ( !otherIsNoData )
308 const double dz = val - otherVal;
309 if ( dz > 0 && dz > distance )
319 if ( outputGroundBlock )
320 outputGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
321 if ( outputNonGroundBlock )
322 outputNonGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, val );
326 if ( outputGroundBlock )
327 outputGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, val );
328 if ( outputNonGroundBlock )
329 outputNonGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
334 if ( groundDestProvider )
335 groundDestProvider->writeBlock( outputGroundBlock.get(), mBand, tileLeft, tileTop );
336 if ( nonGroundDestProvider )
337 nonGroundDestProvider->writeBlock( outputNonGroundBlock.get(), mBand, tileLeft, tileTop );
339 if ( groundDestProvider )
340 groundDestProvider->setEditable(
false );
341 if ( nonGroundDestProvider )
342 nonGroundDestProvider->setEditable(
false );
345 outputs.insert( QStringLiteral(
"OUTPUT_GROUND" ), groundOutputFile );
346 outputs.insert( QStringLiteral(
"OUTPUT_NONGROUND" ), nonGroundOutputFile );
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.
@ Double
Double/float values.
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.