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