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