22 #include "cpl_string.h"
23 #include <QProgressDialog>
25 #include "gdalwarper.h"
26 #include <ogr_srs_api.h>
28 #if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800
29 #define TO8(x) (x).toUtf8().constData()
31 #define TO8(x) (x).toLocal8Bit().constData()
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 )
54 double targetGeoTransform[6];
58 QMap< QString, GDALRasterBandH > mInputRasterBands;
59 QMap< QString, QgsRasterMatrix* > inputScanLineData;
60 QVector< GDALDatasetH > mInputDatasets;
62 QVector<QgsRasterCalculatorEntry>::const_iterator it =
mRasterEntries.constBegin();
69 GDALDatasetH inputDataset = GDALOpen(
TO8( it->raster->source() ), GA_ReadOnly );
70 if ( inputDataset == NULL )
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 ) )
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;
90 mInputDatasets.push_back( inputDataset );
94 GDALRasterBandH inputRasterBand = GDALGetRasterBand( inputDataset, it->bandNumber );
95 if ( inputRasterBand == NULL )
101 double nodataValue = GDALGetRasterNoDataValue( inputRasterBand, &nodataSuccess );
103 mInputRasterBands.insert( it->ref, inputRasterBand );
109 if ( outputDriver == NULL )
122 OGRSpatialReferenceH ogrSRS = OSRNewSpatialReference( NULL );
123 if ( OSRSetFromUserInput( ogrSRS, rl->
crs().
authid().toUtf8().constData() ) == OGRERR_NONE )
125 OSRExportToWkt( ogrSRS, &crsWKT );
126 GDALSetProjection( outputDataset, crsWKT );
130 GDALSetProjection( outputDataset,
TO8( rl->
crs().
toWkt() ) );
132 OSRDestroySpatialReference( ogrSRS );
138 GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset, 1 );
140 float outputNodataValue = -FLT_MAX;
141 GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
143 float* resultScanLine = (
float * ) CPLMalloc(
sizeof(
float ) *
mNumOutputColumns );
160 if ( p && p->wasCanceled() )
166 QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin();
167 for ( ; bufferIt != inputScanLineData.end(); ++bufferIt )
169 double sourceTransformation[6];
170 GDALRasterBandH sourceRasterBand = mInputRasterBands[bufferIt.key()];
171 GDALGetGeoTransform( GDALGetBandDataset( sourceRasterBand ), sourceTransformation );
176 if ( calcNode->
calculate( inputScanLineData, resultMatrix ) )
178 bool resultIsNumber = resultMatrix.
isNumber();
181 if ( resultIsNumber )
186 calcData[j] = resultMatrix.
number();
191 calcData = resultMatrix.
data();
199 calcData[j] = outputNodataValue;
204 if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
206 qWarning(
"RasterIO error!" );
209 if ( resultIsNumber )
219 p->setValue( mNumOutputRows );
224 QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin();
225 for ( ; bufferIt != inputScanLineData.end(); ++bufferIt )
227 delete bufferIt.value();
229 inputScanLineData.clear();
231 QVector< GDALDatasetH >::iterator datasetIt = mInputDatasets.begin();
232 for ( ; datasetIt != mInputDatasets.end(); ++ datasetIt )
234 GDALClose( *datasetIt );
237 if ( p && p->wasCanceled() )
240 GDALDeleteDataset( outputDriver,
mOutputFile.toLocal8Bit().data() );
243 GDALClose( outputDataset );
244 CPLFree( resultScanLine );
254 char **driverMetadata;
257 GDALDriverH outputDriver = GDALGetDriverByName(
mOutputFormat.toLocal8Bit().data() );
259 if ( outputDriver == NULL )
264 driverMetadata = GDALGetMetadata( outputDriver, NULL );
265 if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE,
false ) )
276 char **papszOptions = NULL;
278 if ( outputDataset == NULL )
280 return outputDataset;
284 double geotransform[6];
286 GDALSetGeoTransform( outputDataset, geotransform );
288 return outputDataset;
291 void QgsRasterCalculator::readRasterPart(
double* targetGeotransform,
int xOffset,
int yOffset,
int nCols,
int nRows,
double* sourceTransform, GDALRasterBandH sourceBand,
float* rasterBuffer )
296 GDALRasterIO( sourceBand, GF_Read, xOffset, yOffset, nCols, nRows, rasterBuffer, nCols, nRows, GDT_Float32, 0, 0 );
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] );
312 int nPixels = nCols * nRows;
313 for (
int i = 0; i < nPixels; ++i )
315 rasterBuffer[i] = nodataValue;
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 );
335 double targetPixelXMin = targetGeotransform[0] + targetGeotransform[1] * xOffset + targetGeotransform[1] / 2.0;
336 double targetPixelY = targetGeotransform[3] + targetGeotransform[5] * yOffset + targetGeotransform[5] / 2.0;
337 int sourceIndexX, sourceIndexY;
339 for (
int i = 0; i < nRows; ++i )
341 targetPixelX = targetPixelXMin;
342 for (
int j = 0; j < nCols; ++j )
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 )
351 rasterBuffer[j + i*nRows] = sourceRaster[ sourceIndexX + nSourcePixelsX * sourceIndexY ];
355 rasterBuffer[j + i*j] = nodataValue;
357 targetPixelX += targetGeotransform[1];
359 targetPixelY += targetGeotransform[5];
362 CPLFree( sourceRaster );
368 for (
int i = 0; i < 6; ++i )