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