[rastercalc] Switch all internal calculations to doubles

...for more accurate calculations (fix #9081)

Also:
- Add test suite for raster calculator
- Fix errors with log/log10 calculation and inputs <= 0
- Fix loss of nodata values in certain circumstances

Author: nyalldawson

python/analysis/raster/qgsrastermatrix.sip

@@ -48,12 +48,12 @@ class QgsRasterMatrix
 
     /**Returns data array (but not ownership)*/
     //! @note not available in python bindings
-    // float* data();
+    // double* data();
     /**Returns data and ownership. Sets data and nrows, ncols of this matrix to 0*/
     //! @note not available in python bindings
-    // float* takeData();
+    // double* takeData();
 
-    void setData( int cols, int rows, float* data, double nodataValue );
+    void setData( int cols, int rows, double* data, double nodataValue );
 
     int nColumns() const;
     int nRows() const;

src/analysis/raster/qgsrastercalcnode.cpp

@@ -80,15 +80,17 @@ bool QgsRasterCalcNode::calculate( QMap<QString, QgsRasterMatrix*>& rasterData,
     }
 
     int nEntries = ( *it )->nColumns() * ( *it )->nRows();
-    float* data = new float[nEntries];
-    memcpy( data, ( *it )->data(), nEntries * sizeof( float ) );
+    double* data = new double[nEntries];
+    memcpy( data, ( *it )->data(), nEntries * sizeof( double ) );
     result.setData(( *it )->nColumns(), ( *it )->nRows(), data, ( *it )->nodataValue() );
     return true;
   }
   else if ( mType == tOperator )
   {
     QgsRasterMatrix leftMatrix, rightMatrix;
-    QgsRasterMatrix resultMatrix;
+    leftMatrix.setNodataValue( result.nodataValue() );
+    rightMatrix.setNodataValue( result.nodataValue() );
+
     if ( !mLeft || !mLeft->calculate( rasterData, leftMatrix ) )
     {
       return false;
@@ -179,9 +181,9 @@ bool QgsRasterCalcNode::calculate( QMap<QString, QgsRasterMatrix*>& rasterData,
   }
   else if ( mType == tNumber )
   {
-    float* data = new float[1];
+    double* data = new double[1];
     data[0] = mNumber;
-    result.setData( 1, 1, data, -FLT_MAX );
+    result.setData( 1, 1, data, result.nodataValue() );
     return true;
   }
   return false;

src/analysis/raster/qgsrastercalculator.cpp

@@ -111,7 +111,7 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
     double nodataValue = GDALGetRasterNoDataValue( inputRasterBand, &nodataSuccess );
 
     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 double[mNumOutputColumns], nodataValue ) );
   }
 
   //open output dataset for writing
@@ -183,35 +183,25 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
         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, bufferIt.value()->data() );
+      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 ) )
     {
       bool resultIsNumber = resultMatrix.isNumber();
-      float* calcData;
-
-      if ( resultIsNumber ) //scalar result. Insert number for every pixel
-      {
-        calcData = new float[mNumOutputColumns];
-        for ( int j = 0; j < mNumOutputColumns; ++j )
-        {
-          calcData[j] = resultMatrix.number();
-        }
-      }
-      else //result is real matrix
-      {
-        calcData = resultMatrix.data();
-      }
+      float* calcData = new float[mNumOutputColumns];
 
-      //replace all matrix nodata values with output nodatas
       for ( int j = 0; j < mNumOutputColumns; ++j )
       {
-        if ( calcData[j] == resultMatrix.nodataValue() )
-        {
-          calcData[j] = outputNodataValue;
-        }
+        double result = resultIsNumber ? resultMatrix.number() : resultMatrix.data()[j];
+        calcData[j] = ( calcData[j] == resultMatrix.nodataValue() ? outputNodataValue : ( float ) result );
       }
 
       //write scanline to the dataset