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 QObject::tr(
"Input layer" ) ) );
56 QObject::tr(
"Band number" ), 1, QStringLiteral(
"INPUT" ) ) );
58 QObject::tr(
"Zones layer" ) ) );
60 QObject::tr(
"Zones band number" ), 1, QStringLiteral(
"ZONES" ) ) );
62 std::unique_ptr< QgsProcessingParameterEnum > refParam = std::make_unique< QgsProcessingParameterEnum >( QStringLiteral(
"REF_LAYER" ), QObject::tr(
"Reference layer" ),
63 QStringList() << QObject::tr(
"Input layer" ) << QObject::tr(
"Zones layer" ),
false, 0 );
65 addParameter( refParam.release() );
75 addOutput(
new QgsProcessingOutputNumber( QStringLiteral(
"NODATA_PIXEL_COUNT" ), QObject::tr(
"NoData pixel count" ) ) );
78QString QgsRasterLayerZonalStatsAlgorithm::shortDescription()
const
80 return QObject::tr(
"Calculates statistics for a raster layer's values, categorized by zones defined in another raster layer." );
83QString QgsRasterLayerZonalStatsAlgorithm::shortHelpString()
const
85 return QObject::tr(
"This algorithm calculates statistics for a raster layer's values, categorized by zones defined in another raster layer.\n\n"
86 "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"
87 "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"
88 "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." );
91QgsRasterLayerZonalStatsAlgorithm *QgsRasterLayerZonalStatsAlgorithm::createInstance()
const
93 return new QgsRasterLayerZonalStatsAlgorithm();
98 mRefLayer =
static_cast< RefLayer
>( parameterAsEnum( parameters, QStringLiteral(
"REF_LAYER" ), context ) );
100 QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral(
"INPUT" ), context );
101 const int band = parameterAsInt( parameters, QStringLiteral(
"BAND" ), context );
106 mBand = parameterAsInt( parameters, QStringLiteral(
"BAND" ), context );
107 if ( mBand < 1 || mBand > layer->
bandCount() )
108 throw QgsProcessingException( QObject::tr(
"Invalid band number for BAND (%1): Valid values for input raster are 1 to %2" ).arg( mBand )
113 QgsRasterLayer *zonesLayer = parameterAsRasterLayer( parameters, QStringLiteral(
"ZONES" ), context );
118 mZonesBand = parameterAsInt( parameters, QStringLiteral(
"ZONES_BAND" ), context );
119 if ( mZonesBand < 1 || mZonesBand > zonesLayer->
bandCount() )
120 throw QgsProcessingException( QObject::tr(
"Invalid band number for ZONES_BAND (%1): Valid values for input raster are 1 to %2" ).arg( mZonesBand )
125 mSourceInterface = mSourceDataProvider.get();
127 mZonesInterface = mZonesDataProvider.get();
135 mLayerWidth = layer->
width();
136 mLayerHeight = layer->
height();
137 mExtent = layer->
extent();
140 if ( layer->
crs() != zonesLayer->
crs() )
142 mProjector = std::make_unique< QgsRasterProjector >();
143 mProjector->setInput( mZonesDataProvider.get() );
145 mZonesInterface = mProjector.get();
150 mCrs = zonesLayer->
crs();
153 mLayerWidth = zonesLayer->
width();
154 mLayerHeight = zonesLayer->
height();
155 mExtent = zonesLayer->
extent();
158 if ( layer->
crs() != zonesLayer->
crs() )
160 mProjector = std::make_unique< QgsRasterProjector >();
161 mProjector->setInput( mSourceDataProvider.get() );
163 mSourceInterface = mProjector.get();
176 std::unique_ptr< QgsFeatureSink > sink;
177 if ( parameters.contains( QStringLiteral(
"OUTPUT_TABLE" ) ) && parameters.value( QStringLiteral(
"OUTPUT_TABLE" ) ).isValid() )
180 outFields.
append(
QgsField( QStringLiteral(
"zone" ), QMetaType::Type::Double, QString(), 20, 8 ) );
181 outFields.
append(
QgsField( areaUnit.isEmpty() ?
"area" : areaUnit.replace( QStringLiteral(
"²" ), QStringLiteral(
"2" ) ), QMetaType::Type::Double, QString(), 20, 8 ) );
182 outFields.
append(
QgsField( QStringLiteral(
"sum" ), QMetaType::Type::Double, QString(), 20, 8 ) );
183 outFields.
append(
QgsField( QStringLiteral(
"count" ), QMetaType::Type::LongLong, QString(), 20 ) );
184 outFields.
append(
QgsField( QStringLiteral(
"min" ), QMetaType::Type::Double, QString(), 20, 8 ) );
185 outFields.
append(
QgsField( QStringLiteral(
"max" ), QMetaType::Type::Double, QString(), 20, 8 ) );
186 outFields.
append(
QgsField( QStringLiteral(
"mean" ), QMetaType::Type::Double, QString(), 20, 8 ) );
193 struct StatCalculator
199 QHash<double, StatCalculator > zoneStats;
202 const qgssize layerSize =
static_cast< qgssize >( mLayerWidth ) *
static_cast< qgssize >( mLayerHeight );
205 const int nbBlocksWidth =
static_cast< int>( std::ceil( 1.0 * mLayerWidth / maxWidth ) );
206 const int nbBlocksHeight =
static_cast< int >( std::ceil( 1.0 * mLayerHeight / maxHeight ) );
207 const int nbBlocks = nbBlocksWidth * nbBlocksHeight;
211 iter.
startRasterRead( mRefLayer == Source ? mBand : mZonesBand, mLayerWidth, mLayerHeight, mExtent );
218 std::unique_ptr< QgsRasterBlock > rasterBlock;
219 std::unique_ptr< QgsRasterBlock > zonesRasterBlock;
220 bool isNoData =
false;
223 if ( mRefLayer == Source )
225 if ( !iter.
readNextRasterPart( mBand, iterCols, iterRows, rasterBlock, iterLeft, iterTop, &blockExtent ) )
228 zonesRasterBlock.reset( mZonesInterface->block( mZonesBand, blockExtent, iterCols, iterRows ) );
232 if ( !iter.
readNextRasterPart( mZonesBand, iterCols, iterRows, zonesRasterBlock, iterLeft, iterTop, &blockExtent ) )
235 rasterBlock.reset( mSourceInterface->block( mBand, blockExtent, iterCols, iterRows ) );
237 if ( !zonesRasterBlock || !rasterBlock )
240 feedback->
setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
241 if ( !rasterBlock->isValid() || rasterBlock->isEmpty() || !zonesRasterBlock->isValid() || zonesRasterBlock->isEmpty() )
244 for (
int row = 0; row < iterRows; row++ )
249 for (
int column = 0; column < iterCols; column++ )
251 const double value = rasterBlock->valueAndNoData( row, column, isNoData );
252 if ( mHasNoDataValue && isNoData )
257 const double zone = zonesRasterBlock->valueAndNoData( row, column, isNoData );
258 if ( mZonesHasNoDataValue && isNoData )
263 zoneStats[ zone ].s.addValue( value );
269 outputs.insert( QStringLiteral(
"EXTENT" ), mExtent.toString() );
270 outputs.insert( QStringLiteral(
"CRS_AUTHID" ), mCrs.authid() );
271 outputs.insert( QStringLiteral(
"WIDTH_IN_PIXELS" ), mLayerWidth );
272 outputs.insert( QStringLiteral(
"HEIGHT_IN_PIXELS" ), mLayerHeight );
273 outputs.insert( QStringLiteral(
"TOTAL_PIXEL_COUNT" ), layerSize );
274 outputs.insert( QStringLiteral(
"NODATA_PIXEL_COUNT" ), noDataCount );
276 const double pixelArea = mRasterUnitsPerPixelX * mRasterUnitsPerPixelY;
278 for (
auto it = zoneStats.begin(); it != zoneStats.end(); ++it )
283 it->s.min() << it->s.max() << it->s.mean() );
285 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral(
"OUTPUT_TABLE" ) ) );
287 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...