QGIS API Documentation  3.1.0-Master (f35745f)
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  QgsRasterBlockFeedback *rasterBlockFeedback = new QgsRasterBlockFeedback();
92  QObject::connect( feedback, &QgsFeedback::canceled, rasterBlockFeedback, &QgsRasterBlockFeedback::cancel );
93  block = proj.block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows, rasterBlockFeedback );
94  if ( rasterBlockFeedback->isCanceled() )
95  {
96  delete block;
97  delete calcNode;
98  qDeleteAll( inputBlocks );
99  return static_cast< int >( Canceled );
100  }
101  }
102  else
103  {
104  block = it->raster->dataProvider()->block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows );
105  }
106  if ( block->isEmpty() )
107  {
108  delete block;
109  delete calcNode;
110  qDeleteAll( inputBlocks );
111  return static_cast<int>( MemoryError );
112  }
113  inputBlocks.insert( it->ref, block );
114  }
115 
116  //open output dataset for writing
117  GDALDriverH outputDriver = openOutputDriver();
118  if ( !outputDriver )
119  {
120  return static_cast< int >( CreateOutputError );
121  }
122 
123  gdal::dataset_unique_ptr outputDataset( openOutputFile( outputDriver ) );
124  if ( !outputDataset )
125  {
126  return static_cast< int >( CreateOutputError );
127  }
128 
129  GDALSetProjection( outputDataset.get(), mOutputCrs.toWkt().toLocal8Bit().data() );
130  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
131 
132  float outputNodataValue = -FLT_MAX;
133  GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
134 
135  QgsRasterMatrix resultMatrix;
136  resultMatrix.setNodataValue( outputNodataValue );
137 
138  //read / write line by line
139  for ( int i = 0; i < mNumOutputRows; ++i )
140  {
141  if ( feedback )
142  {
143  feedback->setProgress( 100.0 * static_cast< double >( i ) / mNumOutputRows );
144  }
145 
146  if ( feedback && feedback->isCanceled() )
147  {
148  break;
149  }
150 
151  if ( calcNode->calculate( inputBlocks, resultMatrix, i ) )
152  {
153  bool resultIsNumber = resultMatrix.isNumber();
154  float *calcData = new float[mNumOutputColumns];
155 
156  for ( int j = 0; j < mNumOutputColumns; ++j )
157  {
158  calcData[j] = ( float )( resultIsNumber ? resultMatrix.number() : resultMatrix.data()[j] );
159  }
160 
161  //write scanline to the dataset
162  if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
163  {
164  QgsDebugMsg( "RasterIO error!" );
165  }
166 
167  delete[] calcData;
168  }
169 
170  }
171 
172  if ( feedback )
173  {
174  feedback->setProgress( 100.0 );
175  }
176 
177  //close datasets and release memory
178  delete calcNode;
179  qDeleteAll( inputBlocks );
180  inputBlocks.clear();
181 
182  if ( feedback && feedback->isCanceled() )
183  {
184  //delete the dataset without closing (because it is faster)
185  gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
186  return static_cast< int >( Canceled );
187  }
188  return static_cast< int >( Success );
189 }
190 
191 GDALDriverH QgsRasterCalculator::openOutputDriver()
192 {
193  char **driverMetadata = nullptr;
194 
195  //open driver
196  GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
197 
198  if ( !outputDriver )
199  {
200  return outputDriver; //return nullptr, driver does not exist
201  }
202 
203  driverMetadata = GDALGetMetadata( outputDriver, nullptr );
204  if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
205  {
206  return nullptr; //driver exist, but it does not support the create operation
207  }
208 
209  return outputDriver;
210 }
211 
212 gdal::dataset_unique_ptr QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver )
213 {
214  //open output file
215  char **papszOptions = nullptr;
216  gdal::dataset_unique_ptr outputDataset( GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions ) );
217  if ( !outputDataset )
218  {
219  return nullptr;
220  }
221 
222  //assign georef information
223  double geotransform[6];
224  outputGeoTransform( geotransform );
225  GDALSetGeoTransform( outputDataset.get(), geotransform );
226 
227  return outputDataset;
228 }
229 
230 void QgsRasterCalculator::outputGeoTransform( double *transform ) const
231 {
232  transform[0] = mOutputRectangle.xMinimum();
233  transform[1] = mOutputRectangle.width() / mNumOutputColumns;
234  transform[2] = 0;
235  transform[3] = mOutputRectangle.yMaximum();
236  transform[4] = 0;
237  transform[5] = -mOutputRectangle.height() / mNumOutputRows;
238 }
A rectangle specified with double values.
Definition: qgsrectangle.h:39
bool isNumber() const
Returns true if matrix is 1x1 (=scalar number)
void cancel()
Tells the internal routines that the current operation should be canceled. This should be run by the ...
Definition: qgsfeedback.h:85
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
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
void canceled()
Internal routines can connect to this signal if they use event loop.
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:142
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:120
static QgsRasterCalcNode * parseRasterCalcString(const QString &str, QString &parserErrorMsg)
double yMaximum() const
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:125
Exact, precise but slow.
Feedback object tailored for raster block reading.
double height() const
Returns the height of the rectangle.
Definition: qgsrectangle.h:149