From 202ec2bca478200fb59dd3c6ff154f131901ae88 Mon Sep 17 00:00:00 2001 From: mhugent Date: Sun, 5 Dec 2010 19:22:04 +0000 Subject: [PATCH] Use GDALAutoCreateWarpedVRT to handle south-up (or rotated) rasters also for calculator. Fixes bug #3281 git-svn-id: http://svn.osgeo.org/qgis/trunk@14842 c8812cc2-4d05-0410-92ff-de0c093fc19c --- src/analysis/raster/qgsrastercalculator.cpp | 93 +++++++++++++-------- 1 file changed, 57 insertions(+), 36 deletions(-) diff --git a/src/analysis/raster/qgsrastercalculator.cpp b/src/analysis/raster/qgsrastercalculator.cpp index bf2e70d6263..68894417021 100644 --- a/src/analysis/raster/qgsrastercalculator.cpp +++ b/src/analysis/raster/qgsrastercalculator.cpp @@ -22,9 +22,11 @@ #include "cpl_string.h" #include +#include "gdalwarper.h" + QgsRasterCalculator::QgsRasterCalculator( const QString& formulaString, const QString& outputFile, const QString& outputFormat, const QgsRectangle& outputExtent, int nOutputColumns, int nOutputRows, const QVector& rasterEntries ): mFormulaString( formulaString ), mOutputFile( outputFile ), mOutputFormat( outputFormat ), - mOutputRectangle( outputExtent ), mNumOutputColumns( nOutputColumns ), mNumOutputRows( nOutputRows ), mRasterEntries( rasterEntries ) + mOutputRectangle( outputExtent ), mNumOutputColumns( nOutputColumns ), mNumOutputRows( nOutputRows ), mRasterEntries( rasterEntries ) { } @@ -37,7 +39,7 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p ) //prepare search string / tree QString errorString; QgsRasterCalcNode* calcNode = QgsRasterCalcNode::parseRasterCalcString( mFormulaString, errorString ); - if( !calcNode ) + if ( !calcNode ) { //error } @@ -51,19 +53,39 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p ) QVector< GDALDatasetH > mInputDatasets; //raster references and corresponding dataset QVector::const_iterator it = mRasterEntries.constBegin(); - for( ; it != mRasterEntries.constEnd(); ++it ) + for ( ; it != mRasterEntries.constEnd(); ++it ) { - if( !it->raster ) // no raster layer in entry + if ( !it->raster ) // no raster layer in entry { return 2; } GDALDatasetH inputDataset = GDALOpen( it->raster->source().toLocal8Bit().data(), GA_ReadOnly ); - if( inputDataset == NULL ) + if ( inputDataset == NULL ) { return 2; } + + //check if the input dataset is south up or rotated. If yes, use GDALAutoCreateWarpedVRT to create a north up raster + double inputGeoTransform[6]; + if ( GDALGetGeoTransform( inputDataset, inputGeoTransform ) == CE_None + && ( inputGeoTransform[1] < 0.0 + || inputGeoTransform[2] != 0.0 + || inputGeoTransform[4] != 0.0 + || inputGeoTransform[5] > 0.0 ) ) + { + GDALDatasetH vDataset = GDALAutoCreateWarpedVRT( inputDataset, NULL, NULL, GRA_NearestNeighbour, 0.2, NULL ); + mInputDatasets.push_back( vDataset ); + mInputDatasets.push_back( inputDataset ); + inputDataset = vDataset; + } + else + { + mInputDatasets.push_back( inputDataset ); + } + + GDALRasterBandH inputRasterBand = GDALGetRasterBand( inputDataset, it->bandNumber ); - if( inputRasterBand == NULL ) + if ( inputRasterBand == NULL ) { return 2; } @@ -71,14 +93,13 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p ) int nodataSuccess; double nodataValue = GDALGetRasterNoDataValue( inputRasterBand, &nodataSuccess ); - mInputDatasets.push_back( inputDataset ); mInputRasterBands.insert( it->ref, inputRasterBand ); inputScanLineData.insert( it->ref, new QgsRasterMatrix( mNumOutputColumns, 1, new float[mNumOutputColumns], nodataValue ) ); } //open output dataset for writing GDALDriverH outputDriver = openOutputDriver(); - if( outputDriver == NULL ) + if ( outputDriver == NULL ) { return 1; } @@ -90,7 +111,7 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p ) float* resultScanLine = ( float * ) CPLMalloc( sizeof( float ) * mNumOutputColumns ); - if( p ) + if ( p ) { p->setMaximum( mNumOutputRows ); } @@ -98,21 +119,21 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p ) QgsRasterMatrix resultMatrix; //read / write line by line - for( int i = 0; i < mNumOutputRows; ++i ) + for ( int i = 0; i < mNumOutputRows; ++i ) { - if( p ) + if ( p ) { p->setValue( i ); } - if( p && p->wasCanceled() ) + if ( p && p->wasCanceled() ) { break; } //fill buffers QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin(); - for( ; bufferIt != inputScanLineData.end(); ++bufferIt ) + for ( ; bufferIt != inputScanLineData.end(); ++bufferIt ) { double sourceTransformation[6]; GDALRasterBandH sourceRasterBand = mInputRasterBands[bufferIt.key()]; @@ -121,15 +142,15 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p ) readRasterPart( targetGeoTransform, 0, i, mNumOutputColumns, 1, sourceTransformation, sourceRasterBand, bufferIt.value()->data() ); } - if( calcNode->calculate( inputScanLineData, resultMatrix ) ) + if ( calcNode->calculate( inputScanLineData, resultMatrix ) ) { bool resultIsNumber = resultMatrix.isNumber(); float* calcData; - if( resultIsNumber ) //scalar result. Insert number for every pixel + if ( resultIsNumber ) //scalar result. Insert number for every pixel { calcData = new float[mNumOutputColumns]; - for( int j = 0; j < mNumOutputColumns; ++j ) + for ( int j = 0; j < mNumOutputColumns; ++j ) { calcData[j] = resultMatrix.number(); } @@ -140,21 +161,21 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p ) } //replace all matrix nodata values with output nodatas - for( int j = 0; j < mNumOutputColumns; ++j ) + for ( int j = 0; j < mNumOutputColumns; ++j ) { - if( calcData[j] == resultMatrix.nodataValue() ) + if ( calcData[j] == resultMatrix.nodataValue() ) { calcData[j] = outputNodataValue; } } //write scanline to the dataset - if( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None ) + if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None ) { qWarning( "RasterIO error!" ); } - if( resultIsNumber ) + if ( resultIsNumber ) { delete[] calcData; } @@ -162,7 +183,7 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p ) } - if( p ) + if ( p ) { p->setValue( mNumOutputRows ); } @@ -170,19 +191,19 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p ) //close datasets and release memory delete calcNode; QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin(); - for( ; bufferIt != inputScanLineData.end(); ++bufferIt ) + for ( ; bufferIt != inputScanLineData.end(); ++bufferIt ) { delete bufferIt.value(); } inputScanLineData.clear(); QVector< GDALDatasetH >::iterator datasetIt = mInputDatasets.begin(); - for( ; datasetIt != mInputDatasets.end(); ++ datasetIt ) + for ( ; datasetIt != mInputDatasets.end(); ++ datasetIt ) { GDALClose( *datasetIt ); } - if( p && p->wasCanceled() ) + if ( p && p->wasCanceled() ) { //delete the dataset without closing (because it is faster) GDALDeleteDataset( outputDriver, mOutputFile.toLocal8Bit().data() ); @@ -204,13 +225,13 @@ GDALDriverH QgsRasterCalculator::openOutputDriver() //open driver GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() ); - if( outputDriver == NULL ) + if ( outputDriver == NULL ) { return outputDriver; //return NULL, driver does not exist } driverMetadata = GDALGetMetadata( outputDriver, NULL ); - if( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) ) + if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) ) { return NULL; //driver exist, but it does not support the create operation } @@ -223,7 +244,7 @@ GDALDatasetH QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver ) //open output file char **papszOptions = NULL; GDALDatasetH outputDataset = GDALCreate( outputDriver, mOutputFile.toLocal8Bit().data(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions ); - if( outputDataset == NULL ) + if ( outputDataset == NULL ) { return outputDataset; } @@ -239,7 +260,7 @@ GDALDatasetH QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver ) void QgsRasterCalculator::readRasterPart( double* targetGeotransform, int xOffset, int yOffset, int nCols, int nRows, double* sourceTransform, GDALRasterBandH sourceBand, float* rasterBuffer ) { //If dataset transform is the same as the requested transform, do a normal GDAL raster io - if( transformationsEqual( targetGeotransform, sourceTransform ) ) + if ( transformationsEqual( targetGeotransform, sourceTransform ) ) { GDALRasterIO( sourceBand, GF_Read, xOffset, yOffset, nCols, nRows, rasterBuffer, nCols, nRows, GDT_Float32, 0, 0 ); return; @@ -255,10 +276,10 @@ void QgsRasterCalculator::readRasterPart( double* targetGeotransform, int xOffse QgsRectangle intersection = targetRect.intersect( &sourceRect ); //no intersection, fill all the pixels with nodata values - if( intersection.isEmpty() ) + if ( intersection.isEmpty() ) { int nPixels = nCols * nRows; - for( int i = 0; i < nPixels; ++i ) + for ( int i = 0; i < nPixels; ++i ) { rasterBuffer[i] = nodataValue; } @@ -284,17 +305,17 @@ void QgsRasterCalculator::readRasterPart( double* targetGeotransform, int xOffse double targetPixelY = targetGeotransform[3] + targetGeotransform[5] * yOffset + targetGeotransform[5] / 2.0; //coordinates of current target pixel int sourceIndexX, sourceIndexY; //current raster index in source pixels double sx, sy; - for( int i = 0; i < nRows; ++i ) + for ( int i = 0; i < nRows; ++i ) { targetPixelX = targetPixelXMin; - for( int j = 0; j < nCols; ++j ) + for ( int j = 0; j < nCols; ++j ) { sx = ( targetPixelX - sourceRasterXMin ) / sourceTransform[1]; sourceIndexX = sx > 0 ? sx : floor( sx ); sy = ( targetPixelY - sourceRasterYMax ) / sourceTransform[5]; sourceIndexY = sy > 0 ? sy : floor( sy ); - if( sourceIndexX >= 0 && sourceIndexX < nSourcePixelsX - && sourceIndexY >= 0 && sourceIndexY < nSourcePixelsY ) + if ( sourceIndexX >= 0 && sourceIndexX < nSourcePixelsX + && sourceIndexY >= 0 && sourceIndexY < nSourcePixelsY ) { rasterBuffer[j + i*nRows] = sourceRaster[ sourceIndexX + nSourcePixelsX * sourceIndexY ]; } @@ -313,9 +334,9 @@ void QgsRasterCalculator::readRasterPart( double* targetGeotransform, int xOffse bool QgsRasterCalculator::transformationsEqual( double* t1, double* t2 ) const { - for( int i = 0; i < 6; ++i ) + for ( int i = 0; i < 6; ++i ) { - if( !doubleNear( t1[i], t2[i], 0.00001 ) ) + if ( !doubleNear( t1[i], t2[i], 0.00001 ) ) { return false; }