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