QGIS API Documentation  2.11.0-Master
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 
64 
66 {
67 }
68 
70 {
71  //prepare search string / tree
72  QString errorString;
73  QgsRasterCalcNode* calcNode = QgsRasterCalcNode::parseRasterCalcString( mFormulaString, errorString );
74  if ( !calcNode )
75  {
76  //error
77  return 4;
78  }
79 
82  for ( ; it != mRasterEntries.constEnd(); ++it )
83  {
84  if ( !it->raster ) // no raster layer in entry
85  {
86  return 2;
87  }
88 
89  QgsRasterBlock* block = 0;
90  // if crs transform needed
91  if ( it->raster->crs() != mOutputCrs )
92  {
94  proj->setCRS( it->raster->crs(), mOutputCrs );
95  proj->setInput( it->raster->dataProvider()->clone() );
97 
98  block = proj->block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows );
99  }
100  else
101  {
102  block = it->raster->dataProvider()->block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows );
103  }
104  inputBlocks.insert( it->ref, block );
105  }
106 
107  //open output dataset for writing
108  GDALDriverH outputDriver = openOutputDriver();
109  if ( outputDriver == NULL )
110  {
111  return 1;
112  }
113 
114  GDALDatasetH outputDataset = openOutputFile( outputDriver );
115  GDALSetProjection( outputDataset, mOutputCrs.toWkt().toLocal8Bit().data() );
116  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset, 1 );
117 
118  float outputNodataValue = -FLT_MAX;
119  GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
120 
121  if ( p )
122  {
123  p->setMaximum( mNumOutputRows );
124  }
125 
126  QgsRasterMatrix resultMatrix;
127  resultMatrix.setNodataValue( outputNodataValue );
128 
129  //read / write line by line
130  for ( int i = 0; i < mNumOutputRows; ++i )
131  {
132  if ( p )
133  {
134  p->setValue( i );
135  }
136 
137  if ( p && p->wasCanceled() )
138  {
139  break;
140  }
141 
142  if ( calcNode->calculate( inputBlocks, resultMatrix, i ) )
143  {
144  bool resultIsNumber = resultMatrix.isNumber();
145  float* calcData = new float[mNumOutputColumns];
146 
147  for ( int j = 0; j < mNumOutputColumns; ++j )
148  {
149  calcData[j] = ( float )( resultIsNumber ? resultMatrix.number() : resultMatrix.data()[j] );
150  }
151 
152  //write scanline to the dataset
153  if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
154  {
155  qWarning( "RasterIO error!" );
156  }
157 
158  if ( resultIsNumber )
159  {
160  delete[] calcData;
161  }
162  }
163 
164  }
165 
166  if ( p )
167  {
168  p->setValue( mNumOutputRows );
169  }
170 
171  //close datasets and release memory
172  delete calcNode;
173  qDeleteAll( inputBlocks );
174  inputBlocks.clear();
175 
176  if ( p && p->wasCanceled() )
177  {
178  //delete the dataset without closing (because it is faster)
179  GDALDeleteDataset( outputDriver, TO8F( mOutputFile ) );
180  return 3;
181  }
182  GDALClose( outputDataset );
183 
184  return 0;
185 }
186 
187 QgsRasterCalculator::QgsRasterCalculator()
188  : mNumOutputColumns( 0 )
189  , mNumOutputRows( 0 )
190 {
191 }
192 
193 GDALDriverH QgsRasterCalculator::openOutputDriver()
194 {
195  char **driverMetadata;
196 
197  //open driver
198  GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
199 
200  if ( outputDriver == NULL )
201  {
202  return outputDriver; //return NULL, driver does not exist
203  }
204 
205  driverMetadata = GDALGetMetadata( outputDriver, NULL );
206  if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
207  {
208  return NULL; //driver exist, but it does not support the create operation
209  }
210 
211  return outputDriver;
212 }
213 
214 GDALDatasetH QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver )
215 {
216  //open output file
217  char **papszOptions = NULL;
218  GDALDatasetH outputDataset = GDALCreate( outputDriver, TO8F( mOutputFile ), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions );
219  if ( outputDataset == NULL )
220  {
221  return outputDataset;
222  }
223 
224  //assign georef information
225  double geotransform[6];
226  outputGeoTransform( geotransform );
227  GDALSetGeoTransform( outputDataset, geotransform );
228 
229  return outputDataset;
230 }
231 
232 void QgsRasterCalculator::outputGeoTransform( double* transform ) const
233 {
234  transform[0] = mOutputRectangle.xMinimum();
235  transform[1] = mOutputRectangle.width() / mNumOutputColumns;
236  transform[2] = 0;
237  transform[3] = mOutputRectangle.yMaximum();
238  transform[4] = 0;
239  transform[5] = -mOutputRectangle.height() / mNumOutputRows;
240 }
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:192
void setNodataValue(double d)
const_iterator constEnd() const
int processCalculation(QProgressDialog *p=0)
Starts the calculation and writes new raster.
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)
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:202
double xMinimum() const
Get the x minimum value (left side of rectangle)
Definition: qgsrectangle.h:187
double height() const
Height of the rectangle.
Definition: qgsrectangle.h:207