Skip to content

Commit

Permalink
Scanline implementation
Browse files Browse the repository at this point in the history
  • Loading branch information
@elpaso
elpaso committed Nov 28, 2018
1 parent 10517c8 commit e329fb5
Show file tree
Hide file tree
Showing 5 changed files with 217 additions and 78 deletions.
16 changes: 16 additions & 0 deletions python/analysis/auto_generated/raster/qgsrastercalcnode.sip.in
View file
Expand Up @@ -74,6 +74,22 @@ QgsRasterCalcNode cannot be copied
void setRight( QgsRasterCalcNode *right );


QString toString( bool cStyle = false ) const;
%Docstring
Returns a string representation of the expression

:param cStyle: if true operators will follow C syntax

.. versionadded:: 3.6
%End

QList<const QgsRasterCalcNode *> findNodes( const QgsRasterCalcNode::Type type ) const;
%Docstring
Populates a list of nodes of a specific type.

.. versionadded:: 3.6
%End

static QgsRasterCalcNode *parseRasterCalcString( const QString &str, QString &parserErrorMsg ) /Factory/;

private:
Expand Down
10 changes: 7 additions & 3 deletions src/analysis/raster/qgsrastercalcnode.cpp
View file
Expand Up @@ -16,6 +16,7 @@
#include "qgsrasterblock.h"
#include "qgsrastermatrix.h"
#include <cfloat>
#include <QtDebug>

QgsRasterCalcNode::QgsRasterCalcNode( double number )
: mNumber( number )
Expand Down Expand Up @@ -78,6 +79,7 @@ bool QgsRasterCalcNode::calculate( QMap<QString, QgsRasterBlock * > &rasterData,
{
for ( int dataCol = 0; dataCol < nCols; ++dataCol )
{
//qDebug() << "Reading value:" << dataRow << dataCol << ( *it )->value( dataRow, dataCol );
data[ dataCol + nCols * outRow] = ( *it )->isNoData( dataRow, dataCol ) ? result.nodataValue() : ( *it )->value( dataRow, dataCol );
}
}
Expand Down Expand Up @@ -308,14 +310,16 @@ QString QgsRasterCalcNode::toString( bool cStyle ) const
return result;
}

void QgsRasterCalcNode::findNodes( const QgsRasterCalcNode::Type type, QList<const QgsRasterCalcNode *> &nodeList ) const
QList<const QgsRasterCalcNode *> QgsRasterCalcNode::findNodes( const QgsRasterCalcNode::Type type ) const
{
QList<const QgsRasterCalcNode *> nodeList;
if ( mType == type )
nodeList.push_back( this );
if ( mLeft )
mLeft->findNodes( type, nodeList );
nodeList.append( mLeft->findNodes( type ) );
if ( mRight )
mRight->findNodes( type, nodeList );
nodeList.append( mRight->findNodes( type ) );
return nodeList;
}

QgsRasterCalcNode *QgsRasterCalcNode::parseRasterCalcString( const QString &str, QString &parserErrorMsg )
Expand Down
2 changes: 1 addition & 1 deletion src/analysis/raster/qgsrastercalcnode.h
View file
Expand Up @@ -117,7 +117,7 @@ class ANALYSIS_EXPORT QgsRasterCalcNode
* Populates a list of nodes of a specific type.
* \since QGIS 3.6
*/
void findNodes( const QgsRasterCalcNode::Type type, QList<const QgsRasterCalcNode *> &nodeList ) const;
QList<const QgsRasterCalcNode *> findNodes( const QgsRasterCalcNode::Type type ) const;

static QgsRasterCalcNode *parseRasterCalcString( const QString &str, QString &parserErrorMsg ) SIP_FACTORY;

Expand Down
236 changes: 172 additions & 64 deletions src/analysis/raster/qgsrastercalculator.cpp
View file
Expand Up @@ -70,53 +70,19 @@ QgsRasterCalculator::Result QgsRasterCalculator::processCalculation( QgsFeedback
return ParserError;
}

