Branch data Line data Source code
1 : : /***************************************************************************
2 : : qgsalgorithmzonalhistogram.cpp
3 : : ---------------------
4 : : begin : May, 2018
5 : : copyright : (C) 2018 by Mathieu Pellerin
6 : : email : nirvn dot asia 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 "qgsalgorithmzonalhistogram.h"
19 : : #include "qgsrasteranalysisutils.h"
20 : : #include "qgslogger.h"
21 : :
22 : : ///@cond PRIVATE
23 : :
24 : 0 : QString QgsZonalHistogramAlgorithm::name() const
25 : : {
26 : 0 : return QStringLiteral( "zonalhistogram" );
27 : : }
28 : :
29 : 0 : QString QgsZonalHistogramAlgorithm::displayName() const
30 : : {
31 : 0 : return QObject::tr( "Zonal histogram" );
32 : : }
33 : :
34 : 0 : QStringList QgsZonalHistogramAlgorithm::tags() const
35 : : {
36 : 0 : return QObject::tr( "raster,unique,values,count,area,statistics" ).split( ',' );
37 : 0 : }
38 : :
39 : 0 : QString QgsZonalHistogramAlgorithm::group() const
40 : : {
41 : 0 : return QObject::tr( "Raster analysis" );
42 : : }
43 : :
44 : 0 : QString QgsZonalHistogramAlgorithm::groupId() const
45 : : {
46 : 0 : return QStringLiteral( "rasteranalysis" );
47 : : }
48 : :
49 : 0 : void QgsZonalHistogramAlgorithm::initAlgorithm( const QVariantMap & )
50 : : {
51 : 0 : addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT_RASTER" ),
52 : 0 : QObject::tr( "Raster layer" ) ) );
53 : 0 : addParameter( new QgsProcessingParameterBand( QStringLiteral( "RASTER_BAND" ),
54 : 0 : QObject::tr( "Band number" ), 1, QStringLiteral( "INPUT_RASTER" ) ) );
55 : :
56 : 0 : addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT_VECTOR" ),
57 : 0 : QObject::tr( "Vector layer containing zones" ), QList< int >() << QgsProcessing::TypeVectorPolygon ) );
58 : :
59 : 0 : addParameter( new QgsProcessingParameterString( QStringLiteral( "COLUMN_PREFIX" ), QObject::tr( "Output column prefix" ), QStringLiteral( "HISTO_" ), false, true ) );
60 : :
61 : 0 : addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Output zones" ), QgsProcessing::TypeVectorPolygon ) );
62 : 0 : }
63 : :
64 : 0 : QString QgsZonalHistogramAlgorithm::shortHelpString() const
65 : : {
66 : 0 : return QObject::tr( "This algorithm appends fields representing counts of each unique value from a raster layer contained within zones defined as polygons." );
67 : : }
68 : :
69 : 0 : QgsZonalHistogramAlgorithm *QgsZonalHistogramAlgorithm::createInstance() const
70 : : {
71 : 0 : return new QgsZonalHistogramAlgorithm();
72 : 0 : }
73 : :
74 : 0 : bool QgsZonalHistogramAlgorithm::prepareAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback * )
75 : : {
76 : 0 : QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT_RASTER" ), context );
77 : 0 : if ( !layer )
78 : 0 : throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT_RASTER" ) ) );
79 : :
80 : 0 : mRasterBand = parameterAsInt( parameters, QStringLiteral( "RASTER_BAND" ), context );
81 : 0 : mHasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( mRasterBand );
82 : 0 : mNodataValue = layer->dataProvider()->sourceNoDataValue( mRasterBand );
83 : 0 : mRasterInterface.reset( layer->dataProvider()->clone() );
84 : 0 : mRasterExtent = layer->extent();
85 : 0 : mCrs = layer->crs();
86 : 0 : mCellSizeX = std::abs( layer->rasterUnitsPerPixelX() );
87 : 0 : mCellSizeY = std::abs( layer->rasterUnitsPerPixelX() );
88 : 0 : mNbCellsXProvider = mRasterInterface->xSize();
89 : 0 : mNbCellsYProvider = mRasterInterface->ySize();
90 : :
91 : 0 : return true;
92 : 0 : }
93 : :
94 : 0 : QVariantMap QgsZonalHistogramAlgorithm::processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
95 : : {
96 : :
97 : 0 : std::unique_ptr< QgsFeatureSource > zones( parameterAsSource( parameters, QStringLiteral( "INPUT_VECTOR" ), context ) );
98 : 0 : if ( !zones )
99 : 0 : throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT_VECTOR" ) ) );
100 : :
101 : 0 : long count = zones->featureCount();
102 : 0 : double step = count > 0 ? 100.0 / count : 1;
103 : 0 : long current = 0;
104 : :
105 : 0 : QList< double > uniqueValues;
106 : 0 : QMap< QgsFeatureId, QHash< double, qgssize > > featuresUniqueValues;
107 : :
108 : : // First loop through the zones to build up a list of unique values across all zones to determine sink fields list
109 : 0 : QgsFeatureRequest request;
110 : 0 : request.setNoAttributes();
111 : 0 : if ( zones->sourceCrs() != mCrs )
112 : : {
113 : 0 : request.setDestinationCrs( mCrs, context.transformContext() );
114 : 0 : }
115 : 0 : QgsFeatureIterator it = zones->getFeatures( request );
116 : 0 : QgsFeature f;
117 : 0 : while ( it.nextFeature( f ) )
118 : : {
119 : 0 : if ( feedback->isCanceled() )
120 : : {
121 : 0 : break;
122 : : }
123 : 0 : feedback->setProgress( current * step );
124 : :
125 : 0 : if ( !f.hasGeometry() )
126 : : {
127 : 0 : current++;
128 : 0 : continue;
129 : : }
130 : :
131 : 0 : QgsGeometry featureGeometry = f.geometry();
132 : 0 : QgsRectangle featureRect = featureGeometry.boundingBox().intersect( mRasterExtent );
133 : 0 : if ( featureRect.isEmpty() )
134 : : {
135 : 0 : current++;
136 : 0 : continue;
137 : : }
138 : :
139 : : int nCellsX, nCellsY;
140 : 0 : QgsRectangle rasterBlockExtent;
141 : 0 : QgsRasterAnalysisUtils::cellInfoForBBox( mRasterExtent, featureRect, mCellSizeX, mCellSizeY, nCellsX, nCellsY, mNbCellsXProvider, mNbCellsYProvider, rasterBlockExtent );
142 : :
143 : 0 : QHash< double, qgssize > fUniqueValues;
144 : 0 : QgsRasterAnalysisUtils::statisticsFromMiddlePointTest( mRasterInterface.get(), mRasterBand, featureGeometry, nCellsX, nCellsY, mCellSizeX, mCellSizeY,
145 : 0 : rasterBlockExtent, [ &fUniqueValues]( double value ) { fUniqueValues[value]++; }, false );
146 : :
147 : 0 : if ( fUniqueValues.count() < 1 )
148 : : {
149 : : // The cell resolution is probably larger than the polygon area. We switch to slower precise pixel - polygon intersection in this case
150 : : // TODO: eventually deal with weight if needed
151 : 0 : QgsRasterAnalysisUtils::statisticsFromPreciseIntersection( mRasterInterface.get(), mRasterBand, featureGeometry, nCellsX, nCellsY, mCellSizeX, mCellSizeY,
152 : 0 : rasterBlockExtent, [ &fUniqueValues]( double value, double ) { fUniqueValues[value]++; }, false );
153 : 0 : }
154 : :
155 : 0 : for ( auto it = fUniqueValues.constBegin(); it != fUniqueValues.constEnd(); ++it )
156 : : {
157 : 0 : if ( uniqueValues.indexOf( it.key() ) == -1 )
158 : : {
159 : 0 : uniqueValues << it.key();
160 : 0 : }
161 : 0 : featuresUniqueValues[f.id()][it.key()] += it.value();
162 : 0 : }
163 : :
164 : 0 : current++;
165 : 0 : }
166 : :
167 : 0 : std::sort( uniqueValues.begin(), uniqueValues.end() );
168 : :
169 : 0 : QString fieldPrefix = parameterAsString( parameters, QStringLiteral( "COLUMN_PREFIX" ), context );
170 : 0 : QgsFields newFields;
171 : 0 : for ( auto it = uniqueValues.constBegin(); it != uniqueValues.constEnd(); ++it )
172 : : {
173 : 0 : newFields.append( QgsField( QStringLiteral( "%1%2" ).arg( fieldPrefix, mHasNoDataValue && *it == mNodataValue ? QStringLiteral( "NODATA" ) : QString::number( *it ) ), QVariant::LongLong, QString(), -1, 0 ) );
174 : 0 : }
175 : 0 : QgsFields fields = QgsProcessingUtils::combineFields( zones->fields(), newFields );
176 : :
177 : 0 : QString dest;
178 : 0 : std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, fields,
179 : 0 : zones->wkbType(), zones->sourceCrs() ) );
180 : 0 : if ( !sink )
181 : 0 : throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );
182 : :
183 : 0 : it = zones->getFeatures( QgsFeatureRequest() );
184 : 0 : while ( it.nextFeature( f ) )
185 : : {
186 : 0 : QgsAttributes attributes = f.attributes();
187 : 0 : QHash< double, qgssize > fUniqueValues = featuresUniqueValues.value( f.id() );
188 : 0 : for ( auto it = uniqueValues.constBegin(); it != uniqueValues.constEnd(); ++it )
189 : : {
190 : 0 : attributes += fUniqueValues.value( *it, 0 );
191 : 0 : }
192 : :
193 : 0 : QgsFeature outputFeature;
194 : 0 : outputFeature.setGeometry( f.geometry() );
195 : 0 : outputFeature.setAttributes( attributes );
196 : :
197 : 0 : sink->addFeature( outputFeature, QgsFeatureSink::FastInsert );
198 : 0 : }
199 : :
200 : 0 : QVariantMap outputs;
201 : 0 : outputs.insert( QStringLiteral( "OUTPUT" ), dest );
202 : 0 : return outputs;
203 : 0 : }
204 : :
205 : : ///@endcond
206 : :
207 : :
208 : :
|
