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< QgsProcessingParameterString > createOptsParam = std::make_unique< QgsProcessingParameterString >( QStringLiteral(
"CREATE_OPTIONS" ), QObject::tr(
"Creation options" ), QVariant(),
false,
true );
95 createOptsParam->setMetadata( QVariantMap( {{QStringLiteral(
"widget_wrapper" ), QVariantMap( {{QStringLiteral(
"widget_type" ), QStringLiteral(
"rasteroptions" ) }} ) }} ) );
97 addParameter( createOptsParam.release() );
99 std::unique_ptr< QgsProcessingParameterRasterDestination > outputLayerGroundParam = std::make_unique<QgsProcessingParameterRasterDestination>( QStringLiteral(
"OUTPUT_GROUND" ),
100 QObject::tr(
"Output layer (ground)" ), QVariant(),
true,
true );
101 outputLayerGroundParam->setHelp( QObject::tr(
"The filtered DEM containing only cells classified as ground." ) );
102 addParameter( outputLayerGroundParam.release() );
104 std::unique_ptr< QgsProcessingParameterRasterDestination > outputLayerNonGroundParam = std::make_unique<QgsProcessingParameterRasterDestination>( QStringLiteral(
"OUTPUT_NONGROUND" ),
105 QObject::tr(
"Output layer (non-ground objects)" ), QVariant(),
true,
false );
106 outputLayerNonGroundParam->setHelp( QObject::tr(
"The non-ground objects removed by the filter." ) );
107 addParameter( outputLayerNonGroundParam.release() );
110QgsRasterDtmSlopeBasedFilterAlgorithm *QgsRasterDtmSlopeBasedFilterAlgorithm::createInstance()
const
112 return new QgsRasterDtmSlopeBasedFilterAlgorithm();
117 QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral(
"INPUT" ), context );
121 const int band = parameterAsInt( parameters, QStringLiteral(
"BAND" ), context );
123 mBand = parameterAsInt( parameters, QStringLiteral(
"BAND" ), context );
124 if ( mBand < 1 || mBand > layer->
bandCount() )
125 throw QgsProcessingException( QObject::tr(
"Invalid band number for BAND (%1): Valid values for input raster are 1 to %2" ).arg( mBand )
130 mLayerWidth = layer->
width();
131 mLayerHeight = layer->
height();
132 mExtent = layer->
extent();
143 const QString createOptions = parameterAsString( parameters, QStringLiteral(
"CREATE_OPTIONS" ), context ).trimmed();
144 const QString groundOutputFile = parameterAsOutputLayer( parameters, QStringLiteral(
"OUTPUT_GROUND" ), context );
145 std::unique_ptr< QgsRasterFileWriter > groundWriter;
146 std::unique_ptr<QgsRasterDataProvider > groundDestProvider;
148 if ( !groundOutputFile.isEmpty() )
150 const QFileInfo fi( groundOutputFile );
153 groundWriter = std::make_unique< QgsRasterFileWriter >( groundOutputFile );
154 groundWriter->setOutputProviderKey( QStringLiteral(
"gdal" ) );
155 if ( !createOptions.isEmpty() )
157 groundWriter->setCreateOptions( createOptions.split(
'|' ) );
159 groundWriter->setOutputFormat( outputFormat );
161 groundDestProvider.reset( groundWriter->createOneBandRaster( mDataType, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
163 if ( !groundDestProvider )
164 throw QgsProcessingException( QObject::tr(
"Could not create raster output: %1" ).arg( groundOutputFile ) );
165 if ( !groundDestProvider->isValid() )
168 groundDestProvider->setNoDataValue( 1, mNoData );
169 groundDestProvider->setEditable(
true );
172 const QString nonGroundOutputFile = parameterAsOutputLayer( parameters, QStringLiteral(
"OUTPUT_NONGROUND" ), context );
173 std::unique_ptr< QgsRasterFileWriter > nonGroundWriter;
174 std::unique_ptr<QgsRasterDataProvider > nonGroundDestProvider;
176 if ( !nonGroundOutputFile.isEmpty() )
178 const QFileInfo fi( groundOutputFile );
181 nonGroundWriter = std::make_unique< QgsRasterFileWriter >( nonGroundOutputFile );
182 nonGroundWriter->setOutputProviderKey( QStringLiteral(
"gdal" ) );
183 if ( !createOptions.isEmpty() )
185 nonGroundWriter->setCreateOptions( createOptions.split(
'|' ) );
187 nonGroundWriter->setOutputFormat( outputFormat );
189 nonGroundDestProvider.reset( nonGroundWriter->createOneBandRaster( mDataType, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
191 if ( !nonGroundDestProvider )
192 throw QgsProcessingException( QObject::tr(
"Could not create raster output: %1" ).arg( nonGroundOutputFile ) );
193 if ( !nonGroundDestProvider->isValid() )
196 nonGroundDestProvider->setNoDataValue( 1, mNoData );
197 nonGroundDestProvider->setEditable(
true );
202 const int numBlocksX =
static_cast< int>( std::ceil( 1.0 * mLayerWidth / blockWidth ) );
203 const int numBlocksY =
static_cast< int >( std::ceil( 1.0 * mLayerHeight / blockHeight ) );
204 const int numBlocks = numBlocksX * numBlocksY;
206 const int radius = parameterAsInt( parameters, QStringLiteral(
"RADIUS" ), context );
208 const double terrainSlopePercent = parameterAsDouble( parameters, QStringLiteral(
"TERRAIN_SLOPE" ), context ) / 100;
209 const int filterModification = parameterAsEnum( parameters, QStringLiteral(
"FILTER_MODIFICATION" ), context );
210 const double standardDeviation = parameterAsDouble( parameters, QStringLiteral(
"STANDARD_DEVIATION" ), context );
213 QVector<double> kernel;
214 kernel.reserve( ( radius * 2 ) * ( radius * 2 ) );
216 for (
int y = -radius; y <= radius; y++ )
218 for (
int x = -radius; x <= radius; x++ )
220 const double distance = std::sqrt( x * x + y * y );
221 if ( distance < radius )
224 kernel.push_back( x );
225 kernel.push_back( y );
226 switch ( filterModification )
229 kernel.push_back( distance * terrainSlopePercent );
233 kernel.push_back( distance * terrainSlopePercent + 1.65 * std::sqrt( 2 * standardDeviation ) );
238 const double dz = distance * terrainSlopePercent - 1.65 * std::sqrt( 2 * standardDeviation );
239 kernel.push_back( dz > 0 ? dz : 0 );
248 iter.startRasterRead( 1, mLayerWidth, mLayerHeight, mExtent );
260 std::unique_ptr< QgsRasterBlock > inputBlock;
262 while ( iter.readNextRasterPart( 1, iterCols, iterRows, inputBlock, iterLeft, iterTop, &blockExtent, &tileCols, &tileRows, &tileLeft, &tileTop ) )
264 std::unique_ptr< QgsRasterBlock > outputGroundBlock;
265 if ( groundDestProvider )
266 outputGroundBlock = std::make_unique< QgsRasterBlock >( mDataType, tileCols, tileRows );
268 std::unique_ptr< QgsRasterBlock > outputNonGroundBlock;
269 if ( nonGroundDestProvider )
270 outputNonGroundBlock = std::make_unique< QgsRasterBlock >( mDataType, tileCols, tileRows );
272 double baseProgress =
static_cast< double >( blockIndex ) / numBlocks;
278 const int tileBoundaryLeft = tileLeft - iterLeft;
279 const int tileBoundaryTop = tileTop - iterTop;
281 const double rowProgressStep = 1.0 / numBlocks / tileRows;
282 double rowProgress = 0;
283 for (
int row = tileBoundaryTop; row < tileBoundaryTop + tileRows; row++ )
288 feedback->
setProgress( 100.0 * ( baseProgress + rowProgress ) );
289 rowProgress += rowProgressStep;
291 for (
int col = tileBoundaryLeft; col < tileBoundaryLeft + tileCols; col++ )
296 bool isNoData =
false;
297 const double val = inputBlock->valueAndNoData( row, col, isNoData );
300 if ( outputGroundBlock )
301 outputGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
302 if ( outputNonGroundBlock )
303 outputNonGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
307 bool nonGround =
false;
308 const double *kernelData = kernel.constData();
309 for (
int i = 0; i < kernelSize; ++i )
311 const int dx =
static_cast< int >( *kernelData++ );
312 const int dy =
static_cast< int >( *kernelData++ );
313 const double distance = *kernelData++;
314 const int rCol = col + dx;
315 const int rRow = row + dy;
316 if ( rCol >= 0 && ( rCol < ( iterLeft + iterCols ) ) && rRow >= 0 && ( rRow < ( iterTop + iterRows ) ) )
318 bool otherIsNoData =
false;
319 const double otherVal = inputBlock->valueAndNoData( rRow, rCol, otherIsNoData );
320 if ( !otherIsNoData )
322 const double dz = val - otherVal;
323 if ( dz > 0 && dz > distance )
333 if ( outputGroundBlock )
334 outputGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
335 if ( outputNonGroundBlock )
336 outputNonGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, val );
340 if ( outputGroundBlock )
341 outputGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, val );
342 if ( outputNonGroundBlock )
343 outputNonGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
348 if ( groundDestProvider )
349 groundDestProvider->writeBlock( outputGroundBlock.get(), mBand, tileLeft, tileTop );
350 if ( nonGroundDestProvider )
351 nonGroundDestProvider->writeBlock( outputNonGroundBlock.get(), mBand, tileLeft, tileTop );
353 if ( groundDestProvider )
354 groundDestProvider->setEditable(
false );
355 if ( nonGroundDestProvider )
356 nonGroundDestProvider->setEditable(
false );
359 outputs.insert( QStringLiteral(
"OUTPUT_GROUND" ), groundOutputFile );
360 outputs.insert( QStringLiteral(
"OUTPUT_NONGROUND" ), nonGroundOutputFile );
@ 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.