26QString QgsRasterLayerZonalStatsAlgorithm::name()
const
28 return QStringLiteral(
"rasterlayerzonalstats" );
31QString QgsRasterLayerZonalStatsAlgorithm::displayName()
const
33 return QObject::tr(
"Raster layer zonal statistics" );
36QStringList QgsRasterLayerZonalStatsAlgorithm::tags()
const
38 return QObject::tr(
"count,area,statistics,stats,zones,categories,minimum,maximum,mean,sum,total" ).split(
',' );
41QString QgsRasterLayerZonalStatsAlgorithm::group()
const
43 return QObject::tr(
"Raster analysis" );
46QString QgsRasterLayerZonalStatsAlgorithm::groupId()
const
48 return QStringLiteral(
"rasteranalysis" );
51void QgsRasterLayerZonalStatsAlgorithm::initAlgorithm(
const QVariantMap & )
54 addParameter(
new QgsProcessingParameterBand( QStringLiteral(
"BAND" ), QObject::tr(
"Band number" ), 1, QStringLiteral(
"INPUT" ) ) );
56 addParameter(
new QgsProcessingParameterBand( QStringLiteral(
"ZONES_BAND" ), QObject::tr(
"Zones band number" ), 1, QStringLiteral(
"ZONES" ) ) );
58 std::unique_ptr<QgsProcessingParameterEnum> refParam = std::make_unique<QgsProcessingParameterEnum>( QStringLiteral(
"REF_LAYER" ), QObject::tr(
"Reference layer" ), QStringList() << QObject::tr(
"Input layer" ) << QObject::tr(
"Zones layer" ),
false, 0 );
60 addParameter( refParam.release() );
69 addOutput(
new QgsProcessingOutputNumber( QStringLiteral(
"NODATA_PIXEL_COUNT" ), QObject::tr(
"NoData pixel count" ) ) );
72QString QgsRasterLayerZonalStatsAlgorithm::shortDescription()
const
74 return QObject::tr(
"Calculates statistics for a raster layer's values, categorized by zones defined in another raster layer." );
77QString QgsRasterLayerZonalStatsAlgorithm::shortHelpString()
const
79 return QObject::tr(
"This algorithm calculates statistics for a raster layer's values, categorized by zones defined in another raster layer.\n\n"
80 "If the reference layer parameter is set to \"Input layer\", then zones are determined by sampling the zone raster layer value at the centroid of each pixel from the source raster layer.\n\n"
81 "If the reference layer parameter is set to \"Zones layer\", then the input raster layer will be sampled at the centroid of each pixel from the zones raster layer.\n\n"
82 "If either the source raster layer or the zone raster layer value is NoData for a pixel, that pixel's value will be skipped and not included in the calculated statistics." );
85QgsRasterLayerZonalStatsAlgorithm *QgsRasterLayerZonalStatsAlgorithm::createInstance()
const
87 return new QgsRasterLayerZonalStatsAlgorithm();
92 mRefLayer =
static_cast<RefLayer
>( parameterAsEnum( parameters, QStringLiteral(
"REF_LAYER" ), context ) );
94 QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral(
"INPUT" ), context );
95 const int band = parameterAsInt( parameters, QStringLiteral(
"BAND" ), context );
100 mBand = parameterAsInt( parameters, QStringLiteral(
"BAND" ), context );
101 if ( mBand < 1 || mBand > layer->
bandCount() )
102 throw QgsProcessingException( QObject::tr(
"Invalid band number for BAND (%1): Valid values for input raster are 1 to %2" ).arg( mBand ).arg( layer->
bandCount() ) );
106 QgsRasterLayer *zonesLayer = parameterAsRasterLayer( parameters, QStringLiteral(
"ZONES" ), context );
111 mZonesBand = parameterAsInt( parameters, QStringLiteral(
"ZONES_BAND" ), context );
112 if ( mZonesBand < 1 || mZonesBand > zonesLayer->
bandCount() )
113 throw QgsProcessingException( QObject::tr(
"Invalid band number for ZONES_BAND (%1): Valid values for input raster are 1 to %2" ).arg( mZonesBand ).arg( zonesLayer->
bandCount() ) );
117 mSourceInterface = mSourceDataProvider.get();
119 mZonesInterface = mZonesDataProvider.get();
127 mLayerWidth = layer->
width();
128 mLayerHeight = layer->
height();
129 mExtent = layer->
extent();
132 if ( layer->
crs() != zonesLayer->
crs() )
134 mProjector = std::make_unique<QgsRasterProjector>();
135 mProjector->setInput( mZonesDataProvider.get() );
137 mZonesInterface = mProjector.get();
142 mCrs = zonesLayer->
crs();
145 mLayerWidth = zonesLayer->
width();
146 mLayerHeight = zonesLayer->
height();
147 mExtent = zonesLayer->
extent();
150 if ( layer->
crs() != zonesLayer->
crs() )
152 mProjector = std::make_unique<QgsRasterProjector>();
153 mProjector->setInput( mSourceDataProvider.get() );
155 mSourceInterface = mProjector.get();
168 std::unique_ptr<QgsFeatureSink> sink;
169 if ( parameters.contains( QStringLiteral(
"OUTPUT_TABLE" ) ) && parameters.value( QStringLiteral(
"OUTPUT_TABLE" ) ).isValid() )
172 outFields.
append(
QgsField( QStringLiteral(
"zone" ), QMetaType::Type::Double, QString(), 20, 8 ) );
173 outFields.
append(
QgsField( areaUnit.isEmpty() ?
"area" : areaUnit.replace( QStringLiteral(
"²" ), QStringLiteral(
"2" ) ), QMetaType::Type::Double, QString(), 20, 8 ) );
174 outFields.
append(
QgsField( QStringLiteral(
"sum" ), QMetaType::Type::Double, QString(), 20, 8 ) );
175 outFields.
append(
QgsField( QStringLiteral(
"count" ), QMetaType::Type::LongLong, QString(), 20 ) );
176 outFields.
append(
QgsField( QStringLiteral(
"min" ), QMetaType::Type::Double, QString(), 20, 8 ) );
177 outFields.
append(
QgsField( QStringLiteral(
"max" ), QMetaType::Type::Double, QString(), 20, 8 ) );
178 outFields.
append(
QgsField( QStringLiteral(
"mean" ), QMetaType::Type::Double, QString(), 20, 8 ) );
185 struct StatCalculator
191 QHash<double, StatCalculator> zoneStats;
194 const qgssize layerSize =
static_cast<qgssize>( mLayerWidth ) *
static_cast<qgssize>( mLayerHeight );
197 const int nbBlocksWidth =
static_cast<int>( std::ceil( 1.0 * mLayerWidth / maxWidth ) );
198 const int nbBlocksHeight =
static_cast<int>( std::ceil( 1.0 * mLayerHeight / maxHeight ) );
199 const int nbBlocks = nbBlocksWidth * nbBlocksHeight;
203 iter.
startRasterRead( mRefLayer == Source ? mBand : mZonesBand, mLayerWidth, mLayerHeight, mExtent );
210 std::unique_ptr<QgsRasterBlock> rasterBlock;
211 std::unique_ptr<QgsRasterBlock> zonesRasterBlock;
212 bool isNoData =
false;
215 if ( mRefLayer == Source )
217 if ( !iter.
readNextRasterPart( mBand, iterCols, iterRows, rasterBlock, iterLeft, iterTop, &blockExtent ) )
220 zonesRasterBlock.reset( mZonesInterface->block( mZonesBand, blockExtent, iterCols, iterRows ) );
224 if ( !iter.
readNextRasterPart( mZonesBand, iterCols, iterRows, zonesRasterBlock, iterLeft, iterTop, &blockExtent ) )
227 rasterBlock.reset( mSourceInterface->block( mBand, blockExtent, iterCols, iterRows ) );
229 if ( !zonesRasterBlock || !rasterBlock )
232 feedback->
setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
233 if ( !rasterBlock->isValid() || rasterBlock->isEmpty() || !zonesRasterBlock->isValid() || zonesRasterBlock->isEmpty() )
236 for (
int row = 0; row < iterRows; row++ )
241 for (
int column = 0; column < iterCols; column++ )
243 const double value = rasterBlock->valueAndNoData( row, column, isNoData );
244 if ( mHasNoDataValue && isNoData )
249 const double zone = zonesRasterBlock->valueAndNoData( row, column, isNoData );
250 if ( mZonesHasNoDataValue && isNoData )
255 zoneStats[zone].s.addValue( value );
261 outputs.insert( QStringLiteral(
"EXTENT" ), mExtent.toString() );
262 outputs.insert( QStringLiteral(
"CRS_AUTHID" ), mCrs.authid() );
263 outputs.insert( QStringLiteral(
"WIDTH_IN_PIXELS" ), mLayerWidth );
264 outputs.insert( QStringLiteral(
"HEIGHT_IN_PIXELS" ), mLayerHeight );
265 outputs.insert( QStringLiteral(
"TOTAL_PIXEL_COUNT" ), layerSize );
266 outputs.insert( QStringLiteral(
"NODATA_PIXEL_COUNT" ), noDataCount );
268 const double pixelArea = mRasterUnitsPerPixelX * mRasterUnitsPerPixelY;
270 for (
auto it = zoneStats.begin(); it != zoneStats.end(); ++it )
274 f.
setAttributes(
QgsAttributes() << it.key() << it->s.count() * pixelArea << it->s.sum() << it->s.count() << it->s.min() << it->s.max() << it->s.mean() );
276 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral(
"OUTPUT_TABLE" ) ) );
279 outputs.insert( QStringLiteral(
"OUTPUT_TABLE" ), tableDest );
@ Vector
Tables (i.e. vector layers with or without geometry). When used for a sink this indicates the sink ha...
@ Advanced
Parameter is an advanced parameter which should be hidden from users by default.
This class represents a coordinate reference system (CRS).
@ FastInsert
Use faster inserts, at the cost of updating the passed features to reflect changes made at the provid...
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
void setAttributes(const QgsAttributes &attrs)
Sets the feature's attributes.
bool isCanceled() const
Tells whether the operation has been canceled already.
void setProgress(double progress)
Sets the current progress for the feedback object.
Encapsulate a field in an attribute table or data source.
Container of fields for a vector layer.
bool append(const QgsField &field, Qgis::FieldOrigin origin=Qgis::FieldOrigin::Provider, int originIndex=-1)
Appends a field.
virtual QgsRectangle extent() const
Returns the extent of the layer.
QgsCoordinateReferenceSystem crs
Contains information about the context in which a processing algorithm is executed.
QgsCoordinateTransformContext transformContext() const
Returns the coordinate transform context.
Custom exception class for processing related exceptions.
Base class for providing feedback from a processing algorithm.
A numeric output for processing algorithms.
A string output for processing algorithms.
A raster band parameter for Processing algorithms.
A feature sink output 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.
Iterator for sequentially processing raster cells.
static const int DEFAULT_MAXIMUM_TILE_WIDTH
Default maximum tile width.
bool readNextRasterPart(int bandNumber, int &nCols, int &nRows, QgsRasterBlock **block, int &topLeftCol, int &topLeftRow)
Fetches next part of raster data, caller takes ownership of the block and caller should delete the bl...
static const int DEFAULT_MAXIMUM_TILE_HEIGHT
Default maximum tile height.
void startRasterRead(int bandNumber, qgssize nCols, qgssize nRows, const QgsRectangle &extent, QgsRasterBlockFeedback *feedback=nullptr)
Start reading of raster band.
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.
Calculator for summary statistics for a list of doubles.
static Q_INVOKABLE QString toAbbreviatedString(Qgis::DistanceUnit unit)
Returns a translated abbreviation representing a distance unit.
static Q_INVOKABLE Qgis::AreaUnit distanceToAreaUnit(Qgis::DistanceUnit distanceUnit)
Converts a distance unit to its corresponding area unit, e.g., meters to square meters.
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...