Branch data Line data Source code
1 : : /***************************************************************************
2 : : qgsrasteranalysisutils.cpp
3 : : ---------------------
4 : : Date : June 2018
5 : : Copyright : (C) 2018 by Nyall Dawson
6 : : Email : nyall dot dawson at gmail dot com
7 : : ***************************************************************************
8 : : * *
9 : : * This program is free software; you can redistribute it and/or modify *
10 : : * it under the terms of the GNU General Public License as published by *
11 : : * the Free Software Foundation; either version 2 of the License, or *
12 : : * (at your option) any later version. *
13 : : * *
14 : : ***************************************************************************/
15 : :
16 : : #include "qgsrasteranalysisutils.h"
17 : :
18 : : #include "qgsfeedback.h"
19 : : #include "qgsrasterblock.h"
20 : : #include "qgsrasteriterator.h"
21 : : #include "qgsgeos.h"
22 : : #include "qgsprocessingparameters.h"
23 : : #include <map>
24 : : #include <unordered_map>
25 : : #include <unordered_set>
26 : : #include <cmath>
27 : : ///@cond PRIVATE
28 : :
29 : 0 : void QgsRasterAnalysisUtils::cellInfoForBBox( const QgsRectangle &rasterBBox, const QgsRectangle &featureBBox, double cellSizeX, double cellSizeY,
30 : : int &nCellsX, int &nCellsY, int rasterWidth, int rasterHeight, QgsRectangle &rasterBlockExtent )
31 : : {
32 : : //get intersecting bbox
33 : 0 : QgsRectangle intersectBox = rasterBBox.intersect( featureBBox );
34 : 0 : if ( intersectBox.isEmpty() )
35 : : {
36 : 0 : nCellsX = 0;
37 : 0 : nCellsY = 0;
38 : 0 : rasterBlockExtent = QgsRectangle();
39 : 0 : return;
40 : : }
41 : :
42 : : //get offset in pixels in x- and y- direction
43 : 0 : int offsetX = static_cast< int >( std::floor( ( intersectBox.xMinimum() - rasterBBox.xMinimum() ) / cellSizeX ) );
44 : 0 : int offsetY = static_cast< int >( std::floor( ( rasterBBox.yMaximum() - intersectBox.yMaximum() ) / cellSizeY ) );
45 : :
46 : 0 : int maxColumn = static_cast< int >( std::floor( ( intersectBox.xMaximum() - rasterBBox.xMinimum() ) / cellSizeX ) ) + 1;
47 : 0 : int maxRow = static_cast< int >( std::floor( ( rasterBBox.yMaximum() - intersectBox.yMinimum() ) / cellSizeY ) ) + 1;
48 : :
49 : 0 : nCellsX = maxColumn - offsetX;
50 : 0 : nCellsY = maxRow - offsetY;
51 : :
52 : : //avoid access to cells outside of the raster (may occur because of rounding)
53 : 0 : nCellsX = std::min( offsetX + nCellsX, rasterWidth ) - offsetX;
54 : 0 : nCellsY = std::min( offsetY + nCellsY, rasterHeight ) - offsetY;
55 : :
56 : 0 : rasterBlockExtent = QgsRectangle( rasterBBox.xMinimum() + offsetX * cellSizeX,
57 : 0 : rasterBBox.yMaximum() - offsetY * cellSizeY,
58 : 0 : rasterBBox.xMinimum() + ( nCellsX + offsetX ) * cellSizeX,
59 : 0 : rasterBBox.yMaximum() - ( nCellsY + offsetY ) * cellSizeY );
60 : 0 : }
61 : :
62 : 0 : void QgsRasterAnalysisUtils::statisticsFromMiddlePointTest( QgsRasterInterface *rasterInterface, int rasterBand, const QgsGeometry &poly, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle &rasterBBox, const std::function<void( double )> &addValue, bool skipNodata )
63 : : {
64 : 0 : std::unique_ptr< QgsGeometryEngine > polyEngine( QgsGeometry::createGeometryEngine( poly.constGet( ) ) );
65 : 0 : if ( !polyEngine )
66 : : {
67 : 0 : return;
68 : : }
69 : 0 : polyEngine->prepareGeometry();
70 : :
71 : 0 : QgsRasterIterator iter( rasterInterface );
72 : 0 : iter.startRasterRead( rasterBand, nCellsX, nCellsY, rasterBBox );
73 : :
74 : 0 : std::unique_ptr< QgsRasterBlock > block;
75 : 0 : int iterLeft = 0;
76 : 0 : int iterTop = 0;
77 : 0 : int iterCols = 0;
78 : 0 : int iterRows = 0;
79 : 0 : QgsRectangle blockExtent;
80 : 0 : bool isNoData = false;
81 : 0 : while ( iter.readNextRasterPart( rasterBand, iterCols, iterRows, block, iterLeft, iterTop, &blockExtent ) )
82 : : {
83 : 0 : double cellCenterY = blockExtent.yMaximum() - 0.5 * cellSizeY;
84 : :
85 : 0 : for ( int row = 0; row < iterRows; ++row )
86 : : {
87 : 0 : double cellCenterX = blockExtent.xMinimum() + 0.5 * cellSizeX;
88 : 0 : for ( int col = 0; col < iterCols; ++col )
89 : : {
90 : 0 : const double pixelValue = block->valueAndNoData( row, col, isNoData );
91 : 0 : if ( validPixel( pixelValue ) && ( !skipNodata || !isNoData ) )
92 : : {
93 : 0 : QgsPoint cellCenter( cellCenterX, cellCenterY );
94 : 0 : if ( polyEngine->contains( &cellCenter ) )
95 : : {
96 : 0 : addValue( pixelValue );
97 : 0 : }
98 : 0 : }
99 : 0 : cellCenterX += cellSizeX;
100 : 0 : }
101 : 0 : cellCenterY -= cellSizeY;
102 : 0 : }
103 : : }
104 : 0 : }
105 : :
106 : 0 : void QgsRasterAnalysisUtils::statisticsFromPreciseIntersection( QgsRasterInterface *rasterInterface, int rasterBand, const QgsGeometry &poly, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle &rasterBBox, const std::function<void( double, double )> &addValue, bool skipNodata )
107 : : {
108 : 0 : QgsGeometry pixelRectGeometry;
109 : :
110 : 0 : double hCellSizeX = cellSizeX / 2.0;
111 : 0 : double hCellSizeY = cellSizeY / 2.0;
112 : 0 : double pixelArea = cellSizeX * cellSizeY;
113 : 0 : double weight = 0;
114 : :
115 : 0 : std::unique_ptr< QgsGeometryEngine > polyEngine( QgsGeometry::createGeometryEngine( poly.constGet( ) ) );
116 : 0 : if ( !polyEngine )
117 : : {
118 : 0 : return;
119 : : }
120 : 0 : polyEngine->prepareGeometry();
121 : :
122 : 0 : QgsRasterIterator iter( rasterInterface );
123 : 0 : iter.startRasterRead( rasterBand, nCellsX, nCellsY, rasterBBox );
124 : :
125 : 0 : std::unique_ptr< QgsRasterBlock > block;
126 : 0 : int iterLeft = 0;
127 : 0 : int iterTop = 0;
128 : 0 : int iterCols = 0;
129 : 0 : int iterRows = 0;
130 : 0 : QgsRectangle blockExtent;
131 : 0 : bool isNoData = false;
132 : 0 : while ( iter.readNextRasterPart( rasterBand, iterCols, iterRows, block, iterLeft, iterTop, &blockExtent ) )
133 : : {
134 : 0 : double currentY = blockExtent.yMaximum() - 0.5 * cellSizeY;
135 : 0 : for ( int row = 0; row < iterRows; ++row )
136 : : {
137 : 0 : double currentX = blockExtent.xMinimum() + 0.5 * cellSizeX;
138 : 0 : for ( int col = 0; col < iterCols; ++col )
139 : : {
140 : 0 : const double pixelValue = block->valueAndNoData( row, col, isNoData );
141 : 0 : if ( validPixel( pixelValue ) && ( !skipNodata || !isNoData ) )
142 : : {
143 : 0 : pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) );
144 : : // GEOS intersects tests on prepared geometry is MAGNITUDES faster than calculating the intersection itself,
145 : : // so we first test to see if there IS an intersection before doing the actual calculation
146 : 0 : if ( !pixelRectGeometry.isNull() && polyEngine->intersects( pixelRectGeometry.constGet() ) )
147 : : {
148 : : //intersection
149 : 0 : QgsGeometry intersectGeometry = pixelRectGeometry.intersection( poly );
150 : 0 : if ( !intersectGeometry.isEmpty() )
151 : : {
152 : 0 : double intersectionArea = intersectGeometry.area();
153 : 0 : if ( intersectionArea > 0.0 )
154 : : {
155 : 0 : weight = intersectionArea / pixelArea;
156 : 0 : addValue( pixelValue, weight );
157 : 0 : }
158 : 0 : }
159 : 0 : }
160 : 0 : }
161 : 0 : currentX += cellSizeX;
162 : 0 : }
163 : 0 : currentY -= cellSizeY;
164 : 0 : }
165 : : }
166 : 0 : }
167 : :
168 : 0 : bool QgsRasterAnalysisUtils::validPixel( double value )
169 : : {
170 : 0 : return !std::isnan( value );
171 : : }
172 : :
173 : 0 : void QgsRasterAnalysisUtils::mapToPixel( const double x, const double y, const QgsRectangle bounds, const double unitsPerPixelX, const double unitsPerPixelY, int &px, int &py )
174 : : {
175 : 0 : px = trunc( ( x - bounds.xMinimum() ) / unitsPerPixelX );
176 : 0 : py = trunc( ( y - bounds.yMaximum() ) / -unitsPerPixelY );
177 : 0 : }
178 : :
179 : 0 : void QgsRasterAnalysisUtils::pixelToMap( const int px, const int py, const QgsRectangle bounds, const double unitsPerPixelX, const double unitsPerPixelY, double &x, double &y )
180 : : {
181 : 0 : x = bounds.xMinimum() + ( px + 0.5 ) * unitsPerPixelX;
182 : 0 : y = bounds.yMaximum() - ( py + 0.5 ) * unitsPerPixelY;
183 : 0 : }
184 : :
185 : 2 : static QVector< QPair< QString, Qgis::DataType > > sDataTypes;
186 : :
187 : 0 : void populateDataTypes()
188 : : {
189 : 0 : if ( sDataTypes.empty() )
190 : : {
191 : 0 : sDataTypes.append( qMakePair( QStringLiteral( "Byte" ), Qgis::Byte ) );
192 : 0 : sDataTypes.append( qMakePair( QStringLiteral( "Int16" ), Qgis::Int16 ) );
193 : 0 : sDataTypes.append( qMakePair( QStringLiteral( "UInt16" ), Qgis::UInt16 ) );
194 : 0 : sDataTypes.append( qMakePair( QStringLiteral( "Int32" ), Qgis::Int32 ) );
195 : 0 : sDataTypes.append( qMakePair( QStringLiteral( "UInt32" ), Qgis::UInt32 ) );
196 : 0 : sDataTypes.append( qMakePair( QStringLiteral( "Float32" ), Qgis::Float32 ) );
197 : 0 : sDataTypes.append( qMakePair( QStringLiteral( "Float64" ), Qgis::Float64 ) );
198 : 0 : sDataTypes.append( qMakePair( QStringLiteral( "CInt16" ), Qgis::CInt16 ) );
199 : 0 : sDataTypes.append( qMakePair( QStringLiteral( "CInt32" ), Qgis::CInt32 ) );
200 : 0 : sDataTypes.append( qMakePair( QStringLiteral( "CFloat32" ), Qgis::CFloat32 ) );
201 : 0 : sDataTypes.append( qMakePair( QStringLiteral( "CFloat64" ), Qgis::CFloat64 ) );
202 : 0 : }
203 : 0 : }
204 : :
205 : 0 : std::unique_ptr<QgsProcessingParameterDefinition> QgsRasterAnalysisUtils::createRasterTypeParameter( const QString &name, const QString &description, Qgis::DataType defaultType )
206 : : {
207 : 0 : populateDataTypes();
208 : :
209 : 0 : QStringList names;
210 : 0 : int defaultChoice = 0;
211 : 0 : int i = 0;
212 : 0 : for ( auto it = sDataTypes.constBegin(); it != sDataTypes.constEnd(); ++it )
213 : : {
214 : 0 : names.append( it->first );
215 : 0 : if ( it->second == defaultType )
216 : 0 : defaultChoice = i;
217 : 0 : i++;
218 : 0 : }
219 : :
220 : 0 : return std::make_unique< QgsProcessingParameterEnum >( name, description, names, false, defaultChoice );
221 : 0 : }
222 : :
223 : 0 : Qgis::DataType QgsRasterAnalysisUtils::rasterTypeChoiceToDataType( int choice )
224 : : {
225 : 0 : if ( choice < 0 || choice >= sDataTypes.count() )
226 : 0 : return Qgis::Float32;
227 : :
228 : 0 : return sDataTypes.value( choice ).second;
229 : 0 : }
230 : :
231 : 0 : void QgsRasterAnalysisUtils::applyRasterLogicOperator( const std::vector< QgsRasterAnalysisUtils::RasterLogicInput > &inputs, QgsRasterDataProvider *destinationRaster, double outputNoDataValue, const bool treatNoDataAsFalse,
232 : : int width, int height, const QgsRectangle &extent, QgsFeedback *feedback,
233 : : std::function<void( const std::vector< std::unique_ptr< QgsRasterBlock > > &, bool &, bool &, int, int, bool )> &applyLogicFunc,
234 : : qgssize &noDataCount, qgssize &trueCount, qgssize &falseCount )
235 : : {
236 : 0 : int maxWidth = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_WIDTH;
237 : 0 : int maxHeight = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_HEIGHT;
238 : 0 : int nbBlocksWidth = static_cast< int>( std::ceil( 1.0 * width / maxWidth ) );
239 : 0 : int nbBlocksHeight = static_cast< int >( std::ceil( 1.0 * height / maxHeight ) );
240 : 0 : int nbBlocks = nbBlocksWidth * nbBlocksHeight;
241 : :
242 : 0 : destinationRaster->setEditable( true );
243 : 0 : QgsRasterIterator outputIter( destinationRaster );
244 : 0 : outputIter.startRasterRead( 1, width, height, extent );
245 : :
246 : 0 : int iterLeft = 0;
247 : 0 : int iterTop = 0;
248 : 0 : int iterCols = 0;
249 : 0 : int iterRows = 0;
250 : 0 : QgsRectangle blockExtent;
251 : 0 : std::unique_ptr< QgsRasterBlock > outputBlock;
252 : 0 : while ( outputIter.readNextRasterPart( 1, iterCols, iterRows, outputBlock, iterLeft, iterTop, &blockExtent ) )
253 : : {
254 : 0 : std::vector< std::unique_ptr< QgsRasterBlock > > inputBlocks;
255 : 0 : for ( const QgsRasterAnalysisUtils::RasterLogicInput &i : inputs )
256 : : {
257 : 0 : for ( int band : i.bands )
258 : : {
259 : 0 : std::unique_ptr< QgsRasterBlock > b( i.interface->block( band, blockExtent, iterCols, iterRows ) );
260 : 0 : inputBlocks.emplace_back( std::move( b ) );
261 : 0 : }
262 : : }
263 : :
264 : 0 : feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
265 : 0 : for ( int row = 0; row < iterRows; row++ )
266 : : {
267 : 0 : if ( feedback->isCanceled() )
268 : 0 : break;
269 : :
270 : 0 : for ( int column = 0; column < iterCols; column++ )
271 : : {
272 : 0 : bool res = false;
273 : 0 : bool resIsNoData = false;
274 : 0 : applyLogicFunc( inputBlocks, res, resIsNoData, row, column, treatNoDataAsFalse );
275 : 0 : if ( resIsNoData )
276 : 0 : noDataCount++;
277 : 0 : else if ( res )
278 : 0 : trueCount++;
279 : : else
280 : 0 : falseCount++;
281 : :
282 : 0 : outputBlock->setValue( row, column, resIsNoData ? outputNoDataValue : ( res ? 1 : 0 ) );
283 : 0 : }
284 : 0 : }
285 : 0 : destinationRaster->writeBlock( outputBlock.get(), 1, iterLeft, iterTop );
286 : 0 : }
287 : 0 : destinationRaster->setEditable( false );
288 : 0 : }
289 : :
290 : 0 : std::vector<double> QgsRasterAnalysisUtils::getCellValuesFromBlockStack( const std::vector< std::unique_ptr< QgsRasterBlock > > &inputBlocks, int &row, int &col, bool &noDataInStack )
291 : : {
292 : : //get all values from inputBlocks
293 : 0 : std::vector<double> cellValues;
294 : 0 : bool hasNoData = false;
295 : 0 : cellValues.reserve( inputBlocks.size() );
296 : :
297 : 0 : for ( auto &block : inputBlocks )
298 : : {
299 : 0 : double value = 0;
300 : 0 : if ( !block || !block->isValid() )
301 : : {
302 : 0 : noDataInStack = true;
303 : 0 : break;
304 : : }
305 : : else
306 : : {
307 : 0 : value = block->valueAndNoData( row, col, hasNoData );
308 : 0 : if ( hasNoData )
309 : : {
310 : 0 : noDataInStack = true;
311 : 0 : continue; //NoData is not included in the cell value vector
312 : : }
313 : : else
314 : : {
315 : 0 : cellValues.push_back( value );
316 : : }
317 : : }
318 : : }
319 : 0 : return cellValues;
320 : 0 : }
321 : :
322 : 0 : double QgsRasterAnalysisUtils::meanFromCellValues( std::vector<double> &cellValues, int stackSize )
323 : : {
324 : 0 : double sum = std::accumulate( cellValues.begin(), cellValues.end(), 0.0 );
325 : 0 : double mean = sum / static_cast<double>( stackSize );
326 : 0 : return mean;
327 : : }
328 : :
329 : 0 : double QgsRasterAnalysisUtils::medianFromCellValues( std::vector<double> &cellValues, int stackSize )
330 : : {
331 : 0 : std::sort( cellValues.begin(), cellValues.end() );
332 : 0 : bool even = ( stackSize % 2 ) < 1;
333 : 0 : if ( even )
334 : : {
335 : 0 : return ( cellValues[stackSize / 2 - 1] + cellValues[stackSize / 2] ) / 2.0;
336 : : }
337 : : else //odd
338 : : {
339 : 0 : return cellValues[( stackSize + 1 ) / 2 - 1];
340 : : }
341 : 0 : }
342 : :
343 : :
344 : 0 : double QgsRasterAnalysisUtils::stddevFromCellValues( std::vector<double> &cellValues, int stackSize )
345 : : {
346 : 0 : double variance = varianceFromCellValues( cellValues, stackSize );
347 : 0 : double stddev = std::sqrt( variance );
348 : 0 : return stddev;
349 : : }
350 : :
351 : 0 : double QgsRasterAnalysisUtils::varianceFromCellValues( std::vector<double> &cellValues, int stackSize )
352 : : {
353 : 0 : double mean = meanFromCellValues( cellValues, stackSize );
354 : 0 : double accum = 0.0;
355 : 0 : for ( int i = 0; i < stackSize; i++ )
356 : : {
357 : 0 : accum += std::pow( ( cellValues.at( i ) - mean ), 2.0 );
358 : 0 : }
359 : 0 : double variance = accum / static_cast<double>( stackSize );
360 : 0 : return variance;
361 : : }
362 : :
363 : 0 : double QgsRasterAnalysisUtils::maximumFromCellValues( std::vector<double> &cellValues )
364 : : {
365 : 0 : return *std::max_element( cellValues.begin(), cellValues.end() );
366 : : }
367 : :
368 : 0 : double QgsRasterAnalysisUtils::minimumFromCellValues( std::vector<double> &cellValues )
369 : : {
370 : 0 : return *std::min_element( cellValues.begin(), cellValues.end() );
371 : : }
372 : :
373 : 0 : double QgsRasterAnalysisUtils::majorityFromCellValues( std::vector<double> &cellValues, const double noDataValue, int stackSize )
374 : : {
375 : 0 : if ( stackSize == 1 )
376 : : {
377 : : //output will be same as input if only one layer is entered
378 : 0 : return cellValues[0];
379 : : }
380 : 0 : else if ( stackSize == 2 )
381 : : {
382 : : //if only two layers are input, return NoData if values are not the same (eg. no Majority could be found)
383 : 0 : return ( qgsDoubleNear( cellValues[0], cellValues[1] ) ) ? cellValues[0] : noDataValue;
384 : : }
385 : 0 : else if ( std::adjacent_find( cellValues.begin(), cellValues.end(), std::not_equal_to<double>() ) == cellValues.end() )
386 : : {
387 : : //check if all values in cellValues are equal
388 : : //output will be same as input if all cellValues of the stack are the same
389 : 0 : return cellValues[0];
390 : : }
391 : : else
392 : : {
393 : : //search for majority using hash map [O(n)]
394 : 0 : std::unordered_map<double, int> map;
395 : :
396 : 0 : for ( int i = 0; i < stackSize; i++ )
397 : : {
398 : 0 : map[cellValues[i]]++;
399 : 0 : }
400 : :
401 : 0 : int maxCount = 0;
402 : 0 : bool multipleMajorities = false;
403 : 0 : double result = noDataValue;
404 : 0 : for ( auto pair : map )
405 : : {
406 : 0 : if ( maxCount < pair.second )
407 : : {
408 : 0 : result = pair.first;
409 : 0 : maxCount = pair.second;
410 : 0 : multipleMajorities = false;
411 : 0 : }
412 : 0 : else if ( maxCount == pair.second )
413 : : {
414 : 0 : multipleMajorities = true;
415 : 0 : }
416 : : }
417 : 0 : return multipleMajorities ? noDataValue : result;
418 : 0 : }
419 : 0 : }
420 : :
421 : 0 : double QgsRasterAnalysisUtils::minorityFromCellValues( std::vector<double> &cellValues, const double noDataValue, int stackSize )
422 : : {
423 : 0 : if ( stackSize == 1 )
424 : : {
425 : : //output will be same as input if only one layer is entered
426 : 0 : return cellValues[0];
427 : : }
428 : 0 : else if ( stackSize == 2 )
429 : : {
430 : : //if only two layers are input, return NoData if values are not the same (eg. no minority could be found)
431 : 0 : return ( qgsDoubleNear( cellValues[0], cellValues[1] ) ) ? cellValues[0] : noDataValue;
432 : : }
433 : 0 : else if ( std::adjacent_find( cellValues.begin(), cellValues.end(), std::not_equal_to<double>() ) == cellValues.end() )
434 : : {
435 : : //check if all values in cellValues are equal
436 : : //output will be same as input if all cellValues of the stack are the same
437 : 0 : return cellValues[0];
438 : : }
439 : : else
440 : : {
441 : : //search for minority using hash map [O(n)]
442 : 0 : std::unordered_map<double, int> map;
443 : :
444 : 0 : for ( int i = 0; i < stackSize; i++ )
445 : : {
446 : 0 : map[cellValues[i]]++;
447 : 0 : }
448 : :
449 : 0 : int minCount = stackSize;
450 : 0 : bool multipleMinorities = false;
451 : 0 : double result = noDataValue; //result will stay NoData if no minority value exists
452 : 0 : for ( auto pair : map )
453 : : {
454 : 0 : if ( minCount > pair.second )
455 : : {
456 : 0 : result = pair.first;
457 : 0 : minCount = pair.second;
458 : 0 : multipleMinorities = false;
459 : 0 : }
460 : 0 : else if ( minCount == pair.second )
461 : : {
462 : 0 : multipleMinorities = true;
463 : 0 : }
464 : : }
465 : 0 : return multipleMinorities ? noDataValue : result;
466 : 0 : }
467 : 0 : }
468 : :
469 : 0 : double QgsRasterAnalysisUtils::rangeFromCellValues( std::vector<double> &cellValues )
470 : : {
471 : 0 : double max = *std::max_element( cellValues.begin(), cellValues.end() );
472 : 0 : double min = *std::min_element( cellValues.begin(), cellValues.end() );
473 : 0 : return max - min;
474 : : }
475 : :
476 : 0 : double QgsRasterAnalysisUtils::varietyFromCellValues( std::vector<double> &cellValues )
477 : : {
478 : 0 : std::unordered_set<double> uniqueValues( cellValues.begin(), cellValues.end() );
479 : 0 : return uniqueValues.size();
480 : 0 : }
481 : :
482 : 0 : double QgsRasterAnalysisUtils::nearestRankPercentile( std::vector<double> &cellValues, int stackSize, double percentile )
483 : : {
484 : : //if percentile equals 0 -> pick the first element of the ordered list
485 : 0 : std::sort( cellValues.begin(), cellValues.end() );
486 : :
487 : 0 : int i = 0;
488 : 0 : if ( percentile > 0 )
489 : : {
490 : 0 : i = std::ceil( percentile * static_cast<double>( stackSize ) ) - 1;
491 : 0 : }
492 : :
493 : 0 : return cellValues[i];
494 : : }
495 : :
496 : 0 : double QgsRasterAnalysisUtils::interpolatedPercentileInc( std::vector<double> &cellValues, int stackSize, double percentile )
497 : : {
498 : 0 : std::sort( cellValues.begin(), cellValues.end() );
499 : 0 : double x = ( percentile * ( stackSize - 1 ) );
500 : :
501 : 0 : int i = static_cast<int>( std::floor( x ) );
502 : 0 : double xFraction = std::fmod( x, 1 );
503 : :
504 : 0 : if ( stackSize == 1 )
505 : : {
506 : 0 : return cellValues[0];
507 : : }
508 : 0 : else if ( stackSize == 2 )
509 : : {
510 : 0 : return cellValues[0] + ( cellValues[1] - cellValues[0] ) * percentile;
511 : : }
512 : : else
513 : : {
514 : 0 : return cellValues[i] + ( cellValues[i + 1] - cellValues[i] ) * xFraction;
515 : : }
516 : 0 : }
517 : :
518 : 0 : double QgsRasterAnalysisUtils::interpolatedPercentileExc( std::vector<double> &cellValues, int stackSize, double percentile, double noDataValue )
519 : : {
520 : 0 : std::sort( cellValues.begin(), cellValues.end() );
521 : 0 : double x = ( percentile * ( stackSize + 1 ) );
522 : :
523 : 0 : int i = static_cast<int>( std::floor( x ) ) - 1;
524 : 0 : double xFraction = std::fmod( x, 1 );
525 : 0 : double lowerExcValue = 1.0 / ( static_cast<double>( stackSize ) + 1.0 );
526 : 0 : double upperExcValue = static_cast<double>( stackSize ) / ( static_cast<double>( stackSize ) + 1.0 );
527 : :
528 : 0 : if ( stackSize < 2 || ( ( percentile < lowerExcValue || percentile > upperExcValue ) ) )
529 : : {
530 : 0 : return noDataValue;
531 : : }
532 : : else
533 : : {
534 : 0 : return cellValues[i] + ( cellValues[i + 1] - cellValues[i] ) * xFraction;
535 : : }
536 : 0 : }
537 : :
538 : 0 : double QgsRasterAnalysisUtils::interpolatedPercentRankInc( std::vector<double> &cellValues, int stackSize, double value, double noDataValue )
539 : : {
540 : 0 : std::sort( cellValues.begin(), cellValues.end() );
541 : :
542 : 0 : if ( value < cellValues[0] || value > cellValues[stackSize - 1] )
543 : : {
544 : 0 : return noDataValue;
545 : : }
546 : : else
547 : : {
548 : 0 : for ( int i = 0; i < stackSize - 1; i++ )
549 : : {
550 : 0 : if ( cellValues[i] <= value && cellValues[i + 1] >= value )
551 : : {
552 : 0 : double fraction = 0.0;
553 : :
554 : : //make sure that next number in the distribution is not the same to prevent NaN fractions
555 : 0 : if ( !qgsDoubleNear( cellValues[i], cellValues[i + 1] ) )
556 : 0 : fraction = ( value - cellValues[i] ) / ( cellValues[i + 1] - cellValues[i] );
557 : :
558 : 0 : return ( fraction + i ) / ( stackSize - 1 );
559 : : }
560 : 0 : }
561 : 0 : return noDataValue;
562 : : }
563 : 0 : }
564 : :
565 : 0 : double QgsRasterAnalysisUtils::interpolatedPercentRankExc( std::vector<double> &cellValues, int stackSize, double value, double noDataValue )
566 : : {
567 : 0 : std::sort( cellValues.begin(), cellValues.end() );
568 : :
569 : 0 : if ( value < cellValues[0] || value > cellValues[stackSize - 1] )
570 : : {
571 : 0 : return noDataValue;
572 : : }
573 : : else
574 : : {
575 : 0 : for ( int i = 0; i < stackSize - 1; i++ )
576 : : {
577 : 0 : if ( cellValues[i] <= value && cellValues[i + 1] >= value )
578 : : {
579 : 0 : double fraction = 0.0;
580 : :
581 : : //make sure that next number in the distribution is not the same to prevent NaN fractions
582 : 0 : if ( !qgsDoubleNear( cellValues[i], cellValues[i + 1] ) )
583 : 0 : fraction = ( value - cellValues[i] ) / ( cellValues[i + 1] - cellValues[i] );
584 : :
585 : 0 : return ( ( i + 1 ) + fraction ) / ( stackSize + 1 );
586 : : }
587 : 0 : }
588 : 0 : return noDataValue;
589 : : }
590 : 0 : }
591 : :
592 : :
593 : : ///@endcond PRIVATE
594 : :
|
