00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include "qgsrastercalculator.h"
00019 #include "qgsrastercalcnode.h"
00020 #include "qgsrasterlayer.h"
00021 #include "qgsrastermatrix.h"
00022 #include "cpl_string.h"
00023 #include <QProgressDialog>
00024
00025 #if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800
00026 #define TO8(x) (x).toUtf8().constData()
00027 #else
00028 #define TO8(x) (x).toLocal8Bit().constData()
00029 #endif
00030
00031 QgsRasterCalculator::QgsRasterCalculator( const QString& formulaString, const QString& outputFile, const QString& outputFormat,
00032 const QgsRectangle& outputExtent, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry>& rasterEntries ): mFormulaString( formulaString ), mOutputFile( outputFile ), mOutputFormat( outputFormat ),
00033 mOutputRectangle( outputExtent ), mNumOutputColumns( nOutputColumns ), mNumOutputRows( nOutputRows ), mRasterEntries( rasterEntries )
00034 {
00035 }
00036
00037 QgsRasterCalculator::~QgsRasterCalculator()
00038 {
00039 }
00040
00041 int QgsRasterCalculator::processCalculation( QProgressDialog* p )
00042 {
00043
00044 QString errorString;
00045 QgsRasterCalcNode* calcNode = QgsRasterCalcNode::parseRasterCalcString( mFormulaString, errorString );
00046 if( !calcNode )
00047 {
00048
00049 }
00050
00051 double targetGeoTransform[6];
00052 outputGeoTransform( targetGeoTransform );
00053
00054
00055 QMap< QString, GDALRasterBandH > mInputRasterBands;
00056 QMap< QString, QgsRasterMatrix* > inputScanLineData;
00057 QVector< GDALDatasetH > mInputDatasets;
00058
00059 QVector<QgsRasterCalculatorEntry>::const_iterator it = mRasterEntries.constBegin();
00060 for( ; it != mRasterEntries.constEnd(); ++it )
00061 {
00062 if( !it->raster )
00063 {
00064 return 2;
00065 }
00066 GDALDatasetH inputDataset = GDALOpen( TO8( it->raster->source() ), GA_ReadOnly );
00067 if( inputDataset == NULL )
00068 {
00069 return 2;
00070 }
00071 GDALRasterBandH inputRasterBand = GDALGetRasterBand( inputDataset, it->bandNumber );
00072 if( inputRasterBand == NULL )
00073 {
00074 return 2;
00075 }
00076
00077 int nodataSuccess;
00078 double nodataValue = GDALGetRasterNoDataValue( inputRasterBand, &nodataSuccess );
00079
00080 mInputDatasets.push_back( inputDataset );
00081 mInputRasterBands.insert( it->ref, inputRasterBand );
00082 inputScanLineData.insert( it->ref, new QgsRasterMatrix( mNumOutputColumns, 1, new float[mNumOutputColumns], nodataValue ) );
00083 }
00084
00085
00086 GDALDriverH outputDriver = openOutputDriver();
00087 if( outputDriver == NULL )
00088 {
00089 return 1;
00090 }
00091 GDALDatasetH outputDataset = openOutputFile( outputDriver );
00092 GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset, 1 );
00093
00094 float outputNodataValue = -FLT_MAX;
00095 GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
00096
00097 float* resultScanLine = ( float * ) CPLMalloc( sizeof( float ) * mNumOutputColumns );
00098
00099 if( p )
00100 {
00101 p->setMaximum( mNumOutputRows );
00102 }
00103
00104 QgsRasterMatrix resultMatrix;
00105
00106
00107 for( int i = 0; i < mNumOutputRows; ++i )
00108 {
00109 if( p )
00110 {
00111 p->setValue( i );
00112 }
00113
00114 if( p && p->wasCanceled() )
00115 {
00116 break;
00117 }
00118
00119
00120 QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin();
00121 for( ; bufferIt != inputScanLineData.end(); ++bufferIt )
00122 {
00123 double sourceTransformation[6];
00124 GDALRasterBandH sourceRasterBand = mInputRasterBands[bufferIt.key()];
00125 GDALGetGeoTransform( GDALGetBandDataset( sourceRasterBand ), sourceTransformation );
00126
00127 readRasterPart( targetGeoTransform, 0, i, mNumOutputColumns, 1, sourceTransformation, sourceRasterBand, bufferIt.value()->data() );
00128 }
00129
00130 if( calcNode->calculate( inputScanLineData, resultMatrix ) )
00131 {
00132 bool resultIsNumber = resultMatrix.isNumber();
00133 float* calcData;
00134
00135 if( resultIsNumber )
00136 {
00137 calcData = new float[mNumOutputColumns];
00138 for( int j = 0; j < mNumOutputColumns; ++j )
00139 {
00140 calcData[j] = resultMatrix.number();
00141 }
00142 }
00143 else
00144 {
00145 calcData = resultMatrix.data();
00146 }
00147
00148
00149 for( int j = 0; j < mNumOutputColumns; ++j )
00150 {
00151 if( calcData[j] == resultMatrix.nodataValue() )
00152 {
00153 calcData[j] = outputNodataValue;
00154 }
00155 }
00156
00157
00158 if( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
00159 {
00160 qWarning( "RasterIO error!" );
00161 }
00162
00163 if( resultIsNumber )
00164 {
00165 delete[] calcData;
00166 }
00167 }
00168
00169 }
00170
00171 if( p )
00172 {
00173 p->setValue( mNumOutputRows );
00174 }
00175
00176
00177 delete calcNode;
00178 QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin();
00179 for( ; bufferIt != inputScanLineData.end(); ++bufferIt )
00180 {
00181 delete bufferIt.value();
00182 }
00183 inputScanLineData.clear();
00184
00185 QVector< GDALDatasetH >::iterator datasetIt = mInputDatasets.begin();
00186 for( ; datasetIt != mInputDatasets.end(); ++ datasetIt )
00187 {
00188 GDALClose( *datasetIt );
00189 }
00190
00191 if( p && p->wasCanceled() )
00192 {
00193
00194 GDALDeleteDataset( outputDriver, mOutputFile.toLocal8Bit().data() );
00195 return 3;
00196 }
00197 GDALClose( outputDataset );
00198 CPLFree( resultScanLine );
00199 return 0;
00200 }
00201
00202 QgsRasterCalculator::QgsRasterCalculator()
00203 {
00204 }
00205
00206 GDALDriverH QgsRasterCalculator::openOutputDriver()
00207 {
00208 char **driverMetadata;
00209
00210
00211 GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
00212
00213 if( outputDriver == NULL )
00214 {
00215 return outputDriver;
00216 }
00217
00218 driverMetadata = GDALGetMetadata( outputDriver, NULL );
00219 if( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
00220 {
00221 return NULL;
00222 }
00223
00224 return outputDriver;
00225 }
00226
00227 GDALDatasetH QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver )
00228 {
00229
00230 char **papszOptions = NULL;
00231 GDALDatasetH outputDataset = GDALCreate( outputDriver, mOutputFile.toLocal8Bit().data(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions );
00232 if( outputDataset == NULL )
00233 {
00234 return outputDataset;
00235 }
00236
00237
00238 double geotransform[6];
00239 outputGeoTransform( geotransform );
00240 GDALSetGeoTransform( outputDataset, geotransform );
00241
00242 return outputDataset;
00243 }
00244
00245 void QgsRasterCalculator::readRasterPart( double* targetGeotransform, int xOffset, int yOffset, int nCols, int nRows, double* sourceTransform, GDALRasterBandH sourceBand, float* rasterBuffer )
00246 {
00247
00248 if( transformationsEqual( targetGeotransform, sourceTransform ) )
00249 {
00250 GDALRasterIO( sourceBand, GF_Read, xOffset, yOffset, nCols, nRows, rasterBuffer, nCols, nRows, GDT_Float32, 0, 0 );
00251 return;
00252 }
00253
00254
00255 int nodataSuccess;
00256 double nodataValue = GDALGetRasterNoDataValue( sourceBand, &nodataSuccess );
00257 QgsRectangle targetRect( targetGeotransform[0] + targetGeotransform[1] * xOffset, targetGeotransform[3] + yOffset * targetGeotransform[5] + nRows * targetGeotransform[5]
00258 , targetGeotransform[0] + targetGeotransform[1] * xOffset + targetGeotransform[1] * nCols, targetGeotransform[3] + yOffset * targetGeotransform[5] );
00259 QgsRectangle sourceRect( sourceTransform[0], sourceTransform[3] + GDALGetRasterBandYSize( sourceBand ) * sourceTransform[5],
00260 sourceTransform[0] + GDALGetRasterBandXSize( sourceBand )* sourceTransform[1], sourceTransform[3] );
00261 QgsRectangle intersection = targetRect.intersect( &sourceRect );
00262
00263
00264 if( intersection.isEmpty() )
00265 {
00266 int nPixels = nCols * nRows;
00267 for( int i = 0; i < nPixels; ++i )
00268 {
00269 rasterBuffer[i] = nodataValue;
00270 }
00271 return;
00272 }
00273
00274
00275 int sourcePixelOffsetXMin = floor(( intersection.xMinimum() - sourceTransform[0] ) / sourceTransform[1] );
00276 int sourcePixelOffsetXMax = ceil(( intersection.xMaximum() - sourceTransform[0] ) / sourceTransform[1] );
00277 int nSourcePixelsX = sourcePixelOffsetXMax - sourcePixelOffsetXMin;
00278 int sourcePixelOffsetYMax = floor(( intersection.yMaximum() - sourceTransform[3] ) / sourceTransform[5] );
00279 int sourcePixelOffsetYMin = ceil(( intersection.yMinimum() - sourceTransform[3] ) / sourceTransform[5] );
00280 int nSourcePixelsY = sourcePixelOffsetYMin - sourcePixelOffsetYMax;
00281 float* sourceRaster = ( float * ) CPLMalloc( sizeof( float ) * nSourcePixelsX * nSourcePixelsY );
00282 double sourceRasterXMin = sourceRect.xMinimum() + sourcePixelOffsetXMin * sourceTransform[1];
00283 double sourceRasterYMax = sourceRect.yMaximum() + sourcePixelOffsetYMax * sourceTransform[5];
00284 GDALRasterIO( sourceBand, GF_Read, sourcePixelOffsetXMin, sourcePixelOffsetYMax, nSourcePixelsX, nSourcePixelsY,
00285 sourceRaster, nSourcePixelsX, nSourcePixelsY, GDT_Float32, 0, 0 );
00286
00287
00288 double targetPixelX;
00289 double targetPixelXMin = targetGeotransform[0] + targetGeotransform[1] * xOffset + targetGeotransform[1] / 2.0;
00290 double targetPixelY = targetGeotransform[3] + targetGeotransform[5] * yOffset + targetGeotransform[5] / 2.0;
00291 int sourceIndexX, sourceIndexY;
00292 double sx, sy;
00293 for( int i = 0; i < nRows; ++i )
00294 {
00295 targetPixelX = targetPixelXMin;
00296 for( int j = 0; j < nCols; ++j )
00297 {
00298 sx = ( targetPixelX - sourceRasterXMin ) / sourceTransform[1];
00299 sourceIndexX = sx > 0 ? sx : floor( sx );
00300 sy = ( targetPixelY - sourceRasterYMax ) / sourceTransform[5];
00301 sourceIndexY = sy > 0 ? sy : floor( sy );
00302 if( sourceIndexX >= 0 && sourceIndexX < nSourcePixelsX
00303 && sourceIndexY >= 0 && sourceIndexY < nSourcePixelsY )
00304 {
00305 rasterBuffer[j + i*nRows] = sourceRaster[ sourceIndexX + nSourcePixelsX * sourceIndexY ];
00306 }
00307 else
00308 {
00309 rasterBuffer[j + i*j] = nodataValue;
00310 }
00311 targetPixelX += targetGeotransform[1];
00312 }
00313 targetPixelY += targetGeotransform[5];
00314 }
00315
00316 CPLFree( sourceRaster );
00317 return;
00318 }
00319
00320 bool QgsRasterCalculator::transformationsEqual( double* t1, double* t2 ) const
00321 {
00322 for( int i = 0; i < 6; ++i )
00323 {
00324 if( !doubleNear( t1[i], t2[i], 0.00001 ) )
00325 {
00326 return false;
00327 }
00328 }
00329 return true;
00330 }
00331
00332 void QgsRasterCalculator::outputGeoTransform( double* transform ) const
00333 {
00334 transform[0] = mOutputRectangle.xMinimum();
00335 transform[1] = mOutputRectangle.width() / mNumOutputColumns;
00336 transform[2] = 0;
00337 transform[3] = mOutputRectangle.yMaximum();
00338 transform[4] = 0;
00339 transform[5] = -mOutputRectangle.height() / mNumOutputRows;
00340 }
00341
00342