QGIS API Documentation  2.14.0-Essen
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 "qgsrasterlayer.h"
21 #include "qgsrastermatrix.h"
22 
23 #include <QProgressDialog>
24 #include <QFile>
25 
26 #include <cpl_string.h>
27 #include <gdalwarper.h>
28 
29 #if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800
30 #define TO8(x) (x).toUtf8().constData()
31 #define TO8F(x) (x).toUtf8().constData()
32 #else
33 #define TO8(x) (x).toLocal8Bit().constData()
34 #define TO8F(x) QFile::encodeName( x ).constData()
35 #endif
36 
37 QgsRasterCalculator::QgsRasterCalculator( const QString& formulaString, const QString& outputFile, const QString& outputFormat,
38  const QgsRectangle& outputExtent, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry>& rasterEntries )
39  : mFormulaString( formulaString )
40  , mOutputFile( outputFile )
41  , mOutputFormat( outputFormat )
42  , mOutputRectangle( outputExtent )
43  , mNumOutputColumns( nOutputColumns )
44  , mNumOutputRows( nOutputRows )
45  , mRasterEntries( rasterEntries )
46 {
47  //default to first layer's crs
48  mOutputCrs = mRasterEntries.at( 0 ).raster->crs();
49 }
50 
51 QgsRasterCalculator::QgsRasterCalculator( const QString& formulaString, const QString& outputFile, const QString& outputFormat,
52  const QgsRectangle& outputExtent, const QgsCoordinateReferenceSystem& outputCrs, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry>& rasterEntries )
53  : mFormulaString( formulaString )
54  , mOutputFile( outputFile )
55  , mOutputFormat( outputFormat )
56  , mOutputRectangle( outputExtent )
57  , mOutputCrs( outputCrs )
58  , mNumOutputColumns( nOutputColumns )
59  , mNumOutputRows( nOutputRows )
60  , mRasterEntries( rasterEntries )
61 {
62 }
63 
65 {
66  //prepare search string / tree
67  QString errorString;
68  QgsRasterCalcNode* calcNode = QgsRasterCalcNode::parseRasterCalcString( mFormulaString, errorString );
69  if ( !calcNode )
70  {
71  //error
72  return static_cast<int>( ParserError );
73  }
74 
77  for ( ; it != mRasterEntries.constEnd(); ++it )
78  {
79  if ( !it->raster ) // no raster layer in entry
80  {
81  delete calcNode;
82  qDeleteAll( inputBlocks );
83  return static_cast< int >( InputLayerError );
84  }
85 
86  QgsRasterBlock* block = nullptr;
87  // if crs transform needed
88  if ( it->raster->crs() != mOutputCrs )
89  {
90  QgsRasterProjector proj;
91  proj.setCRS( it->raster->crs(), mOutputCrs );
92  proj.setInput( it->raster->dataProvider() );
94 
95  block = proj.block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows );
96  }
97  else
98  {
99  block = it->raster->dataProvider()->block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows );
100  }
101  if ( block->isEmpty() )
102  {
103  delete block;
104  delete calcNode;
105  qDeleteAll( inputBlocks );
106  return static_cast<int>( MemoryError );
107  }
108  inputBlocks.insert( it->ref, block );
109  }
110 
111  //open output dataset for writing
112  GDALDriverH outputDriver = openOutputDriver();
113  if ( !outputDriver )
114  {
115  return static_cast< int >( CreateOutputError );
116  }
117 
118  GDALDatasetH outputDataset = openOutputFile( outputDriver );
119  GDALSetProjection( outputDataset, mOutputCrs.toWkt().toLocal8Bit().data() );
120  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset, 1 );
121 
122  float outputNodataValue = -FLT_MAX;
123  GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
124 
125  if ( p )
126  {
127  p->setMaximum( mNumOutputRows );
128  }
129 
130  QgsRasterMatrix resultMatrix;
131  resultMatrix.setNodataValue( outputNodataValue );
132 
133  //read / write line by line
134  for ( int i = 0; i < mNumOutputRows; ++i )
135  {
136  if ( p )
137  {
138  p->setValue( i );
139  }
140 
141  if ( p && p->wasCanceled() )
142  {
143  break;
144  }
145 
146  if ( calcNode->calculate( inputBlocks, resultMatrix, i ) )
147  {
148  bool resultIsNumber = resultMatrix.isNumber();
149  float* calcData = new float[mNumOutputColumns];
150 
151  for ( int j = 0; j < mNumOutputColumns; ++j )
152  {
153  calcData[j] = ( float )( resultIsNumber ? resultMatrix.number() : resultMatrix.data()[j] );
154  }
155 
156  //write scanline to the dataset
157  if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
158  {
159  qWarning( "RasterIO error!" );
160  }
161 
162  delete[] calcData;
163  }
164 
165  }
166 
167  if ( p )
168  {
169  p->setValue( mNumOutputRows );
170  }
171 
172  //close datasets and release memory
173  delete calcNode;
174  qDeleteAll( inputBlocks );
175  inputBlocks.clear();
176 
177  if ( p && p->wasCanceled() )
178  {
179  //delete the dataset without closing (because it is faster)
180  GDALDeleteDataset( outputDriver, TO8F( mOutputFile ) );
181  return static_cast< int >( Cancelled );
182  }
183  GDALClose( outputDataset );
184 
185  return static_cast< int >( Success );
186 }
187 
189  : mNumOutputColumns( 0 )
190  , mNumOutputRows( 0 )
191 {
192 }
193 
194 GDALDriverH QgsRasterCalculator::openOutputDriver()
195 {
196  char **driverMetadata;
197 
198  //open driver
199  GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
200 
201  if ( !outputDriver )
202  {
203  return outputDriver; //return nullptr, driver does not exist
204  }
205 
206  driverMetadata = GDALGetMetadata( outputDriver, nullptr );
207  if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
208  {
209  return nullptr; //driver exist, but it does not support the create operation
210  }
211 
212  return outputDriver;
213 }
214 
215 GDALDatasetH QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver )
216 {
217  //open output file
218  char **papszOptions = nullptr;
219  GDALDatasetH outputDataset = GDALCreate( outputDriver, TO8F( mOutputFile ), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions );
220  if ( !outputDataset )
221  {
222  return outputDataset;
223  }
224 
225  //assign georef information
226  double geotransform[6];
227  outputGeoTransform( geotransform );
228  GDALSetGeoTransform( outputDataset, geotransform );
229 
230  return outputDataset;
231 }
232 
233 void QgsRasterCalculator::outputGeoTransform( double* transform ) const
234 {
235  transform[0] = mOutputRectangle.xMinimum();
236  transform[1] = mOutputRectangle.width() / mNumOutputColumns;
237  transform[2] = 0;
238  transform[3] = mOutputRectangle.yMaximum();
239  transform[4] = 0;
240  transform[5] = -mOutputRectangle.height() / mNumOutputRows;
241 }
A rectangle specified with double values.
Definition: qgsrectangle.h:35
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 setCRS(const QgsCoordinateReferenceSystem &theSrcCRS, const QgsCoordinateReferenceSystem &theDestCRS, int srcDatumTransform=-1, int destDatumTransform=-1)
set source and destination CRS
void setMaximum(int maximum)
double yMaximum() const
Get the y maximum value (top side of rectangle)
Definition: qgsrectangle.h:197
void setNodataValue(double d)
const_iterator constEnd() const
QString toWkt() const
Returns a WKT representation of this CRS.
QgsRasterCalculator(const QString &formulaString, const QString &outputFile, const QString &outputFormat, const QgsRectangle &outputExtent, int nOutputColumns, int nOutputRows, const QVector< QgsRasterCalculatorEntry > &rasterEntries)
QgsRasterCalculator constructor.
void clear()
#define TO8F(x)
void setValue(int progress)
double number() const
Raster data container.
QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height) override
Read block of data using given extent and size.
void * GDALDatasetH
QByteArray toLocal8Bit() const
bool isNumber() const
Returns true if matrix is 1x1 (=scalar number)
void setPrecision(Precision precision)
int processCalculation(QProgressDialog *p=nullptr)
Starts the calculation and writes new raster.
virtual bool setInput(QgsRasterInterface *input)
Set input.
const T & at(int i) const
const_iterator constBegin() const
Class for storing a coordinate reference system (CRS)
double * data()
Returns data array (but not ownership)
char * data()
static QgsRasterCalcNode * parseRasterCalcString(const QString &str, QString &parserErrorMsg)
Exact, precise but slow.
iterator insert(const Key &key, const T &value)
double width() const
Width of the rectangle.
Definition: qgsrectangle.h:207
double xMinimum() const
Get the x minimum value (left side of rectangle)
Definition: qgsrectangle.h:192
double height() const
Height of the rectangle.
Definition: qgsrectangle.h:212
bool isEmpty() const
Returns true if block is empty, i.e.