24#include <unordered_map>
25#include <unordered_set>
29void QgsRasterAnalysisUtils::cellInfoForBBox(
const QgsRectangle &rasterBBox,
const QgsRectangle &featureBBox,
double cellSizeX,
double cellSizeY,
30 int &nCellsX,
int &nCellsY,
int rasterWidth,
int rasterHeight,
QgsRectangle &rasterBlockExtent )
43 const int offsetX =
static_cast< int >( std::floor( ( intersectBox.
xMinimum() - rasterBBox.
xMinimum() ) / cellSizeX ) );
44 const int offsetY =
static_cast< int >( std::floor( ( rasterBBox.
yMaximum() - intersectBox.
yMaximum() ) / cellSizeY ) );
46 const int maxColumn =
static_cast< int >( std::floor( ( intersectBox.
xMaximum() - rasterBBox.
xMinimum() ) / cellSizeX ) ) + 1;
47 const int maxRow =
static_cast< int >( std::floor( ( rasterBBox.
yMaximum() - intersectBox.
yMinimum() ) / cellSizeY ) ) + 1;
49 nCellsX = maxColumn - offsetX;
50 nCellsY = maxRow - offsetY;
53 nCellsX = std::min( offsetX + nCellsX, rasterWidth ) - offsetX;
54 nCellsY = std::min( offsetY + nCellsY, rasterHeight ) - offsetY;
57 rasterBBox.
yMaximum() - offsetY * cellSizeY,
58 rasterBBox.
xMinimum() + ( nCellsX + offsetX ) * cellSizeX,
59 rasterBBox.
yMaximum() - ( nCellsY + offsetY ) * cellSizeY );
62void 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 )
64 std::unique_ptr< QgsGeos > polyEngine = std::make_unique< QgsGeos >( poly.
constGet( ) );
69 polyEngine->prepareGeometry();
72 iter.startRasterRead( rasterBand, nCellsX, nCellsY, rasterBBox );
74 std::unique_ptr< QgsRasterBlock > block;
80 bool isNoData =
false;
81 while ( iter.readNextRasterPart( rasterBand, iterCols, iterRows, block, iterLeft, iterTop, &blockExtent ) )
83 double cellCenterY = blockExtent.
yMaximum() - 0.5 * cellSizeY;
85 for (
int row = 0; row < iterRows; ++row )
87 double cellCenterX = blockExtent.
xMinimum() + 0.5 * cellSizeX;
88 for (
int col = 0; col < iterCols; ++col )
90 const double pixelValue = block->valueAndNoData( row, col, isNoData );
91 if ( validPixel( pixelValue ) && ( !skipNodata || !isNoData ) )
93 if ( polyEngine->contains( cellCenterX, cellCenterY ) )
95 addValue( pixelValue );
98 cellCenterX += cellSizeX;
100 cellCenterY -= cellSizeY;
105void 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 )
109 const double hCellSizeX = cellSizeX / 2.0;
110 const double hCellSizeY = cellSizeY / 2.0;
111 const double pixelArea = cellSizeX * cellSizeY;
119 polyEngine->prepareGeometry();
122 iter.startRasterRead( rasterBand, nCellsX, nCellsY, rasterBBox );
124 std::unique_ptr< QgsRasterBlock > block;
130 bool isNoData =
false;
131 while ( iter.readNextRasterPart( rasterBand, iterCols, iterRows, block, iterLeft, iterTop, &blockExtent ) )
133 double currentY = blockExtent.
yMaximum() - 0.5 * cellSizeY;
134 for (
int row = 0; row < iterRows; ++row )
136 double currentX = blockExtent.
xMinimum() + 0.5 * cellSizeX;
137 for (
int col = 0; col < iterCols; ++col )
139 const double pixelValue = block->valueAndNoData( row, col, isNoData );
140 if ( validPixel( pixelValue ) && ( !skipNodata || !isNoData ) )
145 if ( !pixelRectGeometry.
isNull() && polyEngine->intersects( pixelRectGeometry.
constGet() ) )
149 if ( !intersectGeometry.
isEmpty() )
151 const double intersectionArea = intersectGeometry.
area();
152 if ( intersectionArea > 0.0 )
154 weight = intersectionArea / pixelArea;
155 addValue( pixelValue, weight );
160 currentX += cellSizeX;
162 currentY -= cellSizeY;
167bool QgsRasterAnalysisUtils::validPixel(
double value )
169 return !std::isnan( value );
172void QgsRasterAnalysisUtils::mapToPixel(
const double x,
const double y,
const QgsRectangle bounds,
const double unitsPerPixelX,
const double unitsPerPixelY,
int &px,
int &py )
174 px = trunc( ( x - bounds.
xMinimum() ) / unitsPerPixelX );
175 py = trunc( ( y - bounds.
yMaximum() ) / -unitsPerPixelY );
178void QgsRasterAnalysisUtils::pixelToMap(
const int px,
const int py,
const QgsRectangle bounds,
const double unitsPerPixelX,
const double unitsPerPixelY,
double &x,
double &y )
180 x = bounds.
xMinimum() + ( px + 0.5 ) * unitsPerPixelX;
181 y = bounds.
yMaximum() - ( py + 0.5 ) * unitsPerPixelY;
184static QVector< QPair< QString, Qgis::DataType > > sDataTypes;
186void populateDataTypes()
188 if ( sDataTypes.empty() )
205std::unique_ptr<QgsProcessingParameterDefinition> QgsRasterAnalysisUtils::createRasterTypeParameter(
const QString &name,
const QString &description,
Qgis::DataType defaultType )
210 int defaultChoice = 0;
212 for (
auto it = sDataTypes.constBegin(); it != sDataTypes.constEnd(); ++it )
214 names.append( it->first );
215 if ( it->second == defaultType )
220 return std::make_unique< QgsProcessingParameterEnum >( name, description, names,
false, defaultChoice );
223Qgis::DataType QgsRasterAnalysisUtils::rasterTypeChoiceToDataType(
int choice )
225 if ( choice < 0 || choice >= sDataTypes.count() )
228 return sDataTypes.value( choice ).second;
231void QgsRasterAnalysisUtils::applyRasterLogicOperator(
const std::vector< QgsRasterAnalysisUtils::RasterLogicInput > &inputs,
QgsRasterDataProvider *destinationRaster,
double outputNoDataValue,
const bool treatNoDataAsFalse,
233 std::function<
void(
const std::vector< std::unique_ptr< QgsRasterBlock > > &,
bool &,
bool &,
int,
int,
bool )> &applyLogicFunc,
238 const int nbBlocksWidth =
static_cast< int>( std::ceil( 1.0 * width / maxWidth ) );
239 const int nbBlocksHeight =
static_cast< int >( std::ceil( 1.0 * height / maxHeight ) );
240 const int nbBlocks = nbBlocksWidth * nbBlocksHeight;
244 outputIter.startRasterRead( 1, width, height, extent );
251 std::unique_ptr< QgsRasterBlock > outputBlock;
252 while ( outputIter.readNextRasterPart( 1, iterCols, iterRows, outputBlock, iterLeft, iterTop, &blockExtent ) )
254 std::vector< std::unique_ptr< QgsRasterBlock > > inputBlocks;
255 for (
const QgsRasterAnalysisUtils::RasterLogicInput &i : inputs )
257 for (
const int band : i.bands )
259 std::unique_ptr< QgsRasterBlock > b( i.interface->block( band, blockExtent, iterCols, iterRows ) );
260 inputBlocks.emplace_back( std::move( b ) );
264 feedback->
setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
265 for (
int row = 0; row < iterRows; row++ )
270 for (
int column = 0; column < iterCols; column++ )
273 bool resIsNoData =
false;
274 applyLogicFunc( inputBlocks, res, resIsNoData, row, column, treatNoDataAsFalse );
282 outputBlock->setValue( row, column, resIsNoData ? outputNoDataValue : ( res ? 1 : 0 ) );
285 destinationRaster->
writeBlock( outputBlock.get(), 1, iterLeft, iterTop );
290std::vector<double> QgsRasterAnalysisUtils::getCellValuesFromBlockStack(
const std::vector< std::unique_ptr< QgsRasterBlock > > &inputBlocks,
int &row,
int &col,
bool &noDataInStack )
293 std::vector<double> cellValues;
294 bool hasNoData =
false;
295 cellValues.reserve( inputBlocks.size() );
297 for (
auto &block : inputBlocks )
300 if ( !block || !block->isValid() )
302 noDataInStack =
true;
307 value = block->valueAndNoData( row, col, hasNoData );
310 noDataInStack =
true;
315 cellValues.push_back( value );
322double QgsRasterAnalysisUtils::meanFromCellValues( std::vector<double> &cellValues,
int stackSize )
324 const double sum = std::accumulate( cellValues.begin(), cellValues.end(), 0.0 );
325 const double mean = sum /
static_cast<double>( stackSize );
329double QgsRasterAnalysisUtils::medianFromCellValues( std::vector<double> &cellValues,
int stackSize )
331 std::sort( cellValues.begin(), cellValues.end() );
332 const bool even = ( stackSize % 2 ) < 1;
335 return ( cellValues[stackSize / 2 - 1] + cellValues[stackSize / 2] ) / 2.0;
339 return cellValues[( stackSize + 1 ) / 2 - 1];
344double QgsRasterAnalysisUtils::stddevFromCellValues( std::vector<double> &cellValues,
int stackSize )
346 const double variance = varianceFromCellValues( cellValues, stackSize );
347 const double stddev = std::sqrt( variance );
351double QgsRasterAnalysisUtils::varianceFromCellValues( std::vector<double> &cellValues,
int stackSize )
353 const double mean = meanFromCellValues( cellValues, stackSize );
355 for (
int i = 0; i < stackSize; i++ )
357 accum += std::pow( ( cellValues.at( i ) - mean ), 2.0 );
359 const double variance = accum /
static_cast<double>( stackSize );
363double QgsRasterAnalysisUtils::maximumFromCellValues( std::vector<double> &cellValues )
365 return *std::max_element( cellValues.begin(), cellValues.end() );
368double QgsRasterAnalysisUtils::minimumFromCellValues( std::vector<double> &cellValues )
370 return *std::min_element( cellValues.begin(), cellValues.end() );
373double QgsRasterAnalysisUtils::majorityFromCellValues( std::vector<double> &cellValues,
const double noDataValue,
int stackSize )
375 if ( stackSize == 1 )
378 return cellValues[0];
380 else if ( stackSize == 2 )
383 return (
qgsDoubleNear( cellValues[0], cellValues[1] ) ) ? cellValues[0] : noDataValue;
385 else if ( std::adjacent_find( cellValues.begin(), cellValues.end(), std::not_equal_to<double>() ) == cellValues.end() )
389 return cellValues[0];
394 std::unordered_map<double, int> map;
396 for (
int i = 0; i < stackSize; i++ )
398 map[cellValues[i]]++;
402 bool multipleMajorities =
false;
403 double result = noDataValue;
404 for (
const auto &pair : std::as_const( map ) )
406 if ( maxCount < pair.second )
409 maxCount = pair.second;
410 multipleMajorities =
false;
412 else if ( maxCount == pair.second )
414 multipleMajorities =
true;
417 return multipleMajorities ? noDataValue : result;
421double QgsRasterAnalysisUtils::minorityFromCellValues( std::vector<double> &cellValues,
const double noDataValue,
int stackSize )
423 if ( stackSize == 1 )
426 return cellValues[0];
428 else if ( stackSize == 2 )
431 return (
qgsDoubleNear( cellValues[0], cellValues[1] ) ) ? cellValues[0] : noDataValue;
433 else if ( std::adjacent_find( cellValues.begin(), cellValues.end(), std::not_equal_to<double>() ) == cellValues.end() )
437 return cellValues[0];
442 std::unordered_map<double, int> map;
444 for (
int i = 0; i < stackSize; i++ )
446 map[cellValues[i]]++;
449 int minCount = stackSize;
450 bool multipleMinorities =
false;
451 double result = noDataValue;
452 for (
const auto &pair : std::as_const( map ) )
454 if ( minCount > pair.second )
457 minCount = pair.second;
458 multipleMinorities =
false;
460 else if ( minCount == pair.second )
462 multipleMinorities =
true;
465 return multipleMinorities ? noDataValue : result;
469double QgsRasterAnalysisUtils::rangeFromCellValues( std::vector<double> &cellValues )
471 const double max = *std::max_element( cellValues.begin(), cellValues.end() );
472 const double min = *std::min_element( cellValues.begin(), cellValues.end() );
476double QgsRasterAnalysisUtils::varietyFromCellValues( std::vector<double> &cellValues )
478 const std::unordered_set<double> uniqueValues( cellValues.begin(), cellValues.end() );
479 return uniqueValues.size();
482double QgsRasterAnalysisUtils::nearestRankPercentile( std::vector<double> &cellValues,
int stackSize,
double percentile )
485 std::sort( cellValues.begin(), cellValues.end() );
488 if ( percentile > 0 )
490 i = std::ceil( percentile *
static_cast<double>( stackSize ) ) - 1;
493 return cellValues[i];
496double QgsRasterAnalysisUtils::interpolatedPercentileInc( std::vector<double> &cellValues,
int stackSize,
double percentile )
498 std::sort( cellValues.begin(), cellValues.end() );
502 return cellValues[stackSize - 1 ];
506 return cellValues[0];
509 const double x = ( percentile * ( stackSize - 1 ) );
511 const int i =
static_cast<int>( std::floor( x ) );
512 const double xFraction = std::fmod( x, 1 );
514 if ( stackSize == 1 )
516 return cellValues[0];
518 else if ( stackSize == 2 )
520 return cellValues[0] + ( cellValues[1] - cellValues[0] ) * percentile;
524 return cellValues[i] + ( cellValues[i + 1] - cellValues[i] ) * xFraction;
528double QgsRasterAnalysisUtils::interpolatedPercentileExc( std::vector<double> &cellValues,
int stackSize,
double percentile,
double noDataValue )
530 std::sort( cellValues.begin(), cellValues.end() );
531 const double x = ( percentile * ( stackSize + 1 ) );
533 const int i =
static_cast<int>( std::floor( x ) ) - 1;
534 const double xFraction = std::fmod( x, 1 );
535 const double lowerExcValue = 1.0 / (
static_cast<double>( stackSize ) + 1.0 );
536 const double upperExcValue =
static_cast<double>( stackSize ) / (
static_cast<double>( stackSize ) + 1.0 );
538 if ( stackSize < 2 || ( ( percentile < lowerExcValue || percentile > upperExcValue ) ) )
544 return cellValues[i] + ( cellValues[i + 1] - cellValues[i] ) * xFraction;
548double QgsRasterAnalysisUtils::interpolatedPercentRankInc( std::vector<double> &cellValues,
int stackSize,
double value,
double noDataValue )
550 std::sort( cellValues.begin(), cellValues.end() );
552 if ( value < cellValues[0] || value > cellValues[stackSize - 1] )
558 for (
int i = 0; i < stackSize - 1; i++ )
560 if ( cellValues[i] <= value && cellValues[i + 1] >= value )
562 double fraction = 0.0;
566 fraction = ( value - cellValues[i] ) / ( cellValues[i + 1] - cellValues[i] );
568 return ( fraction + i ) / ( stackSize - 1 );
575double QgsRasterAnalysisUtils::interpolatedPercentRankExc( std::vector<double> &cellValues,
int stackSize,
double value,
double noDataValue )
577 std::sort( cellValues.begin(), cellValues.end() );
579 if ( value < cellValues[0] || value > cellValues[stackSize - 1] )
585 for (
int i = 0; i < stackSize - 1; i++ )
587 if ( cellValues[i] <= value && cellValues[i + 1] >= value )
589 double fraction = 0.0;
593 fraction = ( value - cellValues[i] ) / ( cellValues[i + 1] - cellValues[i] );
595 return ( ( i + 1 ) + fraction ) / ( stackSize + 1 );
DataType
Raster data types.
@ Float32
Thirty two bit floating point (float)
@ CFloat64
Complex Float64.
@ Int16
Sixteen bit signed integer (qint16)
@ Int8
Eight bit signed integer (qint8) (added in QGIS 3.30)
@ UInt16
Sixteen bit unsigned integer (quint16)
@ Byte
Eight bit unsigned integer (quint8)
@ Int32
Thirty two bit signed integer (qint32)
@ Float64
Sixty four bit floating point (double)
@ CFloat32
Complex Float32.
@ UInt32
Thirty two bit unsigned integer (quint32)
Base class for feedback objects to be used for cancellation of something running in a worker thread.
bool isCanceled() const
Tells whether the operation has been canceled already.
void setProgress(double progress)
Sets the current progress for the feedback object.
A geometry is the spatial representation of a feature.
static QgsGeometry fromRect(const QgsRectangle &rect)
Creates a new geometry from a QgsRectangle.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
double area() const
Returns the planar, 2-dimensional area of the geometry.
QgsGeometry intersection(const QgsGeometry &geometry, const QgsGeometryParameters ¶meters=QgsGeometryParameters()) const
Returns a geometry representing the points shared by this geometry and other.
bool isEmpty() const
Returns true if the geometry is empty (eg a linestring with no vertices, or a collection with no geom...
static QgsGeometryEngine * createGeometryEngine(const QgsAbstractGeometry *geometry, double precision=0.0)
Creates and returns a new geometry engine representing the specified geometry using precision on a gr...
Base class for raster data providers.
bool writeBlock(QgsRasterBlock *block, int band, int xOffset=0, int yOffset=0)
Writes pixel data from a raster block into the provider data source.
virtual bool setEditable(bool enabled)
Turns on/off editing mode of the provider.
Base class for processing filters like renderers, reprojector, resampler etc.
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.
A rectangle specified with double values.
double xMinimum() const
Returns the x minimum value (left side of rectangle).
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
double xMaximum() const
Returns the x maximum value (right side of rectangle).
double yMaximum() const
Returns the y maximum value (top side of rectangle).
bool isEmpty() const
Returns true if the rectangle has no area.
QgsRectangle intersect(const QgsRectangle &rect) const
Returns the intersection with the given rectangle.
unsigned long long qgssize
Qgssize is used instead of size_t, because size_t is stdlib type, unknown by SIP, and it would be har...
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)