sharpetronics/QGIS
1
0
Fork 0
mirror of https://github.com/qgis/QGIS.git synced 2025-04-16 00:03:12 -04:00
Code Issues Packages Projects Releases Wiki Activity

Use GDALAutoCreateWarpedVRT to handle south-up (or rotated) rasters also for calculator. Fixes bug #3281

Browse Source 
git-svn-id: http://svn.osgeo.org/qgis/trunk@14842 c8812cc2-4d05-0410-92ff-de0c093fc19c
This commit is contained in:
mhugent 2010-12-05 19:22:04 +00:00
parent cfefae12c8
commit 202ec2bca4
1 changed files with 57 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

89
src/analysis/raster/qgsrastercalculator.cpp
View File

@ -22,6 +22,8 @@
#include "cpl_string.h" #include "cpl_string.h"
#include <QProgressDialog> #include <QProgressDialog>
#include "gdalwarper.h"
QgsRasterCalculator::QgsRasterCalculator( const QString& formulaString, const QString& outputFile, const QString& outputFormat, 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 ), 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 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
//prepare search string / tree //prepare search string / tree
QString errorString; QString errorString;
QgsRasterCalcNode* calcNode = QgsRasterCalcNode::parseRasterCalcString( mFormulaString, errorString ); QgsRasterCalcNode* calcNode = QgsRasterCalcNode::parseRasterCalcString( mFormulaString, errorString );
if( !calcNode ) if ( !calcNode )
{ {
//error //error
} }
@ -51,19 +53,39 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
QVector< GDALDatasetH > mInputDatasets; //raster references and corresponding dataset QVector< GDALDatasetH > mInputDatasets; //raster references and corresponding dataset
QVector<QgsRasterCalculatorEntry>::const_iterator it = mRasterEntries.constBegin(); 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; return 2;
} }
GDALDatasetH inputDataset = GDALOpen( it->raster->source().toLocal8Bit().data(), GA_ReadOnly ); GDALDatasetH inputDataset = GDALOpen( it->raster->source().toLocal8Bit().data(), GA_ReadOnly );
if( inputDataset == NULL ) if ( inputDataset == NULL )
{ {
return 2; 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 ); GDALRasterBandH inputRasterBand = GDALGetRasterBand( inputDataset, it->bandNumber );
if( inputRasterBand == NULL ) if ( inputRasterBand == NULL )
{ {
return 2; return 2;
} }
@ -71,14 +93,13 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
int nodataSuccess; int nodataSuccess;
double nodataValue = GDALGetRasterNoDataValue( inputRasterBand, &nodataSuccess ); double nodataValue = GDALGetRasterNoDataValue( inputRasterBand, &nodataSuccess );
mInputDatasets.push_back( inputDataset );
mInputRasterBands.insert( it->ref, inputRasterBand ); mInputRasterBands.insert( it->ref, inputRasterBand );
inputScanLineData.insert( it->ref, new QgsRasterMatrix( mNumOutputColumns, 1, new float[mNumOutputColumns], nodataValue ) ); inputScanLineData.insert( it->ref, new QgsRasterMatrix( mNumOutputColumns, 1, new float[mNumOutputColumns], nodataValue ) );
} }
//open output dataset for writing //open output dataset for writing
GDALDriverH outputDriver = openOutputDriver(); GDALDriverH outputDriver = openOutputDriver();
if( outputDriver == NULL ) if ( outputDriver == NULL )
{ {
return 1; return 1;
} }
@ -90,7 +111,7 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
float* resultScanLine = ( float * ) CPLMalloc( sizeof( float ) * mNumOutputColumns ); float* resultScanLine = ( float * ) CPLMalloc( sizeof( float ) * mNumOutputColumns );
if( p ) if ( p )
{ {
p->setMaximum( mNumOutputRows ); p->setMaximum( mNumOutputRows );
} }
@ -98,21 +119,21 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
QgsRasterMatrix resultMatrix; QgsRasterMatrix resultMatrix;
//read / write line by line //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 ); p->setValue( i );
} }
if( p && p->wasCanceled() ) if ( p && p->wasCanceled() )
{ {
break; break;
} }
//fill buffers //fill buffers
QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin(); QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin();
for( ; bufferIt != inputScanLineData.end(); ++bufferIt ) for ( ; bufferIt != inputScanLineData.end(); ++bufferIt )
{ {
double sourceTransformation[6]; double sourceTransformation[6];
GDALRasterBandH sourceRasterBand = mInputRasterBands[bufferIt.key()]; 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() ); 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(); bool resultIsNumber = resultMatrix.isNumber();
float* calcData; float* calcData;
if( resultIsNumber ) //scalar result. Insert number for every pixel if ( resultIsNumber ) //scalar result. Insert number for every pixel
{ {
calcData = new float[mNumOutputColumns]; calcData = new float[mNumOutputColumns];
for( int j = 0; j < mNumOutputColumns; ++j ) for ( int j = 0; j < mNumOutputColumns; ++j )
{ {
calcData[j] = resultMatrix.number(); calcData[j] = resultMatrix.number();
} }
@ -140,21 +161,21 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
} }
//replace all matrix nodata values with output nodatas //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; calcData[j] = outputNodataValue;
} }
} }
//write scanline to the dataset //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!" ); qWarning( "RasterIO error!" );
} }
if( resultIsNumber ) if ( resultIsNumber )
{ {
delete[] calcData; delete[] calcData;
} }
@ -162,7 +183,7 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
} }
if( p ) if ( p )
{ {
p->setValue( mNumOutputRows ); p->setValue( mNumOutputRows );
} }
@ -170,19 +191,19 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
//close datasets and release memory //close datasets and release memory
delete calcNode; delete calcNode;
QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin(); QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin();
for( ; bufferIt != inputScanLineData.end(); ++bufferIt ) for ( ; bufferIt != inputScanLineData.end(); ++bufferIt )
{ {
delete bufferIt.value(); delete bufferIt.value();
} }
inputScanLineData.clear(); inputScanLineData.clear();
QVector< GDALDatasetH >::iterator datasetIt = mInputDatasets.begin(); QVector< GDALDatasetH >::iterator datasetIt = mInputDatasets.begin();
for( ; datasetIt != mInputDatasets.end(); ++ datasetIt ) for ( ; datasetIt != mInputDatasets.end(); ++ datasetIt )
{ {
GDALClose( *datasetIt ); GDALClose( *datasetIt );
} }
if( p && p->wasCanceled() ) if ( p && p->wasCanceled() )
{ {
//delete the dataset without closing (because it is faster) //delete the dataset without closing (because it is faster)
GDALDeleteDataset( outputDriver, mOutputFile.toLocal8Bit().data() ); GDALDeleteDataset( outputDriver, mOutputFile.toLocal8Bit().data() );
@ -204,13 +225,13 @@ GDALDriverH QgsRasterCalculator::openOutputDriver()
//open driver //open driver
GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() ); GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
if( outputDriver == NULL ) if ( outputDriver == NULL )
{ {
return outputDriver; //return NULL, driver does not exist return outputDriver; //return NULL, driver does not exist
} }
driverMetadata = GDALGetMetadata( outputDriver, NULL ); 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 return NULL; //driver exist, but it does not support the create operation
} }
@ -223,7 +244,7 @@ GDALDatasetH QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver )
//open output file //open output file
char **papszOptions = NULL; char **papszOptions = NULL;
GDALDatasetH outputDataset = GDALCreate( outputDriver, mOutputFile.toLocal8Bit().data(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions ); GDALDatasetH outputDataset = GDALCreate( outputDriver, mOutputFile.toLocal8Bit().data(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions );
if( outputDataset == NULL ) if ( outputDataset == NULL )
{ {
return outputDataset; 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 ) 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 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 ); GDALRasterIO( sourceBand, GF_Read, xOffset, yOffset, nCols, nRows, rasterBuffer, nCols, nRows, GDT_Float32, 0, 0 );
return; return;
@ -255,10 +276,10 @@ void QgsRasterCalculator::readRasterPart( double* targetGeotransform, int xOffse
QgsRectangle intersection = targetRect.intersect( &sourceRect ); QgsRectangle intersection = targetRect.intersect( &sourceRect );
//no intersection, fill all the pixels with nodata values //no intersection, fill all the pixels with nodata values
if( intersection.isEmpty() ) if ( intersection.isEmpty() )
{ {
int nPixels = nCols * nRows; int nPixels = nCols * nRows;
for( int i = 0; i < nPixels; ++i ) for ( int i = 0; i < nPixels; ++i )
{ {
rasterBuffer[i] = nodataValue; rasterBuffer[i] = nodataValue;
} }
@ -284,16 +305,16 @@ void QgsRasterCalculator::readRasterPart( double* targetGeotransform, int xOffse
double targetPixelY = targetGeotransform[3] + targetGeotransform[5] * yOffset + targetGeotransform[5] / 2.0; //coordinates of current target pixel 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 int sourceIndexX, sourceIndexY; //current raster index in source pixels
double sx, sy; double sx, sy;
for( int i = 0; i < nRows; ++i ) for ( int i = 0; i < nRows; ++i )
{ {
targetPixelX = targetPixelXMin; targetPixelX = targetPixelXMin;
for( int j = 0; j < nCols; ++j ) for ( int j = 0; j < nCols; ++j )
{ {
sx = ( targetPixelX - sourceRasterXMin ) / sourceTransform[1]; sx = ( targetPixelX - sourceRasterXMin ) / sourceTransform[1];
sourceIndexX = sx > 0 ? sx : floor( sx ); sourceIndexX = sx > 0 ? sx : floor( sx );
sy = ( targetPixelY - sourceRasterYMax ) / sourceTransform[5]; sy = ( targetPixelY - sourceRasterYMax ) / sourceTransform[5];
sourceIndexY = sy > 0 ? sy : floor( sy ); sourceIndexY = sy > 0 ? sy : floor( sy );
if( sourceIndexX >= 0 && sourceIndexX < nSourcePixelsX if ( sourceIndexX >= 0 && sourceIndexX < nSourcePixelsX
&& sourceIndexY >= 0 && sourceIndexY < nSourcePixelsY ) && sourceIndexY >= 0 && sourceIndexY < nSourcePixelsY )
{ {
rasterBuffer[j + i*nRows] = sourceRaster[ sourceIndexX + nSourcePixelsX * sourceIndexY ]; 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 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; return false;
} }