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

Scanline implementation

Browse Source
This commit is contained in:
Alessandro Pasotti 2018-11-27 21:28:35 +01:00
parent 10517c8932
commit e329fb5b97
5 changed files with 228 additions and 89 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
python/analysis/auto_generated/raster/qgsrastercalcnode.sip.in
View File

@ -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:

10
src/analysis/raster/qgsrastercalcnode.cpp
View File

@ -16,6 +16,7 @@
#include "qgsrasterblock.h"
#include "qgsrastermatrix.h"
#include <cfloat>
#include <QtDebug>
QgsRasterCalcNode::QgsRasterCalcNode( double number )
: mNumber( number )
@ -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 );
}
}
@ -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 )

2
src/analysis/raster/qgsrastercalcnode.h
View File

@ -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;

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

@ -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
@ -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() )
{
break;
}
if ( calcNode->calculate( inputBlocks, resultMatrix, i ) )
{
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;
}
}
if ( feedback )
else
{
feedback->setProgress( 100.0 );
QMap< QString, QgsRasterBlock * > inputBlocks;
QVector<QgsRasterCalculatorEntry>::const_iterator it = mRasterEntries.constBegin();
for ( ; it != mRasterEntries.constEnd(); ++it )
{
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() );
}
QgsRasterMatrix resultMatrix;
resultMatrix.setNodataValue( outputNodataValue );
//read / write line by line
for ( int i = 0; i < mNumOutputRows; ++i )
{
if ( feedback )
{
feedback->setProgress( 100.0 * static_cast< double >( i ) / mNumOutputRows );
}
if ( feedback && feedback->isCanceled() )
{
break;
}
if ( calcNode->calculate( inputBlocks, resultMatrix, i ) )
{
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;
}
}
if ( feedback )
{
feedback->setProgress( 100.0 );
}
//close datasets and release memory
calcNode.reset();
qDeleteAll( inputBlocks );
inputBlocks.clear();
}
//close datasets and release memory
calcNode.reset();
qDeleteAll( inputBlocks );
inputBlocks.clear();
if ( feedback && feedback->isCanceled() )
{

31
tests/src/analysis/testqgsrastercalculator.cpp
View File

@ -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 );
@ -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"