Branch data Line data Source code
1 : : /***************************************************************************
2 : : qgsrasterstackposition.cpp
3 : : ---------------------
4 : : begin : July 2020
5 : : copyright : (C) 2020 by Clemens Raffler
6 : : email : clemens dot raffler 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 : :
18 : : #include "qgsalgorithmrasterstackposition.h"
19 : : #include "qgsrasterprojector.h"
20 : : #include "qgsrasterfilewriter.h"
21 : : #include "qgsrasteranalysisutils.h"
22 : :
23 : : ///@cond PRIVATE
24 : :
25 : : //
26 : : //QgsRasterFrequencyByComparisonOperatorBase
27 : : //
28 : :
29 : 0 : QString QgsRasterStackPositionAlgorithmBase::group() const
30 : : {
31 : 0 : return QObject::tr( "Raster analysis" );
32 : : }
33 : :
34 : 0 : QString QgsRasterStackPositionAlgorithmBase::groupId() const
35 : : {
36 : 0 : return QStringLiteral( "rasteranalysis" );
37 : : }
38 : :
39 : 0 : void QgsRasterStackPositionAlgorithmBase::initAlgorithm( const QVariantMap & )
40 : : {
41 : 0 : addParameter( new QgsProcessingParameterMultipleLayers( QStringLiteral( "INPUT_RASTERS" ),
42 : 0 : QObject::tr( "Input raster layers" ), QgsProcessing::TypeRaster ) );
43 : :
44 : 0 : addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "REFERENCE_LAYER" ), QObject::tr( "Reference layer" ) ) );
45 : :
46 : 0 : addParameter( new QgsProcessingParameterBoolean( QStringLiteral( "IGNORE_NODATA" ), QObject::tr( "Ignore NoData values" ), false ) );
47 : :
48 : 0 : std::unique_ptr< QgsProcessingParameterNumber > output_nodata_parameter = std::make_unique< QgsProcessingParameterNumber >( QStringLiteral( "OUTPUT_NODATA_VALUE" ), QObject::tr( "Output NoData value" ), QgsProcessingParameterNumber::Double, -9999, true );
49 : 0 : output_nodata_parameter->setFlags( output_nodata_parameter->flags() | QgsProcessingParameterDefinition::FlagAdvanced );
50 : 0 : addParameter( output_nodata_parameter.release() );
51 : :
52 : 0 : addParameter( new QgsProcessingParameterRasterDestination( QStringLiteral( "OUTPUT" ),
53 : 0 : QObject::tr( "Output layer" ) ) );
54 : 0 : addOutput( new QgsProcessingOutputString( QStringLiteral( "EXTENT" ), QObject::tr( "Extent" ) ) );
55 : 0 : addOutput( new QgsProcessingOutputString( QStringLiteral( "CRS_AUTHID" ), QObject::tr( "CRS authority identifier" ) ) );
56 : 0 : addOutput( new QgsProcessingOutputNumber( QStringLiteral( "WIDTH_IN_PIXELS" ), QObject::tr( "Width in pixels" ) ) );
57 : 0 : addOutput( new QgsProcessingOutputNumber( QStringLiteral( "HEIGHT_IN_PIXELS" ), QObject::tr( "Height in pixels" ) ) );
58 : 0 : addOutput( new QgsProcessingOutputNumber( QStringLiteral( "TOTAL_PIXEL_COUNT" ), QObject::tr( "Total pixel count" ) ) );
59 : 0 : }
60 : :
61 : 0 : bool QgsRasterStackPositionAlgorithmBase::prepareAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
62 : : {
63 : 0 : QgsRasterLayer *referenceLayer = parameterAsRasterLayer( parameters, QStringLiteral( "REFERENCE_LAYER" ), context );
64 : 0 : if ( !referenceLayer )
65 : 0 : throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "REFERENCE_LAYER" ) ) );
66 : :
67 : 0 : mIgnoreNoData = parameterAsBool( parameters, QStringLiteral( "IGNORE_NODATA" ), context );
68 : 0 : mNoDataValue = parameterAsDouble( parameters, QStringLiteral( "OUTPUT_NODATA_VALUE" ), context );
69 : :
70 : 0 : mCrs = referenceLayer->crs();
71 : 0 : mRasterUnitsPerPixelX = referenceLayer->rasterUnitsPerPixelX();
72 : 0 : mRasterUnitsPerPixelY = referenceLayer->rasterUnitsPerPixelY();
73 : 0 : mLayerWidth = referenceLayer->width();
74 : 0 : mLayerHeight = referenceLayer->height();
75 : 0 : mExtent = referenceLayer->extent();
76 : :
77 : 0 : const QList< QgsMapLayer * > layers = parameterAsLayerList( parameters, QStringLiteral( "INPUT_RASTERS" ), context );
78 : 0 : QList< QgsRasterLayer * > rasterLayers;
79 : 0 : rasterLayers.reserve( layers.count() );
80 : 0 : for ( QgsMapLayer *l : layers )
81 : : {
82 : 0 : if ( feedback->isCanceled() )
83 : 0 : break; //in case some slow data sources are loaded
84 : :
85 : 0 : if ( l->type() == QgsMapLayerType::RasterLayer )
86 : : {
87 : 0 : QgsRasterLayer *layer = qobject_cast< QgsRasterLayer * >( l );
88 : 0 : QgsRasterAnalysisUtils::RasterLogicInput input;
89 : 0 : const int band = 1; //could be made dynamic
90 : 0 : input.hasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( band );
91 : 0 : input.sourceDataProvider.reset( layer->dataProvider()->clone() );
92 : 0 : input.interface = input.sourceDataProvider.get();
93 : : // add projector if necessary
94 : 0 : if ( layer->crs() != mCrs )
95 : : {
96 : 0 : input.projector = std::make_unique< QgsRasterProjector >();
97 : 0 : input.projector->setInput( input.sourceDataProvider.get() );
98 : 0 : input.projector->setCrs( layer->crs(), mCrs, context.transformContext() );
99 : 0 : input.interface = input.projector.get();
100 : 0 : }
101 : 0 : mInputs.emplace_back( std::move( input ) );
102 : 0 : }
103 : : }
104 : :
105 : : return true;
106 : 0 : }
107 : :
108 : 0 : QVariantMap QgsRasterStackPositionAlgorithmBase::processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
109 : : {
110 : 0 : const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
111 : 0 : QFileInfo fi( outputFile );
112 : 0 : const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
113 : :
114 : 0 : std::unique_ptr< QgsRasterFileWriter > writer = std::make_unique< QgsRasterFileWriter >( outputFile );
115 : 0 : writer->setOutputProviderKey( QStringLiteral( "gdal" ) );
116 : 0 : writer->setOutputFormat( outputFormat );
117 : 0 : std::unique_ptr<QgsRasterDataProvider > provider( writer->createOneBandRaster( Qgis::Int32, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
118 : 0 : if ( !provider )
119 : 0 : throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( outputFile ) );
120 : 0 : if ( !provider->isValid() )
121 : 0 : throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( outputFile, provider->error().message( QgsErrorMessage::Text ) ) );
122 : :
123 : 0 : provider->setNoDataValue( 1, mNoDataValue );
124 : 0 : qgssize layerSize = static_cast< qgssize >( mLayerWidth ) * static_cast< qgssize >( mLayerHeight );
125 : :
126 : 0 : int maxWidth = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_WIDTH;
127 : 0 : int maxHeight = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_HEIGHT;
128 : 0 : int nbBlocksWidth = static_cast< int>( std::ceil( 1.0 * mLayerWidth / maxWidth ) );
129 : 0 : int nbBlocksHeight = static_cast< int >( std::ceil( 1.0 * mLayerHeight / maxHeight ) );
130 : 0 : int nbBlocks = nbBlocksWidth * nbBlocksHeight;
131 : 0 : provider->setEditable( true );
132 : :
133 : 0 : QgsRasterIterator iter( provider.get() );
134 : 0 : iter.startRasterRead( 1, mLayerWidth, mLayerHeight, mExtent );
135 : 0 : int iterLeft = 0;
136 : 0 : int iterTop = 0;
137 : 0 : int iterCols = 0;
138 : 0 : int iterRows = 0;
139 : 0 : QgsRectangle blockExtent;
140 : :
141 : 0 : std::unique_ptr< QgsRasterBlock > outputBlock;
142 : 0 : while ( iter.readNextRasterPart( 1, iterCols, iterRows, outputBlock, iterLeft, iterTop, &blockExtent ) )
143 : : {
144 : 0 : std::vector< std::unique_ptr< QgsRasterBlock > > inputBlocks;
145 : 0 : for ( const QgsRasterAnalysisUtils::RasterLogicInput &i : mInputs )
146 : : {
147 : 0 : if ( feedback->isCanceled() )
148 : 0 : break; //in case some slow data sources are loaded
149 : 0 : for ( int band : i.bands )
150 : : {
151 : 0 : if ( feedback->isCanceled() )
152 : 0 : break; //in case some slow data sources are loaded
153 : 0 : std::unique_ptr< QgsRasterBlock > b( i.interface->block( band, blockExtent, iterCols, iterRows ) );
154 : 0 : inputBlocks.emplace_back( std::move( b ) );
155 : 0 : }
156 : : }
157 : :
158 : 0 : feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
159 : 0 : for ( int row = 0; row < iterRows; row++ )
160 : : {
161 : 0 : if ( feedback->isCanceled() )
162 : 0 : break;
163 : :
164 : 0 : for ( int col = 0; col < iterCols; col++ )
165 : : {
166 : 0 : bool noDataInStack = false;
167 : :
168 : 0 : if ( !inputBlocks.empty() )
169 : : {
170 : 0 : int position = findPosition( inputBlocks, row, col, noDataInStack );
171 : :
172 : 0 : if ( position == -1 || ( noDataInStack && !mIgnoreNoData ) )
173 : : {
174 : : //output cell will always be NoData if NoData occurs the current raster cell
175 : : //of the input blocks and NoData is not ignored
176 : : //this saves unnecessary iterations on the cellValueStack
177 : 0 : outputBlock->setValue( row, col, mNoDataValue );
178 : 0 : }
179 : : else
180 : : {
181 : 0 : outputBlock->setValue( row, col, position );
182 : : }
183 : 0 : }
184 : : else
185 : : {
186 : 0 : outputBlock->setValue( row, col, mNoDataValue );
187 : : }
188 : 0 : }
189 : 0 : }
190 : 0 : provider->writeBlock( outputBlock.get(), 1, iterLeft, iterTop );
191 : 0 : }
192 : 0 : provider->setEditable( false );
193 : :
194 : 0 : QVariantMap outputs;
195 : 0 : outputs.insert( QStringLiteral( "EXTENT" ), mExtent.toString() );
196 : 0 : outputs.insert( QStringLiteral( "CRS_AUTHID" ), mCrs.authid() );
197 : 0 : outputs.insert( QStringLiteral( "WIDTH_IN_PIXELS" ), mLayerWidth );
198 : 0 : outputs.insert( QStringLiteral( "HEIGHT_IN_PIXELS" ), mLayerHeight );
199 : 0 : outputs.insert( QStringLiteral( "TOTAL_PIXEL_COUNT" ), layerSize );
200 : 0 : outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );
201 : :
202 : 0 : return outputs;
203 : 0 : }
204 : :
205 : : //
206 : : // QgsRasterStackLowestPositionAlgorithm
207 : : //
208 : 0 : QString QgsRasterStackLowestPositionAlgorithm::displayName() const
209 : : {
210 : 0 : return QObject::tr( "Lowest position in raster stack" );
211 : : }
212 : :
213 : 0 : QString QgsRasterStackLowestPositionAlgorithm::name() const
214 : : {
215 : 0 : return QStringLiteral( "lowestpositioninrasterstack" );
216 : : }
217 : :
218 : 0 : QStringList QgsRasterStackLowestPositionAlgorithm::tags() const
219 : : {
220 : 0 : return QObject::tr( "cell,lowest,position,pixel,stack" ).split( ',' );
221 : 0 : }
222 : :
223 : 0 : QString QgsRasterStackLowestPositionAlgorithm::shortHelpString() const
224 : : {
225 : 0 : return QObject::tr( "The lowest position algorithm evaluates on a cell-by-cell basis the position "
226 : : "of the raster with the lowest value in a stack of rasters. Position counts start "
227 : : "with 1 and range to the total number of input rasters. The order of the input "
228 : : "rasters is relevant for the algorithm. If multiple rasters feature the lowest value, "
229 : : "the first raster will be used for the position value.\n "
230 : : "If multiband rasters are used in the data raster stack, the algorithm will always "
231 : : "perform the analysis on the first band of the rasters - use GDAL to use other bands in the analysis. "
232 : : "Any NoData cells in the raster layer stack will result in a NoData cell "
233 : : "in the output raster unless the \"ignore NoData\" parameter is checked. "
234 : : "The output NoData value can be set manually. The output rasters extent and resolution "
235 : : "is defined by a reference raster layer and is always of int32 type." );
236 : : }
237 : :
238 : 0 : QgsRasterStackLowestPositionAlgorithm *QgsRasterStackLowestPositionAlgorithm::createInstance() const
239 : : {
240 : 0 : return new QgsRasterStackLowestPositionAlgorithm();
241 : 0 : }
242 : :
243 : 0 : int QgsRasterStackLowestPositionAlgorithm::findPosition( std::vector< std::unique_ptr<QgsRasterBlock> > &inputBlocks, int &row, int &col, bool &noDataInRasterBlockStack )
244 : : {
245 : 0 : int lowestPosition = 0;
246 : :
247 : : //auxiliary variables
248 : 0 : int inputBlocksCount = inputBlocks.size();
249 : 0 : int currentPosition = 0;
250 : 0 : int noDataCount = 0;
251 : 0 : double firstValue = mNoDataValue;
252 : 0 : bool firstValueIsNoData = true;
253 : :
254 : 0 : while ( firstValueIsNoData && ( currentPosition < inputBlocksCount ) )
255 : : {
256 : : //check if all blocks are nodata/invalid
257 : 0 : std::unique_ptr<QgsRasterBlock> &firstBlock = inputBlocks.at( currentPosition );
258 : 0 : firstValue = firstBlock->valueAndNoData( row, col, firstValueIsNoData );
259 : :
260 : 0 : if ( !firstBlock->isValid() || firstValueIsNoData )
261 : : {
262 : 0 : noDataInRasterBlockStack = true;
263 : 0 : noDataCount++;
264 : 0 : }
265 : : else
266 : : {
267 : 0 : lowestPosition = currentPosition;
268 : : }
269 : 0 : currentPosition++;
270 : : }
271 : :
272 : 0 : if ( noDataCount == inputBlocksCount )
273 : : {
274 : 0 : noDataInRasterBlockStack = true;
275 : 0 : return -1; //all blocks are NoData
276 : : }
277 : : else
278 : : {
279 : : //scan for the lowest value
280 : 0 : while ( currentPosition < inputBlocksCount )
281 : : {
282 : 0 : std::unique_ptr< QgsRasterBlock > ¤tBlock = inputBlocks.at( currentPosition );
283 : :
284 : 0 : bool currentValueIsNoData = false;
285 : 0 : double currentValue = currentBlock->valueAndNoData( row, col, currentValueIsNoData );
286 : :
287 : 0 : if ( !currentBlock->isValid() || currentValueIsNoData )
288 : : {
289 : 0 : noDataInRasterBlockStack = true;
290 : 0 : noDataCount++;
291 : 0 : }
292 : : else
293 : : {
294 : 0 : if ( currentValue < firstValue )
295 : : {
296 : 0 : firstValue = currentValue;
297 : 0 : lowestPosition = currentPosition;
298 : 0 : }
299 : : }
300 : 0 : currentPosition++;
301 : : }
302 : : }
303 : : //the ArcGIS implementation uses 1 for first position value instead of 0 as in standard c++
304 : 0 : return ++lowestPosition; //therefore ++
305 : 0 : }
306 : :
307 : : //
308 : : // QgsRasterStackHighestPositionAlgorithmAlgorithm
309 : : //
310 : :
311 : 0 : QString QgsRasterStackHighestPositionAlgorithm::displayName() const
312 : : {
313 : 0 : return QObject::tr( "Highest position in raster stack" );
314 : : }
315 : :
316 : 0 : QString QgsRasterStackHighestPositionAlgorithm::name() const
317 : : {
318 : 0 : return QStringLiteral( "highestpositioninrasterstack" );
319 : : }
320 : :
321 : 0 : QStringList QgsRasterStackHighestPositionAlgorithm::tags() const
322 : : {
323 : 0 : return QObject::tr( "cell,highest,position,pixel,stack" ).split( ',' );
324 : 0 : }
325 : :
326 : 0 : QString QgsRasterStackHighestPositionAlgorithm::shortHelpString() const
327 : : {
328 : 0 : return QObject::tr( "The highest position algorithm evaluates on a cell-by-cell basis the position "
329 : : "of the raster with the highest value in a stack of rasters. Position counts start "
330 : : "with 1 and range to the total number of input rasters. The order of the input "
331 : : "rasters is relevant for the algorithm. If multiple rasters feature the highest value, "
332 : : "the first raster will be used for the position value.\n "
333 : : "If multiband rasters are used in the data raster stack, the algorithm will always "
334 : : "perform the analysis on the first band of the rasters - use GDAL to use other bands in the analysis. "
335 : : "Any NoData cells in the raster layer stack will result in a NoData cell "
336 : : "in the output raster unless the \"ignore NoData\" parameter is checked. "
337 : : "The output NoData value can be set manually. The output rasters extent and resolution "
338 : : "is defined by a reference raster layer and is always of int32 type." );
339 : : }
340 : :
341 : 0 : QgsRasterStackHighestPositionAlgorithm *QgsRasterStackHighestPositionAlgorithm::createInstance() const
342 : : {
343 : 0 : return new QgsRasterStackHighestPositionAlgorithm();
344 : 0 : }
345 : :
346 : 0 : int QgsRasterStackHighestPositionAlgorithm::findPosition( std::vector< std::unique_ptr< QgsRasterBlock> > &inputBlocks, int &row, int &col, bool &noDataInRasterBlockStack )
347 : : {
348 : 0 : int highestPosition = 0;
349 : :
350 : : //auxiliary variables
351 : 0 : int inputBlocksCount = inputBlocks.size();
352 : 0 : int currentPosition = 0;
353 : 0 : int noDataCount = 0;
354 : 0 : double firstValue = mNoDataValue;
355 : 0 : bool firstValueIsNoData = true;
356 : :
357 : 0 : while ( firstValueIsNoData && ( currentPosition < inputBlocksCount ) )
358 : : {
359 : : //check if all blocks are nodata/invalid
360 : 0 : std::unique_ptr<QgsRasterBlock> &firstBlock = inputBlocks.at( currentPosition );
361 : 0 : firstValue = firstBlock->valueAndNoData( row, col, firstValueIsNoData );
362 : :
363 : 0 : if ( !firstBlock->isValid() || firstValueIsNoData )
364 : : {
365 : 0 : noDataInRasterBlockStack = true;
366 : 0 : noDataCount++;
367 : 0 : }
368 : : else
369 : : {
370 : 0 : highestPosition = currentPosition;
371 : : }
372 : :
373 : 0 : currentPosition++;
374 : : }
375 : :
376 : 0 : if ( noDataCount == inputBlocksCount )
377 : : {
378 : 0 : noDataInRasterBlockStack = true;
379 : 0 : return -1; //all blocks are NoData
380 : : }
381 : : else
382 : : {
383 : : //scan for the lowest value
384 : 0 : while ( currentPosition < inputBlocksCount )
385 : : {
386 : 0 : std::unique_ptr< QgsRasterBlock > ¤tBlock = inputBlocks.at( currentPosition );
387 : :
388 : 0 : bool currentValueIsNoData = false;
389 : 0 : double currentValue = currentBlock->valueAndNoData( row, col, currentValueIsNoData );
390 : :
391 : 0 : if ( !currentBlock->isValid() || currentValueIsNoData )
392 : : {
393 : 0 : noDataInRasterBlockStack = true;
394 : 0 : noDataCount++;
395 : 0 : }
396 : : else
397 : : {
398 : 0 : if ( currentValue > firstValue )
399 : : {
400 : 0 : firstValue = currentValue;
401 : 0 : highestPosition = currentPosition;
402 : 0 : }
403 : : }
404 : 0 : currentPosition++;
405 : : }
406 : : }
407 : : //the ArcGIS implementation uses 1 for first position value instead of 0 as in standard c++
408 : 0 : return ++highestPosition; //therefore ++
409 : 0 : }
410 : :
411 : : ///@endcond
412 : :
|
