[QGIS Commit] r14842 - trunk/qgis/src/analysis/raster
svn_qgis at osgeo.org
svn_qgis at osgeo.org
Sun Dec 5 14:22:04 EST 2010
Author: mhugent
Date: 2010-12-05 11:22:04 -0800 (Sun, 05 Dec 2010)
New Revision: 14842
Modified:
trunk/qgis/src/analysis/raster/qgsrastercalculator.cpp
Log:
Use GDALAutoCreateWarpedVRT to handle south-up (or rotated) rasters also for calculator. Fixes bug #3281
Modified: trunk/qgis/src/analysis/raster/qgsrastercalculator.cpp
===================================================================
--- trunk/qgis/src/analysis/raster/qgsrastercalculator.cpp 2010-12-05 16:12:51 UTC (rev 14841)
+++ trunk/qgis/src/analysis/raster/qgsrastercalculator.cpp 2010-12-05 19:22:04 UTC (rev 14842)
@@ -22,9 +22,11 @@
#include "cpl_string.h"
#include <QProgressDialog>
+#include "gdalwarper.h"
+
QgsRasterCalculator::QgsRasterCalculator( const QString& formulaString, const QString& outputFile, const QString& outputFormat,
const QgsRectangle& outputExtent, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry>& 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 @@
//prepare search string / tree
QString errorString;
QgsRasterCalcNode* calcNode = QgsRasterCalcNode::parseRasterCalcString( mFormulaString, errorString );
- if( !calcNode )
+ if ( !calcNode )
{
//error
}
@@ -51,19 +53,39 @@
QVector< GDALDatasetH > mInputDatasets; //raster references and corresponding dataset
QVector<QgsRasterCalculatorEntry>::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 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 @@
float* resultScanLine = ( float * ) CPLMalloc( sizeof( float ) * mNumOutputColumns );
- if( p )
+ if ( p )
{
p->setMaximum( mNumOutputRows );
}
@@ -98,21 +119,21 @@
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 @@
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 @@
}
//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 @@
}
- if( p )
+ if ( p )
{
p->setValue( mNumOutputRows );
}
@@ -170,19 +191,19 @@
//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 @@
//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 @@
//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 @@
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 @@
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 @@
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 @@
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;
}
More information about the QGIS-commit
mailing list