QGIS API Documentation 3.39.0-Master (67e056379ed)
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< 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" ) }} ) }} ) );
96 createOptsParam->setFlags( createOptsParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
97 addParameter( createOptsParam.release() );
98
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() );
103
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() );
108}
109
110QgsRasterDtmSlopeBasedFilterAlgorithm *QgsRasterDtmSlopeBasedFilterAlgorithm::createInstance() const
111{
112 return new QgsRasterDtmSlopeBasedFilterAlgorithm();
113}
114
115bool QgsRasterDtmSlopeBasedFilterAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
116{
117 QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT" ), context );
118 if ( !layer )
119 throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT" ) ) );
120
121 const int band = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
122
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 )
126 .arg( layer->bandCount() ) );
127
128 mInterface.reset( layer->dataProvider()->clone() );
129 mHasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( band );
130 mLayerWidth = layer->width();
131 mLayerHeight = layer->height();
132 mExtent = layer->extent();
133 mCrs = layer->crs();
134 mRasterUnitsPerPixelX = layer->rasterUnitsPerPixelX();
135 mRasterUnitsPerPixelY = layer->rasterUnitsPerPixelY();
136 mDataType = layer->dataProvider()->dataType( mBand );
137 mNoData = layer->dataProvider()->sourceNoDataValue( mBand );
138 return true;
139}
140
141QVariantMap QgsRasterDtmSlopeBasedFilterAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
142{
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;
147
148 if ( !groundOutputFile.isEmpty() )
149 {
150 const QFileInfo fi( groundOutputFile );
151 const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
152
153 groundWriter = std::make_unique< QgsRasterFileWriter >( groundOutputFile );
154 groundWriter->setOutputProviderKey( QStringLiteral( "gdal" ) );
155 if ( !createOptions.isEmpty() )
156 {
157 groundWriter->setCreateOptions( createOptions.split( '|' ) );
158 }
159 groundWriter->setOutputFormat( outputFormat );
160
161 groundDestProvider.reset( groundWriter->createOneBandRaster( mDataType, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
162
163 if ( !groundDestProvider )
164 throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( groundOutputFile ) );
165 if ( !groundDestProvider->isValid() )
166 throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( groundOutputFile, groundDestProvider->error().message( QgsErrorMessage::Text ) ) );
167
168 groundDestProvider->setNoDataValue( 1, mNoData );
169 groundDestProvider->setEditable( true );
170 }
171
172 const QString nonGroundOutputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT_NONGROUND" ), context );
173 std::unique_ptr< QgsRasterFileWriter > nonGroundWriter;
174 std::unique_ptr<QgsRasterDataProvider > nonGroundDestProvider;
175
176 if ( !nonGroundOutputFile.isEmpty() )
177 {
178 const QFileInfo fi( groundOutputFile );
179 const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
180
181 nonGroundWriter = std::make_unique< QgsRasterFileWriter >( nonGroundOutputFile );
182 nonGroundWriter->setOutputProviderKey( QStringLiteral( "gdal" ) );
183 if ( !createOptions.isEmpty() )
184 {
185 nonGroundWriter->setCreateOptions( createOptions.split( '|' ) );
186 }
187 nonGroundWriter->setOutputFormat( outputFormat );
188
189 nonGroundDestProvider.reset( nonGroundWriter->createOneBandRaster( mDataType, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
190
191 if ( !nonGroundDestProvider )
192 throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( nonGroundOutputFile ) );
193 if ( !nonGroundDestProvider->isValid() )
194 throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( nonGroundOutputFile, nonGroundDestProvider->error().message( QgsErrorMessage::Text ) ) );
195
196 nonGroundDestProvider->setNoDataValue( 1, mNoData );
197 nonGroundDestProvider->setEditable( true );
198 }
199
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;
205
206 const int radius = parameterAsInt( parameters, QStringLiteral( "RADIUS" ), context );
207
208 const double terrainSlopePercent = parameterAsDouble( parameters, QStringLiteral( "TERRAIN_SLOPE" ), context ) / 100; //20.0 / 100 * 0.143;
209 const int filterModification = parameterAsEnum( parameters, QStringLiteral( "FILTER_MODIFICATION" ), context );
210 const double standardDeviation = parameterAsDouble( parameters, QStringLiteral( "STANDARD_DEVIATION" ), context );
211
212 // create kernel
213 QVector<double> kernel;
214 kernel.reserve( ( radius * 2 ) * ( radius * 2 ) );
215 int kernelSize = 0;
216 for ( int y = -radius; y <= radius; y++ )
217 {
218 for ( int x = -radius; x <= radius; x++ )
219 {
220 const double distance = std::sqrt( x * x + y * y );
221 if ( distance < radius )
222 {
223 kernelSize++;
224 kernel.push_back( x );
225 kernel.push_back( y );
226 switch ( filterModification )
227 {
228 case 0:
229 kernel.push_back( distance * terrainSlopePercent );
230 break;
231
232 case 1:
233 kernel.push_back( distance * terrainSlopePercent + 1.65 * std::sqrt( 2 * standardDeviation ) );
234 break;
235
236 case 2:
237 {
238 const double dz = distance * terrainSlopePercent - 1.65 * std::sqrt( 2 * standardDeviation );
239 kernel.push_back( dz > 0 ? dz : 0 );
240 break;
241 }
242 }
243 }
244 }
245 }
246
247 QgsRasterIterator iter( mInterface.get(), radius );
248 iter.startRasterRead( 1, mLayerWidth, mLayerHeight, mExtent );
249 int iterLeft = 0;
250 int iterTop = 0;
251 int iterCols = 0;
252 int iterRows = 0;
253 int tileLeft = 0;
254 int tileTop = 0;
255 int tileCols = 0;
256 int tileRows = 0;
257
258 QgsRectangle blockExtent;
259
260 std::unique_ptr< QgsRasterBlock > inputBlock;
261 int blockIndex = 0;
262 while ( iter.readNextRasterPart( 1, iterCols, iterRows, inputBlock, iterLeft, iterTop, &blockExtent, &tileCols, &tileRows, &tileLeft, &tileTop ) )
263 {
264 std::unique_ptr< QgsRasterBlock > outputGroundBlock;
265 if ( groundDestProvider )
266 outputGroundBlock = std::make_unique< QgsRasterBlock >( mDataType, tileCols, tileRows );
267
268 std::unique_ptr< QgsRasterBlock > outputNonGroundBlock;
269 if ( nonGroundDestProvider )
270 outputNonGroundBlock = std::make_unique< QgsRasterBlock >( mDataType, tileCols, tileRows );
271
272 double baseProgress = static_cast< double >( blockIndex ) / numBlocks;
273 feedback->setProgress( 100.0 * baseProgress );
274 blockIndex++;
275 if ( feedback->isCanceled() )
276 break;
277
278 const int tileBoundaryLeft = tileLeft - iterLeft;
279 const int tileBoundaryTop = tileTop - iterTop;
280
281 const double rowProgressStep = 1.0 / numBlocks / tileRows;
282 double rowProgress = 0;
283 for ( int row = tileBoundaryTop; row < tileBoundaryTop + tileRows; row++ )
284 {
285 if ( feedback->isCanceled() )
286 break;
287
288 feedback->setProgress( 100.0 * ( baseProgress + rowProgress ) );
289 rowProgress += rowProgressStep;
290
291 for ( int col = tileBoundaryLeft; col < tileBoundaryLeft + tileCols; col++ )
292 {
293 if ( feedback->isCanceled() )
294 break;
295
296 bool isNoData = false;
297 const double val = inputBlock->valueAndNoData( row, col, isNoData );
298 if ( isNoData )
299 {
300 if ( outputGroundBlock )
301 outputGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
302 if ( outputNonGroundBlock )
303 outputNonGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
304 }
305 else
306 {
307 bool nonGround = false;
308 const double *kernelData = kernel.constData();
309 for ( int i = 0; i < kernelSize; ++i )
310 {
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 ) ) )
317 {
318 bool otherIsNoData = false;
319 const double otherVal = inputBlock->valueAndNoData( rRow, rCol, otherIsNoData );
320 if ( !otherIsNoData )
321 {
322 const double dz = val - otherVal;
323 if ( dz > 0 && dz > distance )
324 {
325 nonGround = true;
326 break;
327 }
328 }
329 }
330 }
331 if ( nonGround )
332 {
333 if ( outputGroundBlock )
334 outputGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
335 if ( outputNonGroundBlock )
336 outputNonGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, val );
337 }
338 else
339 {
340 if ( outputGroundBlock )
341 outputGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, val );
342 if ( outputNonGroundBlock )
343 outputNonGroundBlock->setValue( row - tileBoundaryTop, col - tileBoundaryLeft, mNoData );
344 }
345 }
346 }
347 }
348 if ( groundDestProvider )
349 groundDestProvider->writeBlock( outputGroundBlock.get(), mBand, tileLeft, tileTop );
350 if ( nonGroundDestProvider )
351 nonGroundDestProvider->writeBlock( outputNonGroundBlock.get(), mBand, tileLeft, tileTop );
352 }
353 if ( groundDestProvider )
354 groundDestProvider->setEditable( false );
355 if ( nonGroundDestProvider )
356 nonGroundDestProvider->setEditable( false );
357
358 QVariantMap outputs;
359 outputs.insert( QStringLiteral( "OUTPUT_GROUND" ), groundOutputFile );
360 outputs.insert( QStringLiteral( "OUTPUT_NONGROUND" ), nonGroundOutputFile );
361 return outputs;
362}
363
364
366
367
368
@ 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: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.