QGIS API Documentation  2.9.0-Master
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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 #include <ogr_srs_api.h>
29 
30 #if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800
31 #define TO8(x) (x).toUtf8().constData()
32 #define TO8F(x) (x).toUtf8().constData()
33 #else
34 #define TO8(x) (x).toLocal8Bit().constData()
35 #define TO8F(x) QFile::encodeName( x ).constData()
36 #endif
37 
38 QgsRasterCalculator::QgsRasterCalculator( const QString& formulaString, const QString& outputFile, const QString& outputFormat,
39  const QgsRectangle& outputExtent, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry>& rasterEntries )
40  : mFormulaString( formulaString )
41  , mOutputFile( outputFile )
42  , mOutputFormat( outputFormat )
43  , mOutputRectangle( outputExtent )
44  , mNumOutputColumns( nOutputColumns )
45  , mNumOutputRows( nOutputRows )
46  , mRasterEntries( rasterEntries )
47 {
48 }
49 
51 {
52 }
53 
54 int QgsRasterCalculator::processCalculation( QProgressDialog* p )
55 {
56  //prepare search string / tree
57  QString errorString;
58  QgsRasterCalcNode* calcNode = QgsRasterCalcNode::parseRasterCalcString( mFormulaString, errorString );
59  if ( !calcNode )
60  {
61  //error
62  return 4;
63  }
64 
65  double targetGeoTransform[6];
66  outputGeoTransform( targetGeoTransform );
67 
68  //open all input rasters for reading
69  QMap< QString, GDALRasterBandH > mInputRasterBands; //raster references and corresponding scanline data
70  QMap< QString, QgsRasterMatrix* > inputScanLineData; //stores raster references and corresponding scanline data
71  QVector< GDALDatasetH > mInputDatasets; //raster references and corresponding dataset
72 
73  QVector<QgsRasterCalculatorEntry>::const_iterator it = mRasterEntries.constBegin();
74  for ( ; it != mRasterEntries.constEnd(); ++it )
75  {
76  if ( !it->raster ) // no raster layer in entry
77  {
78  return 2;
79  }
80  GDALDatasetH inputDataset = GDALOpen( TO8F( it->raster->source() ), GA_ReadOnly );
81  if ( !inputDataset )
82  {
83  return 2;
84  }
85 
86  //check if the input dataset is south up or rotated. If yes, use GDALAutoCreateWarpedVRT to create a north up raster
87  double inputGeoTransform[6];
88  if ( GDALGetGeoTransform( inputDataset, inputGeoTransform ) == CE_None
89  && ( inputGeoTransform[1] < 0.0
90  || inputGeoTransform[2] != 0.0
91  || inputGeoTransform[4] != 0.0
92  || inputGeoTransform[5] > 0.0 ) )
93  {
94  GDALDatasetH vDataset = GDALAutoCreateWarpedVRT( inputDataset, NULL, NULL, GRA_NearestNeighbour, 0.2, NULL );
95  mInputDatasets.push_back( vDataset );
96  mInputDatasets.push_back( inputDataset );
97  inputDataset = vDataset;
98  }
99  else
100  {
101  mInputDatasets.push_back( inputDataset );
102  }
103 
104  GDALRasterBandH inputRasterBand = GDALGetRasterBand( inputDataset, it->bandNumber );
105  if ( inputRasterBand == NULL )
106  {
107  return 2;
108  }
109 
110  int nodataSuccess;
111  double nodataValue = GDALGetRasterNoDataValue( inputRasterBand, &nodataSuccess );
112 
113  mInputRasterBands.insert( it->ref, inputRasterBand );
114  inputScanLineData.insert( it->ref, new QgsRasterMatrix( mNumOutputColumns, 1, new float[mNumOutputColumns], nodataValue ) );
115  }
116 
117  //open output dataset for writing
118  GDALDriverH outputDriver = openOutputDriver();
119  if ( outputDriver == NULL )
120  {
121  return 1;
122  }
123  GDALDatasetH outputDataset = openOutputFile( outputDriver );
124 
125  //copy the projection info from the first input raster
126  if ( mRasterEntries.size() > 0 )
127  {
128  QgsRasterLayer* rl = mRasterEntries.at( 0 ).raster;
129  if ( rl )
130  {
131  char* crsWKT = 0;
132  OGRSpatialReferenceH ogrSRS = OSRNewSpatialReference( NULL );
133  if ( OSRSetFromUserInput( ogrSRS, rl->crs().authid().toUtf8().constData() ) == OGRERR_NONE )
134  {
135  OSRExportToWkt( ogrSRS, &crsWKT );
136  GDALSetProjection( outputDataset, crsWKT );
137  }
138  else
139  {
140  GDALSetProjection( outputDataset, TO8( rl->crs().toWkt() ) );
141  }
142  OSRDestroySpatialReference( ogrSRS );
143  CPLFree( crsWKT );
144  }
145  }
146 
147 
148  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset, 1 );
149 
150  float outputNodataValue = -FLT_MAX;
151  GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
152 
153  float* resultScanLine = ( float * ) CPLMalloc( sizeof( float ) * mNumOutputColumns );
154 
155  if ( p )
156  {
157  p->setMaximum( mNumOutputRows );
158  }
159 
160  QgsRasterMatrix resultMatrix;
161 
162  //read / write line by line
163  for ( int i = 0; i < mNumOutputRows; ++i )
164  {
165  if ( p )
166  {
167  p->setValue( i );
168  }
169 
170  if ( p && p->wasCanceled() )
171  {
172  break;
173  }
174 
175  //fill buffers
176  QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin();
177  for ( ; bufferIt != inputScanLineData.end(); ++bufferIt )
178  {
179  double sourceTransformation[6];
180  GDALRasterBandH sourceRasterBand = mInputRasterBands[bufferIt.key()];
181  if ( GDALGetGeoTransform( GDALGetBandDataset( sourceRasterBand ), sourceTransformation ) != CE_None )
182  {
183  qWarning( "GDALGetGeoTransform failed!" );
184  }
185 
186  //the function readRasterPart calls GDALRasterIO (and ev. does some conversion if raster transformations are not the same)
187  readRasterPart( targetGeoTransform, 0, i, mNumOutputColumns, 1, sourceTransformation, sourceRasterBand, bufferIt.value()->data() );
188  }
189 
190  if ( calcNode->calculate( inputScanLineData, resultMatrix ) )
191  {
192  bool resultIsNumber = resultMatrix.isNumber();
193  float* calcData;
194 
195  if ( resultIsNumber ) //scalar result. Insert number for every pixel
196  {
197  calcData = new float[mNumOutputColumns];
198  for ( int j = 0; j < mNumOutputColumns; ++j )
199  {
200  calcData[j] = resultMatrix.number();
201  }
202  }
203  else //result is real matrix
204  {
205  calcData = resultMatrix.data();
206  }
207 
208  //replace all matrix nodata values with output nodatas
209  for ( int j = 0; j < mNumOutputColumns; ++j )
210  {
211  if ( calcData[j] == resultMatrix.nodataValue() )
212  {
213  calcData[j] = outputNodataValue;
214  }
215  }
216 
217  //write scanline to the dataset
218  if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
219  {
220  qWarning( "RasterIO error!" );
221  }
222 
223  if ( resultIsNumber )
224  {
225  delete[] calcData;
226  }
227  }
228 
229  }
230 
231  if ( p )
232  {
233  p->setValue( mNumOutputRows );
234  }
235 
236  //close datasets and release memory
237  delete calcNode;
238  QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin();
239  for ( ; bufferIt != inputScanLineData.end(); ++bufferIt )
240  {
241  delete bufferIt.value();
242  }
243  inputScanLineData.clear();
244 
245  QVector< GDALDatasetH >::iterator datasetIt = mInputDatasets.begin();
246  for ( ; datasetIt != mInputDatasets.end(); ++ datasetIt )
247  {
248  GDALClose( *datasetIt );
249  }
250 
251  if ( p && p->wasCanceled() )
252  {
253  //delete the dataset without closing (because it is faster)
254  GDALDeleteDataset( outputDriver, TO8F( mOutputFile ) );
255  return 3;
256  }
257  GDALClose( outputDataset );
258  CPLFree( resultScanLine );
259  return 0;
260 }
261 
262 QgsRasterCalculator::QgsRasterCalculator()
263  : mNumOutputColumns( 0 )
264  , mNumOutputRows( 0 )
265 {
266 }
267 
268 GDALDriverH QgsRasterCalculator::openOutputDriver()
269 {
270  char **driverMetadata;
271 
272  //open driver
273  GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
274 
275  if ( outputDriver == NULL )
276  {
277  return outputDriver; //return NULL, driver does not exist
278  }
279 
280  driverMetadata = GDALGetMetadata( outputDriver, NULL );
281  if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
282  {
283  return NULL; //driver exist, but it does not support the create operation
284  }
285 
286  return outputDriver;
287 }
288 
289 GDALDatasetH QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver )
290 {
291  //open output file
292  char **papszOptions = NULL;
293  GDALDatasetH outputDataset = GDALCreate( outputDriver, TO8F( mOutputFile ), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions );
294  if ( outputDataset == NULL )
295  {
296  return outputDataset;
297  }
298 
299  //assign georef information
300  double geotransform[6];
301  outputGeoTransform( geotransform );
302  GDALSetGeoTransform( outputDataset, geotransform );
303 
304  return outputDataset;
305 }
306 
307 void QgsRasterCalculator::readRasterPart( double* targetGeotransform, int xOffset, int yOffset, int nCols, int nRows, double* sourceTransform, GDALRasterBandH sourceBand, float* rasterBuffer )
308 {
309  //If dataset transform is the same as the requested transform, do a normal GDAL raster io
310  if ( transformationsEqual( targetGeotransform, sourceTransform ) )
311  {
312  GDALRasterIO( sourceBand, GF_Read, xOffset, yOffset, nCols, nRows, rasterBuffer, nCols, nRows, GDT_Float32, 0, 0 );
313  return;
314  }
315 
316  int sourceBandXSize = GDALGetRasterBandXSize( sourceBand );
317  int sourceBandYSize = GDALGetRasterBandYSize( sourceBand );
318 
319  //pixel calculation needed because of different raster position / resolution
320  int nodataSuccess;
321  double nodataValue = GDALGetRasterNoDataValue( sourceBand, &nodataSuccess );
322  QgsRectangle targetRect( targetGeotransform[0] + targetGeotransform[1] * xOffset, targetGeotransform[3] + yOffset * targetGeotransform[5] + nRows * targetGeotransform[5]
323  , targetGeotransform[0] + targetGeotransform[1] * xOffset + targetGeotransform[1] * nCols, targetGeotransform[3] + yOffset * targetGeotransform[5] );
324  QgsRectangle sourceRect( sourceTransform[0], sourceTransform[3] + GDALGetRasterBandYSize( sourceBand ) * sourceTransform[5],
325  sourceTransform[0] + GDALGetRasterBandXSize( sourceBand )* sourceTransform[1], sourceTransform[3] );
326  QgsRectangle intersection = targetRect.intersect( &sourceRect );
327 
328  //no intersection, fill all the pixels with nodata values
329  if ( intersection.isEmpty() )
330  {
331  int nPixels = nCols * nRows;
332  for ( int i = 0; i < nPixels; ++i )
333  {
334  rasterBuffer[i] = nodataValue;
335  }
336  return;
337  }
338 
339  //do raster io in source resolution
340  int sourcePixelOffsetXMin = floor(( intersection.xMinimum() - sourceTransform[0] ) / sourceTransform[1] );
341  int sourcePixelOffsetXMax = ceil(( intersection.xMaximum() - sourceTransform[0] ) / sourceTransform[1] );
342  if ( sourcePixelOffsetXMax > sourceBandXSize )
343  {
344  sourcePixelOffsetXMax = sourceBandXSize;
345  }
346  int nSourcePixelsX = sourcePixelOffsetXMax - sourcePixelOffsetXMin;
347 
348  int sourcePixelOffsetYMax = floor(( intersection.yMaximum() - sourceTransform[3] ) / sourceTransform[5] );
349  int sourcePixelOffsetYMin = ceil(( intersection.yMinimum() - sourceTransform[3] ) / sourceTransform[5] );
350  if ( sourcePixelOffsetYMin > sourceBandYSize )
351  {
352  sourcePixelOffsetYMin = sourceBandYSize;
353  }
354  int nSourcePixelsY = sourcePixelOffsetYMin - sourcePixelOffsetYMax;
355  float* sourceRaster = ( float * ) CPLMalloc( sizeof( float ) * nSourcePixelsX * nSourcePixelsY );
356  double sourceRasterXMin = sourceRect.xMinimum() + sourcePixelOffsetXMin * sourceTransform[1];
357  double sourceRasterYMax = sourceRect.yMaximum() + sourcePixelOffsetYMax * sourceTransform[5];
358  if ( GDALRasterIO( sourceBand, GF_Read, sourcePixelOffsetXMin, sourcePixelOffsetYMax, nSourcePixelsX, nSourcePixelsY,
359  sourceRaster, nSourcePixelsX, nSourcePixelsY, GDT_Float32, 0, 0 ) != CE_None )
360  {
361  //IO error, fill array with nodata values
362  CPLFree( sourceRaster );
363  int npixels = nRows * nCols;
364  for ( int i = 0; i < npixels; ++i )
365  {
366  rasterBuffer[i] = nodataValue;
367  }
368  return;
369  }
370 
371 
372  double targetPixelX;
373  double targetPixelXMin = targetGeotransform[0] + targetGeotransform[1] * xOffset + targetGeotransform[1] / 2.0;
374  double targetPixelY = targetGeotransform[3] + targetGeotransform[5] * yOffset + targetGeotransform[5] / 2.0; //coordinates of current target pixel
375  int sourceIndexX, sourceIndexY; //current raster index in source pixels
376  double sx, sy;
377  for ( int i = 0; i < nRows; ++i )
378  {
379  targetPixelX = targetPixelXMin;
380  for ( int j = 0; j < nCols; ++j )
381  {
382  sx = ( targetPixelX - sourceRasterXMin ) / sourceTransform[1];
383  sourceIndexX = sx > 0 ? sx : floor( sx );
384  sy = ( targetPixelY - sourceRasterYMax ) / sourceTransform[5];
385  sourceIndexY = sy > 0 ? sy : floor( sy );
386  if ( sourceIndexX >= 0 && sourceIndexX < nSourcePixelsX
387  && sourceIndexY >= 0 && sourceIndexY < nSourcePixelsY )
388  {
389  rasterBuffer[j + i*nRows] = sourceRaster[ sourceIndexX + nSourcePixelsX * sourceIndexY ];
390  }
391  else
392  {
393  rasterBuffer[j + i*j] = nodataValue;
394  }
395  targetPixelX += targetGeotransform[1];
396  }
397  targetPixelY += targetGeotransform[5];
398  }
399 
400  CPLFree( sourceRaster );
401  return;
402 }
403 
404 bool QgsRasterCalculator::transformationsEqual( double* t1, double* t2 ) const
405 {
406  for ( int i = 0; i < 6; ++i )
407  {
408  if ( !qgsDoubleNear( t1[i], t2[i], 0.00001 ) )
409  {
410  return false;
411  }
412  }
413  return true;
414 }
415 
416 void QgsRasterCalculator::outputGeoTransform( double* transform ) const
417 {
418  transform[0] = mOutputRectangle.xMinimum();
419  transform[1] = mOutputRectangle.width() / mNumOutputColumns;
420  transform[2] = 0;
421  transform[3] = mOutputRectangle.yMaximum();
422  transform[4] = 0;
423  transform[5] = -mOutputRectangle.height() / mNumOutputRows;
424 }
A rectangle specified with double values.
Definition: qgsrectangle.h:35
bool isEmpty() const
test if rectangle is empty.
#define TO8(x)
double yMaximum() const
Get the y maximum value (top side of rectangle)
Definition: qgsrectangle.h:188
This class provides qgis with the ability to render raster datasets onto the mapcanvas.
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)
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Definition: qgis.h:327
#define TO8F(x)
double number() const
double yMinimum() const
Get the y minimum value (bottom side of rectangle)
Definition: qgsrectangle.h:193
double xMaximum() const
Get the x maximum value (right side of rectangle)
Definition: qgsrectangle.h:178
bool calculate(QMap< QString, QgsRasterMatrix * > &rasterData, QgsRasterMatrix &result) const
Calculates result (might be real matrix or single number)
bool isNumber() const
Returns true if matrix is 1x1 (=scalar number)
double nodataValue() const
QgsRectangle intersect(const QgsRectangle *rect) const
return the intersection with the given rectangle
static QgsRasterCalcNode * parseRasterCalcString(const QString &str, QString &parserErrorMsg)
const QgsCoordinateReferenceSystem & crs() const
Returns layer's spatial reference system.
float * data()
Returns data array (but not ownership)
double width() const
Width of the rectangle.
Definition: qgsrectangle.h:198
double xMinimum() const
Get the x minimum value (left side of rectangle)
Definition: qgsrectangle.h:183
double height() const
Height of the rectangle.
Definition: qgsrectangle.h:203
void * OGRSpatialReferenceH