QGIS API Documentation  2.12.0-Lyon
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  {
93  QgsRasterProjector proj;
94  proj.setCRS( it->raster->crs(), mOutputCrs );
95  proj.setInput( it->raster->dataProvider() );
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  delete[] calcData;
159  }
160 
161  }
162 
163  if ( p )
164  {
165  p->setValue( mNumOutputRows );
166  }
167 
168  //close datasets and release memory
169  delete calcNode;
170  qDeleteAll( inputBlocks );
171  inputBlocks.clear();
172 
173  if ( p && p->wasCanceled() )
174  {
175  //delete the dataset without closing (because it is faster)
176  GDALDeleteDataset( outputDriver, TO8F( mOutputFile ) );
177  return 3;
178  }
179  GDALClose( outputDataset );
180 
181  return 0;
182 }
183 
184 QgsRasterCalculator::QgsRasterCalculator()
185  : mNumOutputColumns( 0 )
186  , mNumOutputRows( 0 )
187 {
188 }
189 
190 GDALDriverH QgsRasterCalculator::openOutputDriver()
191 {
192  char **driverMetadata;
193 
194  //open driver
195  GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
196 
197  if ( outputDriver == NULL )
198  {
199  return outputDriver; //return NULL, driver does not exist
200  }
201 
202  driverMetadata = GDALGetMetadata( outputDriver, NULL );
203  if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
204  {
205  return NULL; //driver exist, but it does not support the create operation
206  }
207 
208  return outputDriver;
209 }
210 
211 GDALDatasetH QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver )
212 {
213  //open output file
214  char **papszOptions = NULL;
215  GDALDatasetH outputDataset = GDALCreate( outputDriver, TO8F( mOutputFile ), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions );
216  if ( outputDataset == NULL )
217  {
218  return outputDataset;
219  }
220 
221  //assign georef information
222  double geotransform[6];
223  outputGeoTransform( geotransform );
224  GDALSetGeoTransform( outputDataset, geotransform );
225 
226  return outputDataset;
227 }
228 
229 void QgsRasterCalculator::outputGeoTransform( double* transform ) const
230 {
231  transform[0] = mOutputRectangle.xMinimum();
232  transform[1] = mOutputRectangle.width() / mNumOutputColumns;
233  transform[2] = 0;
234  transform[3] = mOutputRectangle.yMaximum();
235  transform[4] = 0;
236  transform[5] = -mOutputRectangle.height() / mNumOutputRows;
237 }
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:196
void setNodataValue(double d)
const_iterator constEnd() const
QString toWkt() const
A helper to get an wkt representation of this srs.
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:206
double xMinimum() const
Get the x minimum value (left side of rectangle)
Definition: qgsrectangle.h:191
double height() const
Height of the rectangle.
Definition: qgsrectangle.h:211