[rastercalc] Rework raster calculator to use QGIS raster classes

...rather than reading input layers directly through GDAL.
Benefits include more robust handling of nodata/data type conversions,
less code duplication, also being able to take advantage of features
in QGIS raster code like handling gain/offset in rasters. (fix #12450)

Also, add a choice of output projection to the raster calculator.
Previously the output CRS would be taken from the first raster, with
no guarantees that the output extent matched the output CRS. This
resulted in empty/misplaced rasters. (fix #3649)

python/analysis/raster/qgsrastercalcnode.sip

@@ -54,7 +54,21 @@ class QgsRasterCalcNode
     void setRight( QgsRasterCalcNode* right );
 
     /**Calculates result (might be real matrix or single number)*/
-    bool calculate( QMap<QString, QgsRasterMatrix*>& rasterData, QgsRasterMatrix& result ) const;
+   // bool calculate( QMap<QString, QgsRasterCalculateInputMatrix*>& rasterData, QgsRasterMatrix& result ) const;
+
+    /**Calculates result of raster calculation (might be real matrix or single number).
+     * @param rasterData input raster data references, map of raster name to raster data block
+     * @param result destination raster matrix for calculation results
+     * @param row optional row number to calculate for calculating result by rows, or -1 to
+     * calculate entire result
+     * @note added in QGIS 2.10
+     * @note not available in Python bindings
+     */
+    //bool calculate( QMap<QString, QgsRasterBlock* >& rasterData, QgsRasterMatrix& result, int row = -1 ) const;
+
+    /**@deprecated use method which accepts QgsRasterBlocks instead
+     */
+    // bool calculate( QMap<QString, QgsRasterMatrix*>& rasterData, QgsRasterMatrix& result ) const /Deprecated/;
 
     static QgsRasterCalcNode* parseRasterCalcString( const QString& str, QString& parserErrorMsg ) /Factory/;
 };

python/analysis/raster/qgsrastercalculator.sip

@@ -17,8 +17,33 @@ class QgsRasterCalculator
 %End
 
   public:
+
+    /** QgsRasterCalculator constructor.
+     * @param formulaString formula for raster calculation
+     * @param outputFile output file path
+     * @param outputFormat output file format
+     * @param outputExtent output extent. CRS for output is taken from first entry in rasterEntries.
+     * @param nOutputColumns number of columns in output raster
+     * @param nOutputRows number of rows in output raster
+     * @param rasterEntries list of referenced raster layers
+     */
     QgsRasterCalculator( const QString& formulaString, const QString& outputFile, const QString& outputFormat,
                          const QgsRectangle& outputExtent, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry>& rasterEntries );
+
+    /** QgsRasterCalculator constructor.
+     * @param formulaString formula for raster calculation
+     * @param outputFile output file path
+     * @param outputFormat output file format
+     * @param outputExtent output extent, CRS is specifed by outputCrs parameter
+     * @param outputCrs destination CRS for output raster
+     * @param nOutputColumns number of columns in output raster
+     * @param nOutputRows number of rows in output raster
+     * @param rasterEntries list of referenced raster layers
+     * @note added in QGIS 2.10
+     */
+    QgsRasterCalculator( const QString& formulaString, const QString& outputFile, const QString& outputFormat,
+                         const QgsRectangle& outputExtent, const QgsCoordinateReferenceSystem& outputCrs, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry>& rasterEntries );
+
     ~QgsRasterCalculator();
 
     /**Starts the calculation and writes new raster

python/core/raster/qgsrasterblock.sip

@@ -266,5 +266,17 @@ class QgsRasterBlock
      * @return the rectangle covered by sub extent
      */
     static QRect subRect( const QgsRectangle &theExtent, int theWidth, int theHeight, const QgsRectangle &theSubExtent );
+
+    /** Returns the width (number of columns) of the raster block.
+     * @see height
+     * @note added in QGIS 2.10
+     */
+    int width() const;
+
+    /** Returns the height (number of rows) of the raster block.
+     * @see width
+     * @note added in QGIS 2.10
+     */
+    int height() const;
 };
 

src/analysis/raster/qgsrastercalcnode.cpp

@@ -13,6 +13,7 @@
  *                                                                         *
  ***************************************************************************/
 #include "qgsrastercalcnode.h"
+#include "qgsrasterblock.h"
 #include <cfloat>
 
 QgsRasterCalcNode::QgsRasterCalcNode()
@@ -82,26 +83,58 @@ QgsRasterCalcNode::~QgsRasterCalcNode()
 }
 
 bool QgsRasterCalcNode::calculate( QMap<QString, QgsRasterMatrix*>& rasterData, QgsRasterMatrix& result ) const
+{
+  //deprecated method
+  //convert QgsRasterMatrix to QgsRasterBlock and call replacement method
+  QMap<QString, QgsRasterBlock* > rasterBlockData;
+  QMap<QString, QgsRasterMatrix*>::const_iterator it = rasterData.constBegin();
+  for ( ; it != rasterData.constEnd(); ++it )
+  {
+    QgsRasterBlock* block = new QgsRasterBlock( QGis::Float32, it.value()->nColumns(), it.value()->nRows(), it.value()->nodataValue() );
+    for ( int row = 0; row < it.value()->nRows(); ++row )
+    {
+      for ( int col = 0; col < it.value()->nColumns(); ++col )
+      {
+        block->setValue( row, col, it.value()->data()[ row * it.value()->nColumns() + col ] );
+      }
+    }
+    rasterBlockData.insert( it.key(), block );
+  }
+
+  return calculate( rasterBlockData, result );
+}
+
+bool QgsRasterCalcNode::calculate( QMap<QString, QgsRasterBlock* >& rasterData, QgsRasterMatrix& result, int row ) const
 {
   //if type is raster ref: return a copy of the corresponding matrix
 
   //if type is operator, call the proper matrix operations
   if ( mType == tRasterRef )
   {
-    QMap<QString, QgsRasterMatrix*>::iterator it = rasterData.find( mRasterName );
+    QMap<QString, QgsRasterBlock*>::iterator it = rasterData.find( mRasterName );
     if ( it == rasterData.end() )
     {
       return false;
     }
 
-    int nEntries = ( *it )->nColumns() * ( *it )->nRows();
+    int nRows = ( row >= 0 ? 1 : ( *it )->height() );
+    int startRow = ( row >= 0 ? row : 0 );
+    int endRow = startRow + nRows;
+    int nCols = ( *it )->width();
+    int nEntries = nCols * nRows;
     double* data = new double[nEntries];
 
-    for ( int i = 0; i < nEntries; ++i )
+    //convert input raster values to double, also convert input no data to result no data
+
+    int outRow = 0;
+    for ( int dataRow = startRow; dataRow < endRow ; ++dataRow, ++outRow )
     {
-      data[i] = ( *it )->data()[i] == ( *it )->nodataValue() ? result.nodataValue() : ( *it )->data()[i];
+      for ( int dataCol = 0; dataCol < nCols; ++dataCol )
+      {
+        data[ dataCol + nCols * outRow] = ( *it )->isNoData( dataRow , dataCol ) ? result.nodataValue() : ( *it )->value( dataRow, dataCol );
+      }
     }
-    result.setData(( *it )->nColumns(), ( *it )->nRows(), data, result.nodataValue() );
+    result.setData( nCols, nRows, data, result.nodataValue() );
     return true;
   }
   else if ( mType == tOperator )
@@ -110,11 +143,11 @@ bool QgsRasterCalcNode::calculate( QMap<QString, QgsRasterMatrix*>& rasterData,
     leftMatrix.setNodataValue( result.nodataValue() );
     rightMatrix.setNodataValue( result.nodataValue() );
 
-    if ( !mLeft || !mLeft->calculate( rasterData, leftMatrix ) )
+    if ( !mLeft || !mLeft->calculate( rasterData, leftMatrix, row ) )
     {
       return false;
     }
-    if ( mRight && !mRight->calculate( rasterData, rightMatrix ) )
+    if ( mRight && !mRight->calculate( rasterData, rightMatrix, row ) )
     {
       return false;
     }

src/analysis/raster/qgsrastercalcnode.h

@@ -23,6 +23,8 @@
 #include <QMap>
 #include <QString>
 
+class QgsRasterBlock;
+
 class ANALYSIS_EXPORT QgsRasterCalcNode
 {
   public:
@@ -77,8 +79,19 @@ class ANALYSIS_EXPORT QgsRasterCalcNode
     void setLeft( QgsRasterCalcNode* left ) { delete mLeft; mLeft = left; }
     void setRight( QgsRasterCalcNode* right ) { delete mRight; mRight = right; }
 
-    /**Calculates result (might be real matrix or single number)*/
-    bool calculate( QMap<QString, QgsRasterMatrix*>& rasterData, QgsRasterMatrix& result ) const;
+    /**Calculates result of raster calculation (might be real matrix or single number).
+     * @param rasterData input raster data references, map of raster name to raster data block
+     * @param result destination raster matrix for calculation results
+     * @param row optional row number to calculate for calculating result by rows, or -1 to
+     * calculate entire result
+     * @note added in QGIS 2.10
+     * @note not available in Python bindings
+     */
+    bool calculate( QMap<QString, QgsRasterBlock* >& rasterData, QgsRasterMatrix& result, int row = -1 ) const;
+
+    /**@deprecated use method which accepts QgsRasterBlocks instead
+     */
+    Q_DECL_DEPRECATED bool calculate( QMap<QString, QgsRasterMatrix*>& rasterData, QgsRasterMatrix& result ) const;
 
     static QgsRasterCalcNode* parseRasterCalcString( const QString& str, QString& parserErrorMsg );
 

src/analysis/raster/qgsrastercalculator.cpp

@@ -45,8 +45,24 @@ QgsRasterCalculator::QgsRasterCalculator( const QString& formulaString, const QS
     , mNumOutputRows( nOutputRows )
     , mRasterEntries( rasterEntries )
 {
+  //default to first layer's crs
+  mOutputCrs = mRasterEntries.at( 0 ).raster->crs();
 }
 
+QgsRasterCalculator::QgsRasterCalculator( const QString& formulaString, const QString& outputFile, const QString& outputFormat,
+    const QgsRectangle& outputExtent, const QgsCoordinateReferenceSystem& outputCrs, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry>& rasterEntries )
+    : mFormulaString( formulaString )
+    , mOutputFile( outputFile )
+    , mOutputFormat( outputFormat )
+    , mOutputRectangle( outputExtent )
+    , mOutputCrs( outputCrs )
+    , mNumOutputColumns( nOutputColumns )
+    , mNumOutputRows( nOutputRows )
+    , mRasterEntries( rasterEntries )
+{
+}
+
+
 QgsRasterCalculator::~QgsRasterCalculator()
 {
 }
@@ -62,56 +78,31 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
     return 4;
   }
 
-  double targetGeoTransform[6];
-  outputGeoTransform( targetGeoTransform );
-
-  //open all input rasters for reading
-  QMap< QString, GDALRasterBandH > mInputRasterBands; //raster references and corresponding scanline data
-  QMap< QString, QgsRasterMatrix* > inputScanLineData; //stores raster references and corresponding scanline data
-  QVector< GDALDatasetH > mInputDatasets; //raster references and corresponding dataset
-
+  QMap< QString, QgsRasterBlock* > inputBlocks;
   QVector<QgsRasterCalculatorEntry>::const_iterator it = mRasterEntries.constBegin();
   for ( ; it != mRasterEntries.constEnd(); ++it )
   {
     if ( !it->raster ) // no raster layer in entry
     {
       return 2;
     }
-    GDALDatasetH inputDataset = GDALOpen( TO8F( it->raster->source() ), GA_ReadOnly );
-    if ( !inputDataset )
-    {
-      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 ) )
+    QgsRasterBlock* block = 0;
+    // if crs transform needed
+    if ( it->raster->crs() != mOutputCrs )
     {
-      GDALDatasetH vDataset = GDALAutoCreateWarpedVRT( inputDataset, NULL, NULL, GRA_NearestNeighbour, 0.2, NULL );
-      mInputDatasets.push_back( vDataset );
-      mInputDatasets.push_back( inputDataset );
-      inputDataset = vDataset;
+      QgsRasterProjector* proj = new QgsRasterProjector();
+      proj->setCRS( it->raster->crs(), mOutputCrs );
+      proj->setInput( it->raster->dataProvider()->clone() );
+      proj->setPrecision( QgsRasterProjector::Exact );
+
+      block = proj->block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows );
     }
     else
     {
-      mInputDatasets.push_back( inputDataset );
+      block = it->raster->dataProvider()->block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows );
     }
-
-    GDALRasterBandH inputRasterBand = GDALGetRasterBand( inputDataset, it->bandNumber );
-    if ( inputRasterBand == NULL )
-    {
-      return 2;
-    }
-
-    int nodataSuccess;
-    double nodataValue = GDALGetRasterNoDataValue( inputRasterBand, &nodataSuccess );
-
-    mInputRasterBands.insert( it->ref, inputRasterBand );
-    inputScanLineData.insert( it->ref, new QgsRasterMatrix( mNumOutputColumns, 1, new double[mNumOutputColumns], nodataValue ) );
+    inputBlocks.insert( it->ref, block );
   }
 
   //open output dataset for writing
@@ -120,44 +111,21 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
   {
     return 1;
   }
-  GDALDatasetH outputDataset = openOutputFile( outputDriver );
-
-  //copy the projection info from the first input raster
-  if ( mRasterEntries.size() > 0 )
-  {
-    QgsRasterLayer* rl = mRasterEntries.at( 0 ).raster;
-    if ( rl )
-    {
-      char* crsWKT = 0;
-      OGRSpatialReferenceH ogrSRS = OSRNewSpatialReference( NULL );
-      if ( OSRSetFromUserInput( ogrSRS, rl->crs().authid().toUtf8().constData() ) == OGRERR_NONE )
-      {
-        OSRExportToWkt( ogrSRS, &crsWKT );
-        GDALSetProjection( outputDataset, crsWKT );
-      }
-      else
-      {
-        GDALSetProjection( outputDataset, TO8( rl->crs().toWkt() ) );
-      }
-      OSRDestroySpatialReference( ogrSRS );
-      CPLFree( crsWKT );
-    }
-  }
-
 
+  GDALDatasetH outputDataset = openOutputFile( outputDriver );
+  GDALSetProjection( outputDataset, mOutputCrs.toWkt().toLocal8Bit().data() );
   GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset, 1 );
 
   float outputNodataValue = -FLT_MAX;
   GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
 
-  float* resultScanLine = ( float * ) CPLMalloc( sizeof( float ) * mNumOutputColumns );
-
   if ( p )
   {
     p->setMaximum( mNumOutputRows );
   }
 
   QgsRasterMatrix resultMatrix;
+  resultMatrix.setNodataValue( -FLT_MAX );
 
   //read / write line by line
   for ( int i = 0; i < mNumOutputRows; ++i )
@@ -172,28 +140,7 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
       break;
     }
 
-    //fill buffers
-    QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin();
-    for ( ; bufferIt != inputScanLineData.end(); ++bufferIt )
-    {
-      double sourceTransformation[6];
-      GDALRasterBandH sourceRasterBand = mInputRasterBands[bufferIt.key()];
-      if ( GDALGetGeoTransform( GDALGetBandDataset( sourceRasterBand ), sourceTransformation ) != CE_None )
-      {
-        qWarning( "GDALGetGeoTransform failed!" );
-      }
-
-      float* inputScanLine = ( float * ) CPLMalloc( sizeof( float ) * mNumOutputColumns );
-
-      //the function readRasterPart calls GDALRasterIO (and ev. does some conversion if raster transformations are not the same)
-      readRasterPart( targetGeoTransform, 0, i, mNumOutputColumns, 1, sourceTransformation, sourceRasterBand, inputScanLine );
-      for ( int col = 0; col < mNumOutputColumns; ++col )
-      {
-        bufferIt.value()->data()[col] = ( double )inputScanLine[col];
-      }
-    }
-
-    if ( calcNode->calculate( inputScanLineData, resultMatrix ) )
+    if ( calcNode->calculate( inputBlocks, resultMatrix, i ) )
     {
       bool resultIsNumber = resultMatrix.isNumber();
       float* calcData = new float[mNumOutputColumns];
@@ -225,18 +172,8 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
 
   //close datasets and release memory
   delete calcNode;
-  QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin();
-  for ( ; bufferIt != inputScanLineData.end(); ++bufferIt )
-  {
-    delete bufferIt.value();
-  }
-  inputScanLineData.clear();
-
-  QVector< GDALDatasetH >::iterator datasetIt = mInputDatasets.begin();
-  for ( ; datasetIt != mInputDatasets.end(); ++ datasetIt )
-  {
-    GDALClose( *datasetIt );
-  }
+  qDeleteAll( inputBlocks );
+  inputBlocks.clear();
 
   if ( p && p->wasCanceled() )
   {
@@ -245,7 +182,7 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
     return 3;
   }
   GDALClose( outputDataset );
-  CPLFree( resultScanLine );
+
   return 0;
 }