QGIS API Documentation  2.99.0-Master (19b062c)
qgsrastercalculator.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsrastercalculator.cpp - description
3  -----------------------
4  begin : September 28th, 2010
5  copyright : (C) 2010 by Marco Hugentobler
6  email : marco dot hugentobler at sourcepole dot ch
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 "qgsrastercalculator.h"
19 #include "qgsrastercalcnode.h"
20 #include "qgsrasterdataprovider.h"
21 #include "qgsrasterinterface.h"
22 #include "qgsrasterlayer.h"
23 #include "qgsrastermatrix.h"
24 #include "qgsrasterprojector.h"
25 #include "qgsfeedback.h"
26 #include "qgsogrutils.h"
27 
28 #include <QFile>
29 
30 #include <cpl_string.h>
31 #include <gdalwarper.h>
32 
33 QgsRasterCalculator::QgsRasterCalculator( const QString &formulaString, const QString &outputFile, const QString &outputFormat,
34  const QgsRectangle &outputExtent, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry> &rasterEntries )
35  : mFormulaString( formulaString )
36  , mOutputFile( outputFile )
37  , mOutputFormat( outputFormat )
38  , mOutputRectangle( outputExtent )
39  , mNumOutputColumns( nOutputColumns )
40  , mNumOutputRows( nOutputRows )
41  , mRasterEntries( rasterEntries )
42 {
43  //default to first layer's crs
44  mOutputCrs = mRasterEntries.at( 0 ).raster->crs();
45 }
46 
47 QgsRasterCalculator::QgsRasterCalculator( const QString &formulaString, const QString &outputFile, const QString &outputFormat,
48  const QgsRectangle &outputExtent, const QgsCoordinateReferenceSystem &outputCrs, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry> &rasterEntries )
49  : mFormulaString( formulaString )
50  , mOutputFile( outputFile )
51  , mOutputFormat( outputFormat )
52  , mOutputRectangle( outputExtent )
53  , mOutputCrs( outputCrs )
54  , mNumOutputColumns( nOutputColumns )
55  , mNumOutputRows( nOutputRows )
56  , mRasterEntries( rasterEntries )
57 {
58 }
59 
61 {
62  //prepare search string / tree
63  QString errorString;
64  QgsRasterCalcNode *calcNode = QgsRasterCalcNode::parseRasterCalcString( mFormulaString, errorString );
65  if ( !calcNode )
66  {
67  //error
68  return static_cast<int>( ParserError );
69  }
70 
71  QMap< QString, QgsRasterBlock * > inputBlocks;
72  QVector<QgsRasterCalculatorEntry>::const_iterator it = mRasterEntries.constBegin();
73  for ( ; it != mRasterEntries.constEnd(); ++it )
74  {
75  if ( !it->raster ) // no raster layer in entry
76  {
77  delete calcNode;
78  qDeleteAll( inputBlocks );
79  return static_cast< int >( InputLayerError );
80  }
81 
82  QgsRasterBlock *block = nullptr;
83  // if crs transform needed
84  if ( it->raster->crs() != mOutputCrs )
85  {
86  QgsRasterProjector proj;
87  proj.setCrs( it->raster->crs(), mOutputCrs );
88  proj.setInput( it->raster->dataProvider() );
90 
91  block = proj.block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows );
92  }
93  else
94  {
95  block = it->raster->dataProvider()->block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows );
96  }
97  if ( block->isEmpty() )
98  {
99  delete block;
100  delete calcNode;
101  qDeleteAll( inputBlocks );
102  return static_cast<int>( MemoryError );
103  }
104  inputBlocks.insert( it->ref, block );
105  }
106 
107  //open output dataset for writing
108  GDALDriverH outputDriver = openOutputDriver();
109  if ( !outputDriver )
110  {
111  return static_cast< int >( CreateOutputError );
112  }
113 
114  gdal::dataset_unique_ptr outputDataset( openOutputFile( outputDriver ) );
115  GDALSetProjection( outputDataset.get(), mOutputCrs.toWkt().toLocal8Bit().data() );
116  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
117 
118  float outputNodataValue = -FLT_MAX;
119  GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
120 
121  QgsRasterMatrix resultMatrix;
122  resultMatrix.setNodataValue( outputNodataValue );
123 
124  //read / write line by line
125  for ( int i = 0; i < mNumOutputRows; ++i )
126  {
127  if ( feedback )
128  {
129  feedback->setProgress( 100.0 * static_cast< double >( i ) / mNumOutputRows );
130  }
131 
132  if ( feedback && feedback->isCanceled() )
133  {
134  break;
135  }
136 
137  if ( calcNode->calculate( inputBlocks, resultMatrix, i ) )
138  {
139  bool resultIsNumber = resultMatrix.isNumber();
140  float *calcData = new float[mNumOutputColumns];
141 
142  for ( int j = 0; j < mNumOutputColumns; ++j )
143  {
144  calcData[j] = ( float )( resultIsNumber ? resultMatrix.number() : resultMatrix.data()[j] );
145  }
146 
147  //write scanline to the dataset
148  if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
149  {
150  QgsDebugMsg( "RasterIO error!" );
151  }
152 
153  delete[] calcData;
154  }
155 
156  }
157 
158  if ( feedback )
159  {
160  feedback->setProgress( 100.0 );
161  }
162 
163  //close datasets and release memory
164  delete calcNode;
165  qDeleteAll( inputBlocks );
166  inputBlocks.clear();
167 
168  if ( feedback && feedback->isCanceled() )
169  {
170  //delete the dataset without closing (because it is faster)
171  gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
172  return static_cast< int >( Canceled );
173  }
174  return static_cast< int >( Success );
175 }
176 
177 GDALDriverH QgsRasterCalculator::openOutputDriver()
178 {
179  char **driverMetadata = nullptr;
180 
181  //open driver
182  GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
183 
184  if ( !outputDriver )
185  {
186  return outputDriver; //return nullptr, driver does not exist
187  }
188 
189  driverMetadata = GDALGetMetadata( outputDriver, nullptr );
190  if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
191  {
192  return nullptr; //driver exist, but it does not support the create operation
193  }
194 
195  return outputDriver;
196 }
197 
198 gdal::dataset_unique_ptr QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver )
199 {
200  //open output file
201  char **papszOptions = nullptr;
202  gdal::dataset_unique_ptr outputDataset( GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions ) );
203  if ( !outputDataset )
204  {
205  return nullptr;
206  }
207 
208  //assign georef information
209  double geotransform[6];
210  outputGeoTransform( geotransform );
211  GDALSetGeoTransform( outputDataset.get(), geotransform );
212 
213  return outputDataset;
214 }
215 
216 void QgsRasterCalculator::outputGeoTransform( double *transform ) const
217 {
218  transform[0] = mOutputRectangle.xMinimum();
219  transform[1] = mOutputRectangle.width() / mNumOutputColumns;
220  transform[2] = 0;
221  transform[3] = mOutputRectangle.yMaximum();
222  transform[4] = 0;
223  transform[5] = -mOutputRectangle.height() / mNumOutputRows;
224 }
A rectangle specified with double values.
Definition: qgsrectangle.h:39
bool isNumber() const
Returns true if matrix is 1x1 (=scalar number)
#define QgsDebugMsg(str)
Definition: qgslogger.h:37
Error allocating memory for result.
void setNodataValue(double d)
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
QgsRasterCalculator(const QString &formulaString, const QString &outputFile, const QString &outputFormat, const QgsRectangle &outputExtent, int nOutputColumns, int nOutputRows, const QVector< QgsRasterCalculatorEntry > &rasterEntries)
QgsRasterCalculator constructor.
int processCalculation(QgsFeedback *feedback=nullptr)
Starts the calculation and writes a new raster.
User canceled calculation.
void setCrs(const QgsCoordinateReferenceSystem &srcCRS, const QgsCoordinateReferenceSystem &destCRS, int srcDatumTransform=-1, int destDatumTransform=-1)
set source and destination CRS
Base class for feedback objects to be used for cancelation of something running in a worker thread...
Definition: qgsfeedback.h:44
bool calculate(QMap< QString, QgsRasterBlock * > &rasterData, QgsRasterMatrix &result, int row=-1) const
Calculates result of raster calculation (might be real matrix or single number).
void CORE_EXPORT fast_delete_and_close(dataset_unique_ptr &dataset, GDALDriverH driver, const QString &path)
Performs a fast close of an unwanted GDAL dataset handle by deleting the underlying data store...
Definition: qgsogrutils.cpp:60
Raster data container.
Calculation successful.
double width() const
Returns the width of the rectangle.
Definition: qgsrectangle.h:138
bool isEmpty() const
Returns true if block is empty, i.e.
Error creating output data file.
std::unique_ptr< void, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
Definition: qgsogrutils.h:134
QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback=nullptr) override
Read block of data using given extent and size.
void setPrecision(Precision precision)
QgsRasterProjector implements approximate projection support for it calculates grid of points in sour...
virtual bool setInput(QgsRasterInterface *input)
Set input.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
double number() const
This class represents a coordinate reference system (CRS).
QString toWkt() const
Returns a WKT representation of this CRS.
double * data()
Returns data array (but not ownership)
double xMinimum() const
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:116
static QgsRasterCalcNode * parseRasterCalcString(const QString &str, QString &parserErrorMsg)
double yMaximum() const
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:121
Exact, precise but slow.
double height() const
Returns the height of the rectangle.
Definition: qgsrectangle.h:145