QMap< QString, QgsRasterBlock * > inputBlocks;
QVector<QgsRasterCalculatorEntry>::const_iterator it = mRasterEntries.constBegin();
for ( ; it != mRasterEntries.constEnd(); ++it )
// Check input layers and bands
for ( const auto &entry : qgis::as_const( mRasterEntries ) )
{
if ( !it->raster ) // no raster layer in entry
if ( !entry.raster ) // no raster layer in entry
{
mLastError = QObject::tr( "No raster layer for entry %1" ).arg( it->ref );
qDeleteAll( inputBlocks );
mLastError = QObject::tr( "No raster layer for entry %1" ).arg( entry.ref );
return InputLayerError;
}

if ( it->bandNumber <= 0 || it->bandNumber > it->raster->bandCount() )
if ( entry.bandNumber <= 0 || entry.bandNumber > entry.raster->bandCount() )
{
mLastError = QObject::tr( "Band number %1 is not valid for entry %2" ).arg( it->bandNumber ).arg( it->ref );
qDeleteAll( inputBlocks );
mLastError = QObject::tr( "Band number %1 is not valid for entry %2" ).arg( entry.bandNumber ).arg( entry.ref );
return BandError;
}

std::unique_ptr< QgsRasterBlock > block;
// if crs transform needed
if ( it->raster->crs() != mOutputCrs )
{
QgsRasterProjector proj;
proj.setCrs( it->raster->crs(), mOutputCrs );
proj.setInput( it->raster->dataProvider() );
proj.setPrecision( QgsRasterProjector::Exact );

QgsRasterBlockFeedback *rasterBlockFeedback = new QgsRasterBlockFeedback();
QObject::connect( feedback, &QgsFeedback::canceled, rasterBlockFeedback, &QgsRasterBlockFeedback::cancel );
block.reset( proj.block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows, rasterBlockFeedback ) );
if ( rasterBlockFeedback->isCanceled() )
{
qDeleteAll( inputBlocks );
return Canceled;
}
}
else
{
block.reset( it->raster->dataProvider()->block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows ) );
}
if ( block->isEmpty() )
{
mLastError = QObject::tr( "Could not allocate required memory for %1" ).arg( it->ref );
qDeleteAll( inputBlocks );
return MemoryError;
}
inputBlocks.insert( it->ref, block.release() );
}

//open output dataset for writing
Expand All @@ -140,52 +106,194 @@ QgsRasterCalculator::Result QgsRasterCalculator::processCalculation( QgsFeedback
float outputNodataValue = -FLT_MAX;
GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );

QgsRasterMatrix resultMatrix;
resultMatrix.setNodataValue( outputNodataValue );

//read / write line by line
for ( int i = 0; i < mNumOutputRows; ++i )
// If we need to read the raster as a whole
bool requiresMatrix = ! calcNode->findNodes( QgsRasterCalcNode::Type::tMatrix ).isEmpty();


if ( ! requiresMatrix )
{
// Map of raster names -> blocks
std::map<QString, std::unique_ptr<QgsRasterBlock>> inputBlocks;
std::map<QString, QgsRasterCalculatorEntry> uniqueRasterEntries;
for ( const auto &r : calcNode->findNodes( QgsRasterCalcNode::Type::tRasterRef ) )
{
QString layerRef( r->toString().remove( 0, 1 ) );
layerRef.chop( 1 );
if ( ! inputBlocks.count( layerRef ) )
{
for ( const auto &ref : mRasterEntries )
{
if ( ref.ref == layerRef )
{
uniqueRasterEntries[layerRef] = ref;
inputBlocks[layerRef ] = qgis::make_unique<QgsRasterBlock>();
}
}
}
}

//read / write line by line
QMap<QString, QgsRasterBlock * > _rasterData;
// Cast to float
std::vector<float> castedResult;
castedResult.reserve( static_cast<size_t>( mNumOutputColumns ) );
auto rowHeight = mOutputRectangle.height() / mNumOutputRows;
for ( size_t row = 0; row < mNumOutputRows; ++row )
{
if ( feedback )
{
feedback->setProgress( 100.0 * static_cast< double >( row ) / mNumOutputRows );
}

if ( feedback && feedback->isCanceled() )
{
break;
}

// Read one row
QgsRectangle rect( mOutputRectangle );
rect.setYMaximum( rect.yMaximum() - rowHeight * row );
rect.setYMinimum( rect.yMaximum() - rowHeight );

// Read blocks
for ( auto &layerRef : inputBlocks )
{
QgsRasterCalculatorEntry ref = uniqueRasterEntries[layerRef.first];
if ( uniqueRasterEntries[layerRef.first].raster->crs() != mOutputCrs )
{
QgsRasterProjector proj;
proj.setCrs( ref.raster->crs(), mOutputCrs );
proj.setInput( ref.raster->dataProvider() );
proj.setPrecision( QgsRasterProjector::Exact );
layerRef.second.reset( proj.block( ref.bandNumber, rect, mNumOutputColumns, 1 ) );
}
else
{
inputBlocks[layerRef.first].reset( ref.raster->dataProvider()->block( ref.bandNumber, rect, mNumOutputColumns, 1 ) );
}
}

QgsRasterMatrix resultMatrix;
resultMatrix.setNodataValue( outputNodataValue );

_rasterData.clear();
for ( const auto &layerRef : inputBlocks )
{
_rasterData.insert( layerRef.first, inputBlocks[layerRef.first].get() );
//for ( int i = 0; i < mNumOutputColumns; i++ )
// qDebug() << "Input: " << row << i << " = " << inputBlocks[layerRef.first]->value(0, i);
}

if ( calcNode->calculate( _rasterData, resultMatrix, 0 ) )
{
//write scanline to the dataset
for ( size_t i = 0; i < static_cast<size_t>( mNumOutputColumns ); i++ )
{
castedResult[i] = static_cast<float>( resultMatrix.data()[i] );
// qDebug() << "Calculated: " << row << i << " = " << castedResult[i];
}
if ( GDALRasterIO( outputRasterBand, GF_Write, 0, row, mNumOutputColumns, 1, castedResult.data(), mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
{
QgsDebugMsg( QStringLiteral( "RasterIO error!" ) );
}
}
}

if ( feedback )
{
feedback->setProgress( 100.0 * static_cast< double >( i ) / mNumOutputRows );
feedback->setProgress( 100.0 );
}

if ( feedback && feedback->isCanceled() )

}
else
{
QMap< QString, QgsRasterBlock * > inputBlocks;
QVector<QgsRasterCalculatorEntry>::const_iterator it = mRasterEntries.constBegin();
for ( ; it != mRasterEntries.constEnd(); ++it )
{
break;

std::unique_ptr< QgsRasterBlock > block;
// if crs transform needed
if ( it->raster->crs() != mOutputCrs )
{
QgsRasterProjector proj;
proj.setCrs( it->raster->crs(), mOutputCrs );
proj.setInput( it->raster->dataProvider() );
proj.setPrecision( QgsRasterProjector::Exact );

QgsRasterBlockFeedback *rasterBlockFeedback = new QgsRasterBlockFeedback();
QObject::connect( feedback, &QgsFeedback::canceled, rasterBlockFeedback, &QgsRasterBlockFeedback::cancel );
block.reset( proj.block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows, rasterBlockFeedback ) );
if ( rasterBlockFeedback->isCanceled() )
{
qDeleteAll( inputBlocks );
return Canceled;
}
}
else
{
block.reset( it->raster->dataProvider()->block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows ) );
}
if ( block->isEmpty() )
{
mLastError = QObject::tr( "Could not allocate required memory for %1" ).arg( it->ref );
qDeleteAll( inputBlocks );
return MemoryError;
}
inputBlocks.insert( it->ref, block.release() );
}

if ( calcNode->calculate( inputBlocks, resultMatrix, i ) )
QgsRasterMatrix resultMatrix;
resultMatrix.setNodataValue( outputNodataValue );

//read / write line by line
for ( int i = 0; i < mNumOutputRows; ++i )
{
bool resultIsNumber = resultMatrix.isNumber();
float *calcData = new float[mNumOutputColumns];
if ( feedback )
{
feedback->setProgress( 100.0 * static_cast< double >( i ) / mNumOutputRows );
}

for ( int j = 0; j < mNumOutputColumns; ++j )
if ( feedback && feedback->isCanceled() )
{
calcData[j] = ( float )( resultIsNumber ? resultMatrix.number() : resultMatrix.data()[j] );
break;
}

//write scanline to the dataset
if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
if ( calcNode->calculate( inputBlocks, resultMatrix, i ) )
{
QgsDebugMsg( QStringLiteral( "RasterIO error!" ) );
bool resultIsNumber = resultMatrix.isNumber();
float *calcData = new float[mNumOutputColumns];

for ( int j = 0; j < mNumOutputColumns; ++j )
{
calcData[j] = ( float )( resultIsNumber ? resultMatrix.number() : resultMatrix.data()[j] );
}

//write scanline to the dataset
if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
{
QgsDebugMsg( QStringLiteral( "RasterIO error!" ) );
}

delete[] calcData;
}

delete[] calcData;
}

}
if ( feedback )
{
feedback->setProgress( 100.0 );
}

//close datasets and release memory
calcNode.reset();
qDeleteAll( inputBlocks );
inputBlocks.clear();

if ( feedback )
{
feedback->setProgress( 100.0 );
}

//close datasets and release memory
calcNode.reset();
qDeleteAll( inputBlocks );
inputBlocks.clear();

if ( feedback && feedback->isCanceled() )
{
Expand Down
31 changes: 21 additions & 10 deletions tests/src/analysis/testqgsrastercalculator.cpp
View file
Expand Up @@ -538,19 +538,11 @@ void TestQgsRasterCalculator::findNodes()
std::unique_ptr< QgsRasterCalcNode > calcNode;

auto _test =
[ & ]( QString exp, const QgsRasterCalcNode::Type type ) -> QStringList
[ & ]( QString exp, const QgsRasterCalcNode::Type type ) -> QList<const QgsRasterCalcNode *>
{
QList<const QgsRasterCalcNode *> nodeList;
QString error;
calcNode.reset( QgsRasterCalcNode::parseRasterCalcString( exp, error ) );
if ( error.isEmpty() )
calcNode->findNodes( type, nodeList );
QStringList result;
for ( const auto &n : nodeList )
{
result.push_back( n->toString( true ) );
}
return result;
return calcNode->findNodes( type );
};

QCOMPARE( _test( QStringLiteral( "atan(\"raster@1\") * cos( 3 + 2 )" ), QgsRasterCalcNode::Type::tOperator ).length(), 4 );
Expand Down Expand Up @@ -642,5 +634,24 @@ void TestQgsRasterCalculator::errors( )
QVERIFY( rc.lastError().isEmpty() );
}

void TestQgsRasterCalculator::toString()
{
auto _test = [ ]( QString exp, bool cStyle ) -> QString
{
QString error;
std::unique_ptr< QgsRasterCalcNode > calcNode( QgsRasterCalcNode::parseRasterCalcString( exp, error ) );
if ( ! error.isEmpty() )
return error;
return calcNode->toString( cStyle );
};
QCOMPARE( _test( QStringLiteral( "\"raster@1\" + 2" ), false ), QString( "\"raster@1\" + 2" ) );
QCOMPARE( _test( QStringLiteral( "\"raster@1\" + 2" ), true ), QString( "\"raster@1\" + 2" ) );
QCOMPARE( _test( QStringLiteral( "\"raster@1\" ^ 3 + 2" ), false ), QString( "\"raster@1\"^3 + 2" ) );
QCOMPARE( _test( QStringLiteral( "\"raster@1\" ^ 3 + 2" ), true ), QString( "pow( \"raster@1\", 3 ) + 2" ) );
QCOMPARE( _test( QStringLiteral( "atan(\"raster@1\") * cos( 3 + 2 )" ), false ), QString( "atan( \"raster@1\" ) * cos( 3 + 2 )" ) );
QCOMPARE( _test( QStringLiteral( "atan(\"raster@1\") * cos( 3 + 2 )" ), true ), QString( "atan( \"raster@1\" ) * cos( 3 + 2 )" ) );
}


QGSTEST_MAIN( TestQgsRasterCalculator )
#include "testqgsrastercalculator.moc"

0 comments on commit e329fb5

Please sign in to comment.