Branch data Line data Source code
1 : : /***************************************************************************
2 : : qgsalgorithmrasterzonalstats.cpp
3 : : ---------------------
4 : : begin : December 2018
5 : : copyright : (C) 2018 by Nyall Dawson
6 : : email : nyall dot dawson 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 "qgsalgorithmrasterzonalstats.h"
19 : : #include "qgsstringutils.h"
20 : : #include "qgsstatisticalsummary.h"
21 : : #include "qgsrasterprojector.h"
22 : :
23 : : ///@cond PRIVATE
24 : :
25 : 0 : QString QgsRasterLayerZonalStatsAlgorithm::name() const
26 : : {
27 : 0 : return QStringLiteral( "rasterlayerzonalstats" );
28 : : }
29 : :
30 : 0 : QString QgsRasterLayerZonalStatsAlgorithm::displayName() const
31 : : {
32 : 0 : return QObject::tr( "Raster layer zonal statistics" );
33 : : }
34 : :
35 : 0 : QStringList QgsRasterLayerZonalStatsAlgorithm::tags() const
36 : : {
37 : 0 : return QObject::tr( "count,area,statistics,stats,zones,categories,minimum,maximum,mean,sum,total" ).split( ',' );
38 : 0 : }
39 : :
40 : 0 : QString QgsRasterLayerZonalStatsAlgorithm::group() const
41 : : {
42 : 0 : return QObject::tr( "Raster analysis" );
43 : : }
44 : :
45 : 0 : QString QgsRasterLayerZonalStatsAlgorithm::groupId() const
46 : : {
47 : 0 : return QStringLiteral( "rasteranalysis" );
48 : : }
49 : :
50 : 0 : void QgsRasterLayerZonalStatsAlgorithm::initAlgorithm( const QVariantMap & )
51 : : {
52 : 0 : addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT" ),
53 : 0 : QObject::tr( "Input layer" ) ) );
54 : 0 : addParameter( new QgsProcessingParameterBand( QStringLiteral( "BAND" ),
55 : 0 : QObject::tr( "Band number" ), 1, QStringLiteral( "INPUT" ) ) );
56 : 0 : addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "ZONES" ),
57 : 0 : QObject::tr( "Zones layer" ) ) );
58 : 0 : addParameter( new QgsProcessingParameterBand( QStringLiteral( "ZONES_BAND" ),
59 : 0 : QObject::tr( "Zones band number" ), 1, QStringLiteral( "ZONES" ) ) );
60 : :
61 : 0 : std::unique_ptr< QgsProcessingParameterEnum > refParam = std::make_unique< QgsProcessingParameterEnum >( QStringLiteral( "REF_LAYER" ), QObject::tr( "Reference layer" ),
62 : 0 : QStringList() << QObject::tr( "Input layer" ) << QObject::tr( "Zones layer" ), false, 0 );
63 : 0 : refParam->setFlags( refParam->flags() | QgsProcessingParameterDefinition::FlagAdvanced );
64 : 0 : addParameter( refParam.release() );
65 : :
66 : 0 : addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT_TABLE" ),
67 : 0 : QObject::tr( "Statistics" ), QgsProcessing::TypeVector ) );
68 : :
69 : 0 : addOutput( new QgsProcessingOutputString( QStringLiteral( "EXTENT" ), QObject::tr( "Extent" ) ) );
70 : 0 : addOutput( new QgsProcessingOutputString( QStringLiteral( "CRS_AUTHID" ), QObject::tr( "CRS authority identifier" ) ) );
71 : 0 : addOutput( new QgsProcessingOutputNumber( QStringLiteral( "WIDTH_IN_PIXELS" ), QObject::tr( "Width in pixels" ) ) );
72 : 0 : addOutput( new QgsProcessingOutputNumber( QStringLiteral( "HEIGHT_IN_PIXELS" ), QObject::tr( "Height in pixels" ) ) );
73 : 0 : addOutput( new QgsProcessingOutputNumber( QStringLiteral( "TOTAL_PIXEL_COUNT" ), QObject::tr( "Total pixel count" ) ) );
74 : 0 : addOutput( new QgsProcessingOutputNumber( QStringLiteral( "NODATA_PIXEL_COUNT" ), QObject::tr( "NODATA pixel count" ) ) );
75 : 0 : }
76 : :
77 : 0 : QString QgsRasterLayerZonalStatsAlgorithm::shortDescription() const
78 : : {
79 : 0 : return QObject::tr( "Calculates statistics for a raster layer's values, categorized by zones defined in another raster layer." );
80 : : }
81 : :
82 : 0 : QString QgsRasterLayerZonalStatsAlgorithm::shortHelpString() const
83 : : {
84 : 0 : return QObject::tr( "This algorithm calculates statistics for a raster layer's values, categorized by zones defined in another raster layer.\n\n"
85 : : "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"
86 : : "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"
87 : : "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 including in the calculated statistics." );
88 : : }
89 : :
90 : 0 : QgsRasterLayerZonalStatsAlgorithm *QgsRasterLayerZonalStatsAlgorithm::createInstance() const
91 : : {
92 : 0 : return new QgsRasterLayerZonalStatsAlgorithm();
93 : 0 : }
94 : :
95 : 0 : bool QgsRasterLayerZonalStatsAlgorithm::prepareAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback * )
96 : : {
97 : 0 : mRefLayer = static_cast< RefLayer >( parameterAsEnum( parameters, QStringLiteral( "REF_LAYER" ), context ) );
98 : :
99 : 0 : QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT" ), context );
100 : 0 : int band = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
101 : :
102 : 0 : if ( !layer )
103 : 0 : throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT" ) ) );
104 : :
105 : 0 : mBand = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
106 : 0 : if ( mBand < 1 || mBand > layer->bandCount() )
107 : 0 : throw QgsProcessingException( QObject::tr( "Invalid band number for BAND (%1): Valid values for input raster are 1 to %2" ).arg( mBand )
108 : 0 : .arg( layer->bandCount() ) );
109 : :
110 : 0 : mHasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( band );
111 : :
112 : 0 : QgsRasterLayer *zonesLayer = parameterAsRasterLayer( parameters, QStringLiteral( "ZONES" ), context );
113 : :
114 : 0 : if ( !zonesLayer )
115 : 0 : throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "ZONES" ) ) );
116 : :
117 : 0 : mZonesBand = parameterAsInt( parameters, QStringLiteral( "ZONES_BAND" ), context );
118 : 0 : if ( mZonesBand < 1 || mZonesBand > zonesLayer->bandCount() )
119 : 0 : throw QgsProcessingException( QObject::tr( "Invalid band number for ZONES_BAND (%1): Valid values for input raster are 1 to %2" ).arg( mZonesBand )
120 : 0 : .arg( zonesLayer->bandCount() ) );
121 : 0 : mZonesHasNoDataValue = zonesLayer->dataProvider()->sourceHasNoDataValue( band );
122 : :
123 : 0 : mSourceDataProvider.reset( layer->dataProvider()->clone() );
124 : 0 : mSourceInterface = mSourceDataProvider.get();
125 : 0 : mZonesDataProvider.reset( zonesLayer->dataProvider()->clone() );
126 : 0 : mZonesInterface = mZonesDataProvider.get();
127 : :
128 : 0 : switch ( mRefLayer )
129 : : {
130 : : case Source:
131 : 0 : mCrs = layer->crs();
132 : 0 : mRasterUnitsPerPixelX = layer->rasterUnitsPerPixelX();
133 : 0 : mRasterUnitsPerPixelY = layer->rasterUnitsPerPixelY();
134 : 0 : mLayerWidth = layer->width();
135 : 0 : mLayerHeight = layer->height();
136 : 0 : mExtent = layer->extent();
137 : :
138 : : // add projector if necessary
139 : 0 : if ( layer->crs() != zonesLayer->crs() )
140 : : {
141 : 0 : mProjector = std::make_unique< QgsRasterProjector >();
142 : 0 : mProjector->setInput( mZonesDataProvider.get() );
143 : 0 : mProjector->setCrs( zonesLayer->crs(), layer->crs(), context.transformContext() );
144 : 0 : mZonesInterface = mProjector.get();
145 : 0 : }
146 : 0 : break;
147 : :
148 : : case Zones:
149 : 0 : mCrs = zonesLayer->crs();
150 : 0 : mRasterUnitsPerPixelX = zonesLayer->rasterUnitsPerPixelX();
151 : 0 : mRasterUnitsPerPixelY = zonesLayer->rasterUnitsPerPixelY();
152 : 0 : mLayerWidth = zonesLayer->width();
153 : 0 : mLayerHeight = zonesLayer->height();
154 : 0 : mExtent = zonesLayer->extent();
155 : :
156 : : // add projector if necessary
157 : 0 : if ( layer->crs() != zonesLayer->crs() )
158 : : {
159 : 0 : mProjector = std::make_unique< QgsRasterProjector >();
160 : 0 : mProjector->setInput( mSourceDataProvider.get() );
161 : 0 : mProjector->setCrs( layer->crs(), zonesLayer->crs(), context.transformContext() );
162 : 0 : mSourceInterface = mProjector.get();
163 : 0 : }
164 : 0 : break;
165 : : }
166 : :
167 : 0 : return true;
168 : 0 : }
169 : :
170 : 0 : QVariantMap QgsRasterLayerZonalStatsAlgorithm::processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
171 : : {
172 : 0 : QString areaUnit = QgsUnitTypes::toAbbreviatedString( QgsUnitTypes::distanceToAreaUnit( mCrs.mapUnits() ) );
173 : :
174 : 0 : QString tableDest;
175 : 0 : std::unique_ptr< QgsFeatureSink > sink;
176 : 0 : if ( parameters.contains( QStringLiteral( "OUTPUT_TABLE" ) ) && parameters.value( QStringLiteral( "OUTPUT_TABLE" ) ).isValid() )
177 : : {
178 : 0 : QgsFields outFields;
179 : 0 : outFields.append( QgsField( QStringLiteral( "zone" ), QVariant::Double, QString(), 20, 8 ) );
180 : 0 : outFields.append( QgsField( areaUnit.replace( QStringLiteral( "²" ), QStringLiteral( "2" ) ), QVariant::Double, QString(), 20, 8 ) );
181 : 0 : outFields.append( QgsField( QStringLiteral( "sum" ), QVariant::Double, QString(), 20, 8 ) );
182 : 0 : outFields.append( QgsField( QStringLiteral( "count" ), QVariant::LongLong, QString(), 20 ) );
183 : 0 : outFields.append( QgsField( QStringLiteral( "min" ), QVariant::Double, QString(), 20, 8 ) );
184 : 0 : outFields.append( QgsField( QStringLiteral( "max" ), QVariant::Double, QString(), 20, 8 ) );
185 : 0 : outFields.append( QgsField( QStringLiteral( "mean" ), QVariant::Double, QString(), 20, 8 ) );
186 : :
187 : 0 : sink.reset( parameterAsSink( parameters, QStringLiteral( "OUTPUT_TABLE" ), context, tableDest, outFields, QgsWkbTypes::NoGeometry, QgsCoordinateReferenceSystem() ) );
188 : 0 : if ( !sink )
189 : 0 : throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT_TABLE" ) ) );
190 : 0 : }
191 : :
192 : 0 : struct StatCalculator
193 : : {
194 : : // only calculate cheap stats-- we cannot calculate stats which require holding values in memory -- because otherwise we'll end
195 : : // up trying to store EVERY pixel value from the input in memory
196 : 0 : QgsStatisticalSummary s{ QgsStatisticalSummary::Count | QgsStatisticalSummary::Sum | QgsStatisticalSummary::Min | QgsStatisticalSummary::Max | QgsStatisticalSummary::Mean };
197 : : };
198 : 0 : QHash<double, StatCalculator > zoneStats;
199 : 0 : qgssize noDataCount = 0;
200 : :
201 : 0 : qgssize layerSize = static_cast< qgssize >( mLayerWidth ) * static_cast< qgssize >( mLayerHeight );
202 : 0 : int maxWidth = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_WIDTH;
203 : 0 : int maxHeight = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_HEIGHT;
204 : 0 : int nbBlocksWidth = static_cast< int>( std::ceil( 1.0 * mLayerWidth / maxWidth ) );
205 : 0 : int nbBlocksHeight = static_cast< int >( std::ceil( 1.0 * mLayerHeight / maxHeight ) );
206 : 0 : int nbBlocks = nbBlocksWidth * nbBlocksHeight;
207 : :
208 : 0 : QgsRasterIterator iter = mRefLayer == Source ? QgsRasterIterator( mSourceInterface )
209 : 0 : : QgsRasterIterator( mZonesInterface );
210 : 0 : iter.startRasterRead( mRefLayer == Source ? mBand : mZonesBand, mLayerWidth, mLayerHeight, mExtent );
211 : :
212 : 0 : int iterLeft = 0;
213 : 0 : int iterTop = 0;
214 : 0 : int iterCols = 0;
215 : 0 : int iterRows = 0;
216 : 0 : QgsRectangle blockExtent;
217 : 0 : std::unique_ptr< QgsRasterBlock > rasterBlock;
218 : 0 : std::unique_ptr< QgsRasterBlock > zonesRasterBlock;
219 : 0 : bool isNoData = false;
220 : 0 : while ( true )
221 : : {
222 : 0 : if ( mRefLayer == Source )
223 : : {
224 : 0 : if ( !iter.readNextRasterPart( mBand, iterCols, iterRows, rasterBlock, iterLeft, iterTop, &blockExtent ) )
225 : 0 : break;
226 : :
227 : 0 : zonesRasterBlock.reset( mZonesInterface->block( mZonesBand, blockExtent, iterCols, iterRows ) );
228 : 0 : }
229 : : else
230 : : {
231 : 0 : if ( !iter.readNextRasterPart( mZonesBand, iterCols, iterRows, zonesRasterBlock, iterLeft, iterTop, &blockExtent ) )
232 : 0 : break;
233 : :
234 : 0 : rasterBlock.reset( mSourceInterface->block( mBand, blockExtent, iterCols, iterRows ) );
235 : : }
236 : 0 : if ( !zonesRasterBlock || !rasterBlock )
237 : 0 : continue;
238 : :
239 : 0 : feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
240 : 0 : if ( !rasterBlock->isValid() || rasterBlock->isEmpty() || !zonesRasterBlock->isValid() || zonesRasterBlock->isEmpty() )
241 : 0 : continue;
242 : :
243 : 0 : for ( int row = 0; row < iterRows; row++ )
244 : : {
245 : 0 : if ( feedback->isCanceled() )
246 : 0 : break;
247 : :
248 : 0 : for ( int column = 0; column < iterCols; column++ )
249 : : {
250 : 0 : double value = rasterBlock->valueAndNoData( row, column, isNoData );
251 : 0 : if ( mHasNoDataValue && isNoData )
252 : : {
253 : 0 : noDataCount += 1;
254 : 0 : continue;
255 : : }
256 : 0 : double zone = zonesRasterBlock->valueAndNoData( row, column, isNoData );
257 : 0 : if ( mZonesHasNoDataValue && isNoData )
258 : : {
259 : 0 : noDataCount += 1;
260 : 0 : continue;
261 : : }
262 : 0 : zoneStats[ zone ].s.addValue( value );
263 : 0 : }
264 : 0 : }
265 : : }
266 : :
267 : 0 : QVariantMap outputs;
268 : 0 : outputs.insert( QStringLiteral( "EXTENT" ), mExtent.toString() );
269 : 0 : outputs.insert( QStringLiteral( "CRS_AUTHID" ), mCrs.authid() );
270 : 0 : outputs.insert( QStringLiteral( "WIDTH_IN_PIXELS" ), mLayerWidth );
271 : 0 : outputs.insert( QStringLiteral( "HEIGHT_IN_PIXELS" ), mLayerHeight );
272 : 0 : outputs.insert( QStringLiteral( "TOTAL_PIXEL_COUNT" ), layerSize );
273 : 0 : outputs.insert( QStringLiteral( "NODATA_PIXEL_COUNT" ), noDataCount );
274 : :
275 : 0 : double pixelArea = mRasterUnitsPerPixelX * mRasterUnitsPerPixelY;
276 : :
277 : 0 : for ( auto it = zoneStats.begin(); it != zoneStats.end(); ++it )
278 : : {
279 : 0 : QgsFeature f;
280 : 0 : it->s.finalize();
281 : 0 : f.setAttributes( QgsAttributes() << it.key() << it->s.count() * pixelArea << it->s.sum() << it->s.count() <<
282 : 0 : it->s.min() << it->s.max() << it->s.mean() );
283 : 0 : sink->addFeature( f, QgsFeatureSink::FastInsert );
284 : 0 : }
285 : 0 : outputs.insert( QStringLiteral( "OUTPUT_TABLE" ), tableDest );
286 : :
287 : 0 : return outputs;
288 : 0 : }
289 : :
290 : :
291 : : ///@endcond
292 : :
293 : :
294 : :
|
