Branch data Line data Source code
1 : : /***************************************************************************
2 : : qgsalgorithmrasterfrequencybycomparisonoperator.cpp
3 : : ---------------------
4 : : begin : June 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 "qgsalgorithmrasterfrequencybycomparisonoperator.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 QgsRasterFrequencyByComparisonOperatorBase::group() const
30 : : {
31 : 0 : return QObject::tr( "Raster analysis" );
32 : : }
33 : :
34 : 0 : QString QgsRasterFrequencyByComparisonOperatorBase::groupId() const
35 : : {
36 : 0 : return QStringLiteral( "rasteranalysis" );
37 : : }
38 : :
39 : 0 : void QgsRasterFrequencyByComparisonOperatorBase::initAlgorithm( const QVariantMap & )
40 : : {
41 : 0 : addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT_VALUE_RASTER" ), QObject::tr( "Input value raster" ) ) );
42 : 0 : addParameter( new QgsProcessingParameterBand( QStringLiteral( "INPUT_VALUE_RASTER_BAND" ), QObject::tr( "Value raster band" ), 1, QStringLiteral( "INPUT_VALUE_RASTER" ) ) );
43 : :
44 : :
45 : 0 : addParameter( new QgsProcessingParameterMultipleLayers( QStringLiteral( "INPUT_RASTERS" ),
46 : 0 : QObject::tr( "Input raster layers" ), QgsProcessing::TypeRaster ) );
47 : :
48 : 0 : addParameter( new QgsProcessingParameterBoolean( QStringLiteral( "IGNORE_NODATA" ), QObject::tr( "Ignore NoData values" ), false ) );
49 : :
50 : 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 );
51 : 0 : output_nodata_parameter->setFlags( output_nodata_parameter->flags() | QgsProcessingParameterDefinition::FlagAdvanced );
52 : 0 : addParameter( output_nodata_parameter.release() );
53 : :
54 : 0 : addParameter( new QgsProcessingParameterRasterDestination( QStringLiteral( "OUTPUT" ),
55 : 0 : QObject::tr( "Output layer" ) ) );
56 : 0 : addOutput( new QgsProcessingOutputNumber( QStringLiteral( "OCCURRENCE_COUNT" ), QObject::tr( "Count of value occurrences" ) ) );
57 : 0 : addOutput( new QgsProcessingOutputNumber( QStringLiteral( "FOUND_LOCATIONS_COUNT" ), QObject::tr( "Count of cells with equal value occurrences" ) ) );
58 : 0 : addOutput( new QgsProcessingOutputNumber( QStringLiteral( "MEAN_FREQUENCY_PER_LOCATION" ), QObject::tr( "Mean frequency at valid cell locations" ) ) );
59 : 0 : addOutput( new QgsProcessingOutputString( QStringLiteral( "EXTENT" ), QObject::tr( "Extent" ) ) );
60 : 0 : addOutput( new QgsProcessingOutputString( QStringLiteral( "CRS_AUTHID" ), QObject::tr( "CRS authority identifier" ) ) );
61 : 0 : addOutput( new QgsProcessingOutputNumber( QStringLiteral( "WIDTH_IN_PIXELS" ), QObject::tr( "Width in pixels" ) ) );
62 : 0 : addOutput( new QgsProcessingOutputNumber( QStringLiteral( "HEIGHT_IN_PIXELS" ), QObject::tr( "Height in pixels" ) ) );
63 : 0 : addOutput( new QgsProcessingOutputNumber( QStringLiteral( "TOTAL_PIXEL_COUNT" ), QObject::tr( "Total pixel count" ) ) );
64 : 0 : }
65 : :
66 : 0 : bool QgsRasterFrequencyByComparisonOperatorBase::prepareAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
67 : : {
68 : 0 : QgsRasterLayer *inputValueRaster = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT_VALUE_RASTER" ), context );
69 : 0 : if ( !inputValueRaster )
70 : 0 : throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT_VALUE_RASTER" ) ) );
71 : :
72 : 0 : mInputValueRasterBand = parameterAsInt( parameters, QStringLiteral( "INPUT_VALUE_RASTER_BAND" ), context );
73 : 0 : mIgnoreNoData = parameterAsBool( parameters, QStringLiteral( "IGNORE_NODATA" ), context );
74 : :
75 : 0 : mInputValueRasterInterface.reset( inputValueRaster->dataProvider()->clone() );
76 : 0 : mNoDataValue = parameterAsDouble( parameters, QStringLiteral( "OUTPUT_NODATA_VALUE" ), context );
77 : 0 : mCrs = inputValueRaster->crs();
78 : 0 : mRasterUnitsPerPixelX = inputValueRaster->rasterUnitsPerPixelX();
79 : 0 : mRasterUnitsPerPixelY = inputValueRaster->rasterUnitsPerPixelY();
80 : 0 : mLayerWidth = inputValueRaster->width();
81 : 0 : mLayerHeight = inputValueRaster->height();
82 : 0 : mExtent = inputValueRaster->extent();
83 : :
84 : 0 : const QList< QgsMapLayer * > layers = parameterAsLayerList( parameters, QStringLiteral( "INPUT_RASTERS" ), context );
85 : 0 : QList< QgsRasterLayer * > rasterLayers;
86 : 0 : rasterLayers.reserve( layers.count() );
87 : 0 : for ( QgsMapLayer *l : layers )
88 : : {
89 : 0 : if ( feedback->isCanceled() )
90 : 0 : break; //in case some slow data sources are loaded
91 : :
92 : 0 : if ( l->type() == QgsMapLayerType::RasterLayer )
93 : : {
94 : 0 : QgsRasterLayer *layer = qobject_cast< QgsRasterLayer * >( l );
95 : 0 : QgsRasterAnalysisUtils::RasterLogicInput input;
96 : 0 : const int band = 1; //could be made dynamic
97 : 0 : input.hasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( band );
98 : 0 : input.sourceDataProvider.reset( layer->dataProvider()->clone() );
99 : 0 : input.interface = input.sourceDataProvider.get();
100 : : // add projector if necessary
101 : 0 : if ( layer->crs() != mCrs )
102 : : {
103 : 0 : input.projector = std::make_unique< QgsRasterProjector >();
104 : 0 : input.projector->setInput( input.sourceDataProvider.get() );
105 : 0 : input.projector->setCrs( layer->crs(), mCrs, context.transformContext() );
106 : 0 : input.interface = input.projector.get();
107 : 0 : }
108 : 0 : mInputs.emplace_back( std::move( input ) );
109 : 0 : }
110 : : }
111 : :
112 : : return true;
113 : 0 : }
114 : :
115 : 0 : QVariantMap QgsRasterFrequencyByComparisonOperatorBase::processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
116 : : {
117 : 0 : const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
118 : 0 : QFileInfo fi( outputFile );
119 : 0 : const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
120 : :
121 : 0 : std::unique_ptr< QgsRasterFileWriter > writer = std::make_unique< QgsRasterFileWriter >( outputFile );
122 : 0 : writer->setOutputProviderKey( QStringLiteral( "gdal" ) );
123 : 0 : writer->setOutputFormat( outputFormat );
124 : 0 : std::unique_ptr<QgsRasterDataProvider > provider( writer->createOneBandRaster( Qgis::Int32, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
125 : 0 : if ( !provider )
126 : 0 : throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( outputFile ) );
127 : 0 : if ( !provider->isValid() )
128 : 0 : throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( outputFile, provider->error().message( QgsErrorMessage::Text ) ) );
129 : :
130 : 0 : provider->setNoDataValue( 1, mNoDataValue );
131 : 0 : qgssize layerSize = static_cast< qgssize >( mLayerWidth ) * static_cast< qgssize >( mLayerHeight );
132 : :
133 : 0 : int maxWidth = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_WIDTH;
134 : 0 : int maxHeight = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_HEIGHT;
135 : 0 : int nbBlocksWidth = static_cast< int>( std::ceil( 1.0 * mLayerWidth / maxWidth ) );
136 : 0 : int nbBlocksHeight = static_cast< int >( std::ceil( 1.0 * mLayerHeight / maxHeight ) );
137 : 0 : int nbBlocks = nbBlocksWidth * nbBlocksHeight;
138 : 0 : provider->setEditable( true );
139 : :
140 : 0 : QgsRasterIterator iter( mInputValueRasterInterface.get() );
141 : 0 : iter.startRasterRead( mInputValueRasterBand, mLayerWidth, mLayerHeight, mExtent );
142 : 0 : int iterLeft = 0;
143 : 0 : int iterTop = 0;
144 : 0 : int iterCols = 0;
145 : 0 : int iterRows = 0;
146 : 0 : QgsRectangle blockExtent;
147 : :
148 : 0 : unsigned long long occurrenceCount = 0;
149 : 0 : unsigned long long noDataLocationsCount = 0;
150 : 0 : std::unique_ptr< QgsRasterBlock > inputBlock;
151 : 0 : while ( iter.readNextRasterPart( 1, iterCols, iterRows, inputBlock, iterLeft, iterTop, &blockExtent ) )
152 : : {
153 : 0 : std::vector< std::unique_ptr< QgsRasterBlock > > inputBlocks;
154 : 0 : for ( const QgsRasterAnalysisUtils::RasterLogicInput &i : mInputs )
155 : : {
156 : 0 : if ( feedback->isCanceled() )
157 : 0 : break; //in case some slow data sources are loaded
158 : 0 : for ( int band : i.bands )
159 : : {
160 : 0 : if ( feedback->isCanceled() )
161 : 0 : break; //in case some slow data sources are loaded
162 : 0 : std::unique_ptr< QgsRasterBlock > b( i.interface->block( band, blockExtent, iterCols, iterRows ) );
163 : 0 : inputBlocks.emplace_back( std::move( b ) );
164 : 0 : }
165 : : }
166 : :
167 : 0 : std::unique_ptr< QgsRasterBlock > outputBlock = std::make_unique<QgsRasterBlock>( Qgis::Int32, iterCols, iterRows );
168 : 0 : feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
169 : 0 : for ( int row = 0; row < iterRows; row++ )
170 : : {
171 : 0 : if ( feedback->isCanceled() )
172 : 0 : break;
173 : :
174 : 0 : for ( int col = 0; col < iterCols; col++ )
175 : : {
176 : 0 : bool valueRasterCellIsNoData = false;
177 : 0 : double value = inputBlock->valueAndNoData( row, col, valueRasterCellIsNoData );
178 : :
179 : 0 : if ( valueRasterCellIsNoData && !mIgnoreNoData )
180 : : {
181 : : //output cell will always be NoData if NoData occurs in valueRaster or cellValueStack and NoData is not ignored
182 : : //this saves unnecessary iterations on the cellValueStack
183 : 0 : outputBlock->setValue( row, col, mNoDataValue );
184 : 0 : noDataLocationsCount++;
185 : 0 : }
186 : : else
187 : : {
188 : 0 : bool noDataInStack = false;
189 : 0 : std::vector<double> cellValues = QgsRasterAnalysisUtils::getCellValuesFromBlockStack( inputBlocks, row, col, noDataInStack );
190 : :
191 : 0 : if ( noDataInStack && !mIgnoreNoData )
192 : : {
193 : 0 : outputBlock->setValue( row, col, mNoDataValue );
194 : 0 : noDataLocationsCount++;
195 : 0 : }
196 : : else
197 : : {
198 : 0 : int frequency = applyComparisonOperator( value, cellValues );
199 : 0 : outputBlock->setValue( row, col, frequency );
200 : 0 : occurrenceCount += frequency;
201 : : }
202 : 0 : }
203 : 0 : }
204 : 0 : }
205 : 0 : provider->writeBlock( outputBlock.get(), 1, iterLeft, iterTop );
206 : 0 : }
207 : 0 : provider->setEditable( false );
208 : :
209 : 0 : unsigned long long foundLocationsCount = layerSize - noDataLocationsCount;
210 : 0 : double meanEqualCountPerValidLocation = static_cast<double>( occurrenceCount ) / static_cast<double>( foundLocationsCount * mInputs.size() );
211 : :
212 : 0 : QVariantMap outputs;
213 : 0 : outputs.insert( QStringLiteral( "OCCURRENCE_COUNT" ), occurrenceCount );
214 : 0 : outputs.insert( QStringLiteral( "FOUND_LOCATIONS_COUNT" ), foundLocationsCount );
215 : 0 : outputs.insert( QStringLiteral( "MEAN_FREQUENCY_PER_LOCATION" ), meanEqualCountPerValidLocation );
216 : 0 : outputs.insert( QStringLiteral( "EXTENT" ), mExtent.toString() );
217 : 0 : outputs.insert( QStringLiteral( "CRS_AUTHID" ), mCrs.authid() );
218 : 0 : outputs.insert( QStringLiteral( "WIDTH_IN_PIXELS" ), mLayerWidth );
219 : 0 : outputs.insert( QStringLiteral( "HEIGHT_IN_PIXELS" ), mLayerHeight );
220 : 0 : outputs.insert( QStringLiteral( "TOTAL_PIXEL_COUNT" ), layerSize );
221 : 0 : outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );
222 : :
223 : 0 : return outputs;
224 : 0 : }
225 : :
226 : : //
227 : : // QgsRasterFrequencyByEqualOperatorAlgorithm
228 : : //
229 : :
230 : 0 : QString QgsRasterFrequencyByEqualOperatorAlgorithm::displayName() const
231 : : {
232 : 0 : return QObject::tr( "Equal to frequency" );
233 : : }
234 : :
235 : 0 : QString QgsRasterFrequencyByEqualOperatorAlgorithm::name() const
236 : : {
237 : 0 : return QStringLiteral( "equaltofrequency" );
238 : : }
239 : :
240 : 0 : QStringList QgsRasterFrequencyByEqualOperatorAlgorithm::tags() const
241 : : {
242 : 0 : return QObject::tr( "cell,equal,frequency,pixel,stack" ).split( ',' );
243 : 0 : }
244 : :
245 : 0 : QString QgsRasterFrequencyByEqualOperatorAlgorithm::shortHelpString() const
246 : : {
247 : 0 : return QObject::tr( "The Equal to frequency algorithm evaluates on a cell-by-cell basis the frequency "
248 : : "(number of times) the values of an input stack of rasters are equal "
249 : : "to the value of a value raster. \n "
250 : : "If multiband rasters are used in the data raster stack, the algorithm will always "
251 : : "perform the analysis on the first band of the rasters - use GDAL to use other bands in the analysis. "
252 : : "The input value layer serves as reference layer for the sample layers. "
253 : : "Any NoData cells in the value raster or the data layer stack will result in a NoData cell "
254 : : "in the output raster if the ignore NoData parameter is not checked. "
255 : : "The output NoData value can be set manually. The output rasters extent and resolution "
256 : : "is defined by the input raster layer and is always of int32 type." );
257 : : }
258 : :
259 : 0 : QgsRasterFrequencyByEqualOperatorAlgorithm *QgsRasterFrequencyByEqualOperatorAlgorithm::createInstance() const
260 : : {
261 : 0 : return new QgsRasterFrequencyByEqualOperatorAlgorithm();
262 : 0 : }
263 : :
264 : 0 : int QgsRasterFrequencyByEqualOperatorAlgorithm::applyComparisonOperator( double searchValue, std::vector<double> cellValueStack )
265 : : {
266 : 0 : return static_cast<int>( std::count( cellValueStack.begin(), cellValueStack.end(), searchValue ) );
267 : : }
268 : :
269 : : //
270 : : // QgsRasterFrequencyByGreaterThanOperatorAlgorithm
271 : : //
272 : :
273 : 0 : QString QgsRasterFrequencyByGreaterThanOperatorAlgorithm::displayName() const
274 : : {
275 : 0 : return QObject::tr( "Greater than frequency" );
276 : : }
277 : :
278 : 0 : QString QgsRasterFrequencyByGreaterThanOperatorAlgorithm::name() const
279 : : {
280 : 0 : return QStringLiteral( "greaterthanfrequency" );
281 : : }
282 : :
283 : 0 : QStringList QgsRasterFrequencyByGreaterThanOperatorAlgorithm::tags() const
284 : : {
285 : 0 : return QObject::tr( "cell,greater,frequency,pixel,stack" ).split( ',' );
286 : 0 : }
287 : :
288 : 0 : QString QgsRasterFrequencyByGreaterThanOperatorAlgorithm::shortHelpString() const
289 : : {
290 : 0 : return QObject::tr( "The Greater than frequency algorithm evaluates on a cell-by-cell basis the frequency "
291 : : "(number of times) the values of an input stack of rasters are greater than "
292 : : "the value of a value raster. \n "
293 : : "If multiband rasters are used in the data raster stack, the algorithm will always "
294 : : "perform the analysis on the first band of the rasters - use GDAL to use other bands in the analysis. "
295 : : "The input value layer serves as reference layer for the sample layers. "
296 : : "Any NoData cells in the value raster or the data layer stack will result in a NoData cell "
297 : : "in the output raster if the ignore NoData parameter is not checked. "
298 : : "The output NoData value can be set manually. The output rasters extent and resolution "
299 : : "is defined by the input raster layer and is always of int32 type." );
300 : : }
301 : :
302 : 0 : QgsRasterFrequencyByGreaterThanOperatorAlgorithm *QgsRasterFrequencyByGreaterThanOperatorAlgorithm::createInstance() const
303 : : {
304 : 0 : return new QgsRasterFrequencyByGreaterThanOperatorAlgorithm();
305 : 0 : }
306 : :
307 : 0 : int QgsRasterFrequencyByGreaterThanOperatorAlgorithm::applyComparisonOperator( double searchValue, std::vector<double> cellValueStack )
308 : : {
309 : 0 : return static_cast<int>( std::count_if( cellValueStack.begin(), cellValueStack.end(), [&]( double const & stackValue ) { return stackValue > searchValue; } ) );
310 : : }
311 : :
312 : : //
313 : : // QgsRasterFrequencyByLessThanOperatorAlgorithm
314 : : //
315 : :
316 : 0 : QString QgsRasterFrequencyByLessThanOperatorAlgorithm::displayName() const
317 : : {
318 : 0 : return QObject::tr( "Less than frequency" );
319 : : }
320 : :
321 : 0 : QString QgsRasterFrequencyByLessThanOperatorAlgorithm::name() const
322 : : {
323 : 0 : return QStringLiteral( "lessthanfrequency" );
324 : : }
325 : :
326 : 0 : QStringList QgsRasterFrequencyByLessThanOperatorAlgorithm::tags() const
327 : : {
328 : 0 : return QObject::tr( "cell,less,lower,frequency,pixel,stack" ).split( ',' );
329 : 0 : }
330 : :
331 : 0 : QString QgsRasterFrequencyByLessThanOperatorAlgorithm::shortHelpString() const
332 : : {
333 : 0 : return QObject::tr( "The Less than frequency algorithm evaluates on a cell-by-cell basis the frequency "
334 : : "(number of times) the values of an input stack of rasters are less than "
335 : : "the value of a value raster. \n "
336 : : "If multiband rasters are used in the data raster stack, the algorithm will always "
337 : : "perform the analysis on the first band of the rasters - use GDAL to use other bands in the analysis. "
338 : : "The input value layer serves as reference layer for the sample layers. "
339 : : "Any NoData cells in the value raster or the data layer stack will result in a NoData cell "
340 : : "in the output raster if the ignore NoData parameter is not checked. "
341 : : "The output NoData value can be set manually. The output rasters extent and resolution "
342 : : "is defined by the input raster layer and is always of int32 type." );
343 : : }
344 : :
345 : 0 : QgsRasterFrequencyByLessThanOperatorAlgorithm *QgsRasterFrequencyByLessThanOperatorAlgorithm::createInstance() const
346 : : {
347 : 0 : return new QgsRasterFrequencyByLessThanOperatorAlgorithm();
348 : 0 : }
349 : :
350 : 0 : int QgsRasterFrequencyByLessThanOperatorAlgorithm::applyComparisonOperator( double searchValue, std::vector<double> cellValueStack )
351 : : {
352 : 0 : return static_cast<int>( std::count_if( cellValueStack.begin(), cellValueStack.end(), [&]( double const & stackValue ) { return stackValue < searchValue; } ) );
353 : : }
354 : :
355 : : ///@endcond
356 : :
357 : :
|
