Branch data Line data Source code
1 : : /***************************************************************************
2 : : qgsalgorithmrasterlayeruniquevalues.cpp
3 : : ---------------------
4 : : begin : January 2019
5 : : copyright : (C) 2019 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 "qgsalgorithmrastersurfacevolume.h"
19 : : #include "qgsstringutils.h"
20 : : #include <QTextStream>
21 : :
22 : : ///@cond PRIVATE
23 : :
24 : 0 : QString QgsRasterSurfaceVolumeAlgorithm::name() const
25 : : {
26 : 0 : return QStringLiteral( "rastersurfacevolume" );
27 : : }
28 : :
29 : 0 : QString QgsRasterSurfaceVolumeAlgorithm::displayName() const
30 : : {
31 : 0 : return QObject::tr( "Raster surface volume" );
32 : : }
33 : :
34 : 0 : QStringList QgsRasterSurfaceVolumeAlgorithm::tags() const
35 : : {
36 : 0 : return QObject::tr( "sum,volume,area,height,terrain,dem,elevation" ).split( ',' );
37 : 0 : }
38 : :
39 : 0 : QString QgsRasterSurfaceVolumeAlgorithm::group() const
40 : : {
41 : 0 : return QObject::tr( "Raster analysis" );
42 : : }
43 : :
44 : 0 : QString QgsRasterSurfaceVolumeAlgorithm::groupId() const
45 : : {
46 : 0 : return QStringLiteral( "rasteranalysis" );
47 : : }
48 : :
49 : 0 : void QgsRasterSurfaceVolumeAlgorithm::initAlgorithm( const QVariantMap & )
50 : : {
51 : 0 : addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT" ),
52 : 0 : QObject::tr( "Input layer" ) ) );
53 : 0 : addParameter( new QgsProcessingParameterBand( QStringLiteral( "BAND" ),
54 : 0 : QObject::tr( "Band number" ), 1, QStringLiteral( "INPUT" ) ) );
55 : 0 : addParameter( new QgsProcessingParameterNumber( QStringLiteral( "LEVEL" ),
56 : 0 : QObject::tr( "Base level" ), QgsProcessingParameterNumber::Double, 0 ) );
57 : 0 : addParameter( new QgsProcessingParameterEnum( QStringLiteral( "METHOD" ),
58 : 0 : QObject::tr( "Method" ), QStringList()
59 : 0 : << QObject::tr( "Count Only Above Base Level" )
60 : 0 : << QObject::tr( "Count Only Below Base Level" )
61 : 0 : << QObject::tr( "Subtract Volumes Below Base Level" )
62 : 0 : << QObject::tr( "Add Volumes Below Base Level" ) ) );
63 : :
64 : 0 : addParameter( new QgsProcessingParameterFileDestination( QStringLiteral( "OUTPUT_HTML_FILE" ),
65 : 0 : QObject::tr( "Surface volume report" ), QObject::tr( "HTML files (*.html)" ), QVariant(), true ) );
66 : 0 : addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT_TABLE" ),
67 : 0 : QObject::tr( "Surface volume table" ), QgsProcessing::TypeVector, QVariant(), true, false ) );
68 : :
69 : 0 : addOutput( new QgsProcessingOutputNumber( QStringLiteral( "VOLUME" ), QObject::tr( "Volume" ) ) );
70 : 0 : addOutput( new QgsProcessingOutputNumber( QStringLiteral( "PIXEL_COUNT" ), QObject::tr( "Pixel count" ) ) );
71 : 0 : addOutput( new QgsProcessingOutputNumber( QStringLiteral( "AREA" ), QObject::tr( "Area" ) ) );
72 : 0 : }
73 : :
74 : 0 : QString QgsRasterSurfaceVolumeAlgorithm::shortHelpString() const
75 : : {
76 : 0 : return QObject::tr( "This algorithm calculates the volume under a raster grid's surface.\n\n"
77 : : "Several methods of volume calculation are available, which control whether "
78 : : "only values above or below the specified base level are considered, or "
79 : : "whether volumes below the base level should be added or subtracted from the total volume.\n\n"
80 : : "The algorithm outputs the calculated volume, the total area, and the total number of pixels analysed. "
81 : : "If the 'Count Only Above Base Level' or 'Count Only Below Base Level' methods are used, "
82 : : "then the calculated area and pixel count only includes pixels which are above or below the "
83 : : "specified base level respectively.\n\n"
84 : : "Units of the calculated volume are dependent on the coordinate reference system of "
85 : : "the input raster file. For a CRS in meters, with a DEM height in meters, the calculated "
86 : : "value will be in meters³. If instead the input raster is in a geographic coordinate system "
87 : : "(e.g. latitude/longitude values), then the result will be in degrees² × meters, and an "
88 : : "appropriate scaling factor will need to be applied in order to convert to meters³." );
89 : : }
90 : :
91 : 0 : QString QgsRasterSurfaceVolumeAlgorithm::shortDescription() const
92 : : {
93 : 0 : return QObject::tr( "Calculates the volume under a raster grid's surface." );
94 : : }
95 : :
96 : 0 : QgsRasterSurfaceVolumeAlgorithm *QgsRasterSurfaceVolumeAlgorithm::createInstance() const
97 : : {
98 : 0 : return new QgsRasterSurfaceVolumeAlgorithm();
99 : 0 : }
100 : :
101 : 0 : bool QgsRasterSurfaceVolumeAlgorithm::prepareAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback * )
102 : : {
103 : 0 : QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT" ), context );
104 : 0 : int band = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
105 : :
106 : 0 : if ( !layer )
107 : 0 : throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT" ) ) );
108 : :
109 : 0 : mBand = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
110 : 0 : if ( mBand < 1 || mBand > layer->bandCount() )
111 : 0 : throw QgsProcessingException( QObject::tr( "Invalid band number for BAND (%1): Valid values for input raster are 1 to %2" ).arg( mBand )
112 : 0 : .arg( layer->bandCount() ) );
113 : :
114 : 0 : mInterface.reset( layer->dataProvider()->clone() );
115 : 0 : mHasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( band );
116 : 0 : mLayerWidth = layer->width();
117 : 0 : mLayerHeight = layer->height();
118 : 0 : mExtent = layer->extent();
119 : 0 : mCrs = layer->crs();
120 : 0 : mRasterUnitsPerPixelX = layer->rasterUnitsPerPixelX();
121 : 0 : mRasterUnitsPerPixelY = layer->rasterUnitsPerPixelY();
122 : 0 : mSource = layer->source();
123 : :
124 : 0 : mLevel = parameterAsDouble( parameters, QStringLiteral( "LEVEL" ), context );
125 : 0 : mMethod = static_cast< Method >( parameterAsEnum( parameters, QStringLiteral( "METHOD" ), context ) );
126 : 0 : return true;
127 : 0 : }
128 : :
129 : 0 : QVariantMap QgsRasterSurfaceVolumeAlgorithm::processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
130 : : {
131 : 0 : QString outputFile = parameterAsFileOutput( parameters, QStringLiteral( "OUTPUT_HTML_FILE" ), context );
132 : 0 : QString areaUnit = QgsUnitTypes::toAbbreviatedString( QgsUnitTypes::distanceToAreaUnit( mCrs.mapUnits() ) );
133 : :
134 : 0 : QString tableDest;
135 : 0 : std::unique_ptr< QgsFeatureSink > sink;
136 : 0 : if ( parameters.contains( QStringLiteral( "OUTPUT_TABLE" ) ) && parameters.value( QStringLiteral( "OUTPUT_TABLE" ) ).isValid() )
137 : : {
138 : 0 : QgsFields outFields;
139 : 0 : outFields.append( QgsField( QStringLiteral( "volume" ), QVariant::Double, QString(), 20, 8 ) );
140 : 0 : outFields.append( QgsField( areaUnit.replace( QStringLiteral( "²" ), QStringLiteral( "2" ) ), QVariant::Double, QString(), 20, 8 ) );
141 : 0 : outFields.append( QgsField( QStringLiteral( "pixel_count" ), QVariant::LongLong ) );
142 : 0 : sink.reset( parameterAsSink( parameters, QStringLiteral( "OUTPUT_TABLE" ), context, tableDest, outFields, QgsWkbTypes::NoGeometry, QgsCoordinateReferenceSystem() ) );
143 : 0 : if ( !sink )
144 : 0 : throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT_TABLE" ) ) );
145 : 0 : }
146 : :
147 : 0 : double volume = 0;
148 : 0 : long long count = 0;
149 : :
150 : 0 : int maxWidth = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_WIDTH;
151 : 0 : int maxHeight = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_HEIGHT;
152 : 0 : int nbBlocksWidth = static_cast< int >( std::ceil( 1.0 * mLayerWidth / maxWidth ) );
153 : 0 : int nbBlocksHeight = static_cast< int >( std::ceil( 1.0 * mLayerHeight / maxHeight ) );
154 : 0 : int nbBlocks = nbBlocksWidth * nbBlocksHeight;
155 : :
156 : 0 : QgsRasterIterator iter( mInterface.get() );
157 : 0 : iter.startRasterRead( mBand, mLayerWidth, mLayerHeight, mExtent );
158 : :
159 : 0 : int iterLeft = 0;
160 : 0 : int iterTop = 0;
161 : 0 : int iterCols = 0;
162 : 0 : int iterRows = 0;
163 : 0 : std::unique_ptr< QgsRasterBlock > rasterBlock;
164 : 0 : while ( iter.readNextRasterPart( mBand, iterCols, iterRows, rasterBlock, iterLeft, iterTop ) )
165 : : {
166 : 0 : feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
167 : 0 : for ( int row = 0; row < iterRows; row++ )
168 : : {
169 : 0 : if ( feedback->isCanceled() )
170 : 0 : break;
171 : 0 : for ( int column = 0; column < iterCols; column++ )
172 : : {
173 : 0 : if ( mHasNoDataValue && rasterBlock->isNoData( row, column ) )
174 : : {
175 : 0 : continue;
176 : : }
177 : :
178 : 0 : const double z = rasterBlock->value( row, column ) - mLevel;
179 : :
180 : 0 : switch ( mMethod )
181 : : {
182 : : case CountOnlyAboveBaseLevel:
183 : 0 : if ( z > 0.0 )
184 : : {
185 : 0 : volume += z;
186 : 0 : count++;
187 : 0 : }
188 : 0 : continue;
189 : :
190 : : case CountOnlyBelowBaseLevel:
191 : 0 : if ( z < 0.0 )
192 : : {
193 : 0 : volume += z;
194 : 0 : count++;
195 : 0 : }
196 : 0 : continue;
197 : :
198 : : case SubtractVolumesBelowBaseLevel:
199 : 0 : volume += z;
200 : 0 : count++;
201 : 0 : continue;
202 : :
203 : : case AddVolumesBelowBaseLevel:
204 : 0 : volume += std::fabs( z );
205 : 0 : count++;
206 : 0 : continue;
207 : : }
208 : 0 : }
209 : 0 : }
210 : : }
211 : :
212 : 0 : QVariantMap outputs;
213 : 0 : double pixelArea = mRasterUnitsPerPixelX * mRasterUnitsPerPixelY;
214 : 0 : double area = count * pixelArea;
215 : 0 : volume *= pixelArea;
216 : 0 : if ( !outputFile.isEmpty() )
217 : : {
218 : 0 : QFile file( outputFile );
219 : 0 : if ( file.open( QIODevice::WriteOnly | QIODevice::Text ) )
220 : : {
221 : 0 : const QString encodedAreaUnit = QgsStringUtils::ampersandEncode( areaUnit );
222 : :
223 : 0 : QTextStream out( &file );
224 : 0 : out << QStringLiteral( "<html><head><meta http-equiv=\"Content-Type\" content=\"text/html;charset=utf-8\"/></head><body>\n" );
225 : 0 : out << QStringLiteral( "<p>%1: %2 (%3 %4)</p>\n" ).arg( QObject::tr( "Analyzed file" ), mSource, QObject::tr( "band" ) ).arg( mBand );
226 : 0 : out << QObject::tr( "<p>%1: %2</p>\n" ).arg( QObject::tr( "Volume" ), QString::number( volume, 'g', 16 ) );
227 : 0 : out << QObject::tr( "<p>%1: %2</p>\n" ).arg( QObject::tr( "Pixel count" ) ).arg( count );
228 : 0 : out << QObject::tr( "<p>%1: %2 %3</p>\n" ).arg( QObject::tr( "Area" ), QString::number( area, 'g', 16 ), encodedAreaUnit );
229 : 0 : out << QStringLiteral( "</body></html>" );
230 : 0 : outputs.insert( QStringLiteral( "OUTPUT_HTML_FILE" ), outputFile );
231 : 0 : }
232 : 0 : }
233 : :
234 : 0 : if ( sink )
235 : : {
236 : 0 : QgsFeature f;
237 : 0 : f.setAttributes( QgsAttributes() << volume << area << count );
238 : 0 : sink->addFeature( f, QgsFeatureSink::FastInsert );
239 : 0 : outputs.insert( QStringLiteral( "OUTPUT_TABLE" ), tableDest );
240 : 0 : }
241 : 0 : outputs.insert( QStringLiteral( "VOLUME" ), volume );
242 : 0 : outputs.insert( QStringLiteral( "AREA" ), area );
243 : 0 : outputs.insert( QStringLiteral( "PIXEL_COUNT" ), count );
244 : 0 : return outputs;
245 : 0 : }
246 : :
247 : :
248 : : ///@endcond
249 : :
250 : :
251 : :
|
