QGIS API Documentation  3.14.0-Pi (9f7028fd23)
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 
22 
23 QString QgsRasterSurfaceVolumeAlgorithm::name() const
24 {
25  return QStringLiteral( "rastersurfacevolume" );
26 }
27 
28 QString QgsRasterSurfaceVolumeAlgorithm::displayName() const
29 {
30  return QObject::tr( "Raster surface volume" );
31 }
32 
33 QStringList QgsRasterSurfaceVolumeAlgorithm::tags() const
34 {
35  return QObject::tr( "sum,volume,area,height,terrain,dem,elevation" ).split( ',' );
36 }
37 
38 QString QgsRasterSurfaceVolumeAlgorithm::group() const
39 {
40  return QObject::tr( "Raster analysis" );
41 }
42 
43 QString QgsRasterSurfaceVolumeAlgorithm::groupId() const
44 {
45  return QStringLiteral( "rasteranalysis" );
46 }
47 
48 void QgsRasterSurfaceVolumeAlgorithm::initAlgorithm( const QVariantMap & )
49 {
50  addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT" ),
51  QObject::tr( "Input layer" ) ) );
52  addParameter( new QgsProcessingParameterBand( QStringLiteral( "BAND" ),
53  QObject::tr( "Band number" ), 1, QStringLiteral( "INPUT" ) ) );
54  addParameter( new QgsProcessingParameterNumber( QStringLiteral( "LEVEL" ),
55  QObject::tr( "Base level" ), QgsProcessingParameterNumber::Double, 0 ) );
56  addParameter( new QgsProcessingParameterEnum( QStringLiteral( "METHOD" ),
57  QObject::tr( "Method" ), QStringList()
58  << QObject::tr( "Count Only Above Base Level" )
59  << QObject::tr( "Count Only Below Base Level" )
60  << QObject::tr( "Subtract Volumes Below Base Level" )
61  << QObject::tr( "Add Volumes Below Base Level" ) ) );
62 
63  addParameter( new QgsProcessingParameterFileDestination( QStringLiteral( "OUTPUT_HTML_FILE" ),
64  QObject::tr( "Surface volume report" ), QObject::tr( "HTML files (*.html)" ), QVariant(), true ) );
65  addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT_TABLE" ),
66  QObject::tr( "Surface volume table" ), QgsProcessing::TypeVector, QVariant(), true, false ) );
67 
68  addOutput( new QgsProcessingOutputNumber( QStringLiteral( "VOLUME" ), QObject::tr( "Volume" ) ) );
69  addOutput( new QgsProcessingOutputNumber( QStringLiteral( "PIXEL_COUNT" ), QObject::tr( "Pixel count" ) ) );
70  addOutput( new QgsProcessingOutputNumber( QStringLiteral( "AREA" ), QObject::tr( "Area" ) ) );
71 }
72 
73 QString QgsRasterSurfaceVolumeAlgorithm::shortHelpString() const
74 {
75  return QObject::tr( "This algorithm calculates the volume under a raster grid's surface.\n\n"
76  "Several methods of volume calculation are available, which control whether "
77  "only values above or below the specified base level are considered, or "
78  "whether volumes below the base level should be added or subtracted from the total volume.\n\n"
79  "The algorithm outputs the calculated volume, the total area, and the total number of pixels analysed. "
80  "If the 'Count Only Above Base Level' or 'Count Only Below Base Level' methods are used, "
81  "then the calculated area and pixel count only includes pixels which are above or below the "
82  "specified base level respectively.\n\n"
83  "Units of the calculated volume are dependent on the coordinate reference system of "
84  "the input raster file. For a CRS in meters, with a DEM height in meters, the calculated "
85  "value will be in meters³. If instead the input raster is in a geographic coordinate system "
86  "(e.g. latitude/longitude values), then the result will be in degrees² × meters, and an "
87  "appropriate scaling factor will need to be applied in order to convert to meters³." );
88 }
89 
90 QString QgsRasterSurfaceVolumeAlgorithm::shortDescription() const
91 {
92  return QObject::tr( "Calculates the volume under a raster grid's surface." );
93 }
94 
95 QgsRasterSurfaceVolumeAlgorithm *QgsRasterSurfaceVolumeAlgorithm::createInstance() const
96 {
97  return new QgsRasterSurfaceVolumeAlgorithm();
98 }
99 
100 bool QgsRasterSurfaceVolumeAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
101 {
102  QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT" ), context );
103  int band = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
104 
105  if ( !layer )
106  throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT" ) ) );
107 
108  mBand = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
109  if ( mBand < 1 || mBand > layer->bandCount() )
110  throw QgsProcessingException( QObject::tr( "Invalid band number for BAND (%1): Valid values for input raster are 1 to %2" ).arg( mBand )
111  .arg( layer->bandCount() ) );
112 
113  mInterface.reset( layer->dataProvider()->clone() );
114  mHasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( band );
115  mLayerWidth = layer->width();
116  mLayerHeight = layer->height();
117  mExtent = layer->extent();
118  mCrs = layer->crs();
119  mRasterUnitsPerPixelX = layer->rasterUnitsPerPixelX();
120  mRasterUnitsPerPixelY = layer->rasterUnitsPerPixelY();
121  mSource = layer->source();
122 
123  mLevel = parameterAsDouble( parameters, QStringLiteral( "LEVEL" ), context );
124  mMethod = static_cast< Method >( parameterAsEnum( parameters, QStringLiteral( "METHOD" ), context ) );
125  return true;
126 }
127 
128 QVariantMap QgsRasterSurfaceVolumeAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
129 {
130  QString outputFile = parameterAsFileOutput( parameters, QStringLiteral( "OUTPUT_HTML_FILE" ), context );
131  QString areaUnit = QgsUnitTypes::toAbbreviatedString( QgsUnitTypes::distanceToAreaUnit( mCrs.mapUnits() ) );
132 
133  QString tableDest;
134  std::unique_ptr< QgsFeatureSink > sink;
135  if ( parameters.contains( QStringLiteral( "OUTPUT_TABLE" ) ) && parameters.value( QStringLiteral( "OUTPUT_TABLE" ) ).isValid() )
136  {
137  QgsFields outFields;
138  outFields.append( QgsField( QStringLiteral( "volume" ), QVariant::Double, QString(), 20, 8 ) );
139  outFields.append( QgsField( areaUnit.replace( QStringLiteral( "²" ), QStringLiteral( "2" ) ), QVariant::Double, QString(), 20, 8 ) );
140  outFields.append( QgsField( QStringLiteral( "pixel_count" ), QVariant::LongLong ) );
141  sink.reset( parameterAsSink( parameters, QStringLiteral( "OUTPUT_TABLE" ), context, tableDest, outFields, QgsWkbTypes::NoGeometry, QgsCoordinateReferenceSystem() ) );
142  if ( !sink )
143  throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT_TABLE" ) ) );
144  }
145 
146  double volume = 0;
147  long long count = 0;
148 
151  int nbBlocksWidth = static_cast< int >( std::ceil( 1.0 * mLayerWidth / maxWidth ) );
152  int nbBlocksHeight = static_cast< int >( std::ceil( 1.0 * mLayerHeight / maxHeight ) );
153  int nbBlocks = nbBlocksWidth * nbBlocksHeight;
154 
155  QgsRasterIterator iter( mInterface.get() );
156  iter.startRasterRead( mBand, mLayerWidth, mLayerHeight, mExtent );
157 
158  int iterLeft = 0;
159  int iterTop = 0;
160  int iterCols = 0;
161  int iterRows = 0;
162  std::unique_ptr< QgsRasterBlock > rasterBlock;
163  while ( iter.readNextRasterPart( mBand, iterCols, iterRows, rasterBlock, iterLeft, iterTop ) )
164  {
165  feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
166  for ( int row = 0; row < iterRows; row++ )
167  {
168  if ( feedback->isCanceled() )
169  break;
170  for ( int column = 0; column < iterCols; column++ )
171  {
172  if ( mHasNoDataValue && rasterBlock->isNoData( row, column ) )
173  {
174  continue;
175  }
176 
177  const double z = rasterBlock->value( row, column ) - mLevel;
178 
179  switch ( mMethod )
180  {
181  case CountOnlyAboveBaseLevel:
182  if ( z > 0.0 )
183  {
184  volume += z;
185  count++;
186  }
187  continue;
188 
189  case CountOnlyBelowBaseLevel:
190  if ( z < 0.0 )
191  {
192  volume += z;
193  count++;
194  }
195  continue;
196 
197  case SubtractVolumesBelowBaseLevel:
198  volume += z;
199  count++;
200  continue;
201 
202  case AddVolumesBelowBaseLevel:
203  volume += std::fabs( z );
204  count++;
205  continue;
206  }
207  }
208  }
209  }
210 
211  QVariantMap outputs;
212  double pixelArea = mRasterUnitsPerPixelX * mRasterUnitsPerPixelY;
213  double area = count * pixelArea;
214  volume *= pixelArea;
215  if ( !outputFile.isEmpty() )
216  {
217  QFile file( outputFile );
218  if ( file.open( QIODevice::WriteOnly | QIODevice::Text ) )
219  {
220  const QString encodedAreaUnit = QgsStringUtils::ampersandEncode( areaUnit );
221 
222  QTextStream out( &file );
223  out << QStringLiteral( "<html><head><meta http-equiv=\"Content-Type\" content=\"text/html;charset=utf-8\"/></head><body>\n" );
224  out << QStringLiteral( "<p>%1: %2 (%3 %4)</p>\n" ).arg( QObject::tr( "Analyzed file" ), mSource, QObject::tr( "band" ) ).arg( mBand );
225  out << QObject::tr( "<p>%1: %2</p>\n" ).arg( QObject::tr( "Volume" ), QString::number( volume, 'g', 16 ) );
226  out << QObject::tr( "<p>%1: %2</p>\n" ).arg( QObject::tr( "Pixel count" ) ).arg( count );
227  out << QObject::tr( "<p>%1: %2 %3</p>\n" ).arg( QObject::tr( "Area" ), QString::number( area, 'g', 16 ), encodedAreaUnit );
228  out << QStringLiteral( "</body></html>" );
229  outputs.insert( QStringLiteral( "OUTPUT_HTML_FILE" ), outputFile );
230  }
231  }
232 
233  if ( sink )
234  {
235  QgsFeature f;
236  f.setAttributes( QgsAttributes() << volume << area << count );
237  sink->addFeature( f, QgsFeatureSink::FastInsert );
238  outputs.insert( QStringLiteral( "OUTPUT_TABLE" ), tableDest );
239  }
240  outputs.insert( QStringLiteral( "VOLUME" ), volume );
241  outputs.insert( QStringLiteral( "AREA" ), area );
242  outputs.insert( QStringLiteral( "PIXEL_COUNT" ), count );
243  return outputs;
244 }
245 
246 
248 
249 
250 
QgsMapLayer::crs
QgsCoordinateReferenceSystem crs
Definition: qgsmaplayer.h:88
QgsFeedback::setProgress
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:75
QgsRasterLayer::bandCount
int bandCount() const
Returns the number of bands in this layer.
Definition: qgsrasterlayer.cpp:216
QgsProcessingParameterNumber::Double
@ Double
Double/float values.
Definition: qgsprocessingparameters.h:1846
QgsProcessingParameterNumber
Definition: qgsprocessingparameters.h:1838
QgsProcessingFeedback
Definition: qgsprocessingfeedback.h:37
qgsstringutils.h
QgsFields
Definition: qgsfields.h:44
QgsUnitTypes::distanceToAreaUnit
static Q_INVOKABLE QgsUnitTypes::AreaUnit distanceToAreaUnit(QgsUnitTypes::DistanceUnit distanceUnit)
Converts a distance unit to its corresponding area unit, e.g., meters to square meters.
Definition: qgsunittypes.cpp:1176
QgsUnitTypes::toAbbreviatedString
static Q_INVOKABLE QString toAbbreviatedString(QgsUnitTypes::DistanceUnit unit)
Returns a translated abbreviation representing a distance unit.
Definition: qgsunittypes.cpp:269
QgsStringUtils::ampersandEncode
static QString ampersandEncode(const QString &string)
Makes a raw string safe for inclusion as a HTML/XML string literal.
Definition: qgsstringutils.cpp:112
QgsFields::append
bool append(const QgsField &field, FieldOrigin origin=OriginProvider, int originIndex=-1)
Appends a field. The field must have unique name, otherwise it is rejected (returns false)
Definition: qgsfields.cpp:59
qgsalgorithmrastersurfacevolume.h
QgsProcessingOutputNumber
Definition: qgsprocessingoutputs.h:294
QgsRasterDataProvider::sourceHasNoDataValue
virtual bool sourceHasNoDataValue(int bandNo) const
Returns true if source band has no data value.
Definition: qgsrasterdataprovider.h:243
QgsProcessingParameterFeatureSink
Definition: qgsprocessingparameters.h:2773
QgsRasterLayer::width
int width() const
Returns the width of the (unclipped) raster.
Definition: qgsrasterlayer.cpp:2373
QgsProcessing::TypeVector
@ TypeVector
Tables (i.e. vector layers with or without geometry). When used for a sink this indicates the sink ha...
Definition: qgsprocessing.h:53
QgsProcessingContext
Definition: qgsprocessingcontext.h:43
QgsProcessingParameterFileDestination
Definition: qgsprocessingparameters.h:3005
QgsRasterLayer::height
int height() const
Returns the height of the (unclipped) raster.
Definition: qgsrasterlayer.cpp:2379
QgsMapLayer::extent
virtual QgsRectangle extent() const
Returns the extent of the layer.
Definition: qgsmaplayer.cpp:197
QgsRasterIterator::DEFAULT_MAXIMUM_TILE_HEIGHT
static const int DEFAULT_MAXIMUM_TILE_HEIGHT
Default maximum tile height.
Definition: qgsrasteriterator.h:151
QgsProcessingParameterRasterLayer
Definition: qgsprocessingparameters.h:2101
QgsRasterLayer
Definition: qgsrasterlayer.h:72
QgsRasterIterator
Definition: qgsrasteriterator.h:34
QgsCoordinateReferenceSystem
Definition: qgscoordinatereferencesystem.h:206
QgsRasterLayer::rasterUnitsPerPixelY
double rasterUnitsPerPixelY() const
Returns the number of raster units per each raster pixel in Y axis.
Definition: qgsrasterlayer.cpp:566
QgsFeedback::isCanceled
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:66
QgsMapLayer::source
QString source() const
Returns the source for the layer.
Definition: qgsmaplayer.cpp:192
QgsWkbTypes::NoGeometry
@ NoGeometry
Definition: qgswkbtypes.h:84
QgsRasterDataProvider::clone
QgsRasterInterface * clone() const override=0
Clone itself, create deep copy.
QgsRasterIterator::DEFAULT_MAXIMUM_TILE_WIDTH
static const int DEFAULT_MAXIMUM_TILE_WIDTH
Default maximum tile width.
Definition: qgsrasteriterator.h:148
QgsProcessingParameterBand
Definition: qgsprocessingparameters.h:3103
QgsAttributes
Definition: qgsattributes.h:57
QgsFeature
Definition: qgsfeature.h:55
QgsProcessingParameterEnum
Definition: qgsprocessingparameters.h:2134
QgsFeature::setAttributes
void setAttributes(const QgsAttributes &attrs)
Sets the feature's attributes.
Definition: qgsfeature.cpp:127
QgsProcessingException
Definition: qgsexception.h:82
QgsRasterLayer::dataProvider
QgsRasterDataProvider * dataProvider() override
Returns the source data provider.
Definition: qgsrasterlayer.cpp:233
QgsFeatureSink::FastInsert
@ FastInsert
Use faster inserts, at the cost of updating the passed features to reflect changes made at the provid...
Definition: qgsfeaturesink.h:70
QgsRasterLayer::rasterUnitsPerPixelX
double rasterUnitsPerPixelX() const
Returns the number of raster units per each raster pixel in X axis.
Definition: qgsrasterlayer.cpp:550
QgsField
Definition: qgsfield.h:49