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