QGIS API Documentation  2.99.0-Master (7d4f81d)
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 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 
59 int QgsRasterCalculator::processCalculation( QProgressDialog *p )
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  if ( p )
121  {
122  p->setMaximum( mNumOutputRows );
123  }
124 
125  QgsRasterMatrix resultMatrix;
126  resultMatrix.setNodataValue( outputNodataValue );
127 
128  //read / write line by line
129  for ( int i = 0; i < mNumOutputRows; ++i )
130  {
131  if ( p )
132  {
133  p->setValue( i );
134  }
135 
136  if ( p && p->wasCanceled() )
137  {
138  break;
139  }
140 
141  if ( calcNode->calculate( inputBlocks, resultMatrix, i ) )
142  {
143  bool resultIsNumber = resultMatrix.isNumber();
144  float *calcData = new float[mNumOutputColumns];
145 
146  for ( int j = 0; j < mNumOutputColumns; ++j )
147  {
148  calcData[j] = ( float )( resultIsNumber ? resultMatrix.number() : resultMatrix.data()[j] );
149  }
150 
151  //write scanline to the dataset
152  if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
153  {
154  QgsDebugMsg( "RasterIO error!" );
155  }
156 
157  delete[] calcData;
158  }
159 
160  }
161 
162  if ( p )
163  {
164  p->setValue( mNumOutputRows );
165  }
166 
167  //close datasets and release memory
168  delete calcNode;
169  qDeleteAll( inputBlocks );
170  inputBlocks.clear();
171 
172  if ( p && p->wasCanceled() )
173  {
174  //delete the dataset without closing (because it is faster)
175  GDALDeleteDataset( outputDriver, mOutputFile.toUtf8().constData() );
176  return static_cast< int >( Canceled );
177  }
178  GDALClose( outputDataset );
179 
180  return static_cast< int >( Success );
181 }
182 
184  : mNumOutputColumns( 0 )
185  , mNumOutputRows( 0 )
186 {
187 }
188 
189 GDALDriverH QgsRasterCalculator::openOutputDriver()
190 {
191  char **driverMetadata = nullptr;
192 
193  //open driver
194  GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
195 
196  if ( !outputDriver )
197  {
198  return outputDriver; //return nullptr, driver does not exist
199  }
200 
201  driverMetadata = GDALGetMetadata( outputDriver, nullptr );
202  if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
203  {
204  return nullptr; //driver exist, but it does not support the create operation
205  }
206 
207  return outputDriver;
208 }
209 
210 GDALDatasetH QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver )
211 {
212  //open output file
213  char **papszOptions = nullptr;
214  GDALDatasetH outputDataset = GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions );
215  if ( !outputDataset )
216  {
217  return outputDataset;
218  }
219 
220  //assign georef information
221  double geotransform[6];
222  outputGeoTransform( geotransform );
223  GDALSetGeoTransform( outputDataset, geotransform );
224 
225  return outputDataset;
226 }
227 
228 void QgsRasterCalculator::outputGeoTransform( double *transform ) const
229 {
230  transform[0] = mOutputRectangle.xMinimum();
231  transform[1] = mOutputRectangle.width() / mNumOutputColumns;
232  transform[2] = 0;
233  transform[3] = mOutputRectangle.yMaximum();
234  transform[4] = 0;
235  transform[5] = -mOutputRectangle.height() / mNumOutputRows;
236 }
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)
QgsRasterCalculator(const QString &formulaString, const QString &outputFile, const QString &outputFormat, const QgsRectangle &outputExtent, int nOutputColumns, int nOutputRows, const QVector< QgsRasterCalculatorEntry > &rasterEntries)
QgsRasterCalculator constructor.
User canceled calculation.
void setCrs(const QgsCoordinateReferenceSystem &srcCRS, const QgsCoordinateReferenceSystem &destCRS, int srcDatumTransform=-1, int destDatumTransform=-1)
set source and destination CRS
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:118
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
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:96
static QgsRasterCalcNode * parseRasterCalcString(const QString &str, QString &parserErrorMsg)
double yMaximum() const
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:101
Exact, precise but slow.
double height() const
Returns the height of the rectangle.
Definition: qgsrectangle.h:125