diff --git a/python/analysis/raster/qgsrastercalcnode.sip b/python/analysis/raster/qgsrastercalcnode.sip
index 0db4effdcef..432f56ba289 100644
--- a/python/analysis/raster/qgsrastercalcnode.sip
+++ b/python/analysis/raster/qgsrastercalcnode.sip
@@ -9,7 +9,8 @@ class QgsRasterCalcNode
     {
       tOperator,
       tNumber,
-      tRasterRef
+      tRasterRef,
+      tMatrix
     };
 
     //! possible operators
@@ -41,6 +42,7 @@ class QgsRasterCalcNode
 
     QgsRasterCalcNode();
     QgsRasterCalcNode( double number );
+    QgsRasterCalcNode( QgsRasterMatrix* matrix );
     QgsRasterCalcNode( Operator op, QgsRasterCalcNode* left, QgsRasterCalcNode* right );
     QgsRasterCalcNode( const QString& rasterName );
     ~QgsRasterCalcNode();
diff --git a/src/analysis/raster/qgsrastercalcnode.cpp b/src/analysis/raster/qgsrastercalcnode.cpp
index a31d3c99287..ecb9e1fef1c 100644
--- a/src/analysis/raster/qgsrastercalcnode.cpp
+++ b/src/analysis/raster/qgsrastercalcnode.cpp
@@ -20,6 +20,7 @@ QgsRasterCalcNode::QgsRasterCalcNode()
     , mLeft( 0 )
     , mRight( 0 )
     , mNumber( 0 )
+    , mMatrix( 0 )
     , mOperator( opNONE )
 {
 }
@@ -29,15 +30,28 @@ QgsRasterCalcNode::QgsRasterCalcNode( double number )
     , mLeft( 0 )
     , mRight( 0 )
     , mNumber( number )
+    , mMatrix( 0 )
     , mOperator( opNONE )
 {
 }
 
+QgsRasterCalcNode::QgsRasterCalcNode( QgsRasterMatrix* matrix )
+    : mType( tMatrix )
+    , mLeft( 0 )
+    , mRight( 0 )
+    , mNumber( 0 )
+    , mMatrix( matrix )
+    , mOperator( opNONE )
+{
+
+}
+
 QgsRasterCalcNode::QgsRasterCalcNode( Operator op, QgsRasterCalcNode* left, QgsRasterCalcNode* right )
     : mType( tOperator )
     , mLeft( left )
     , mRight( right )
     , mNumber( 0 )
+    , mMatrix( 0 )
     , mOperator( op )
 {
 }
@@ -48,6 +62,7 @@ QgsRasterCalcNode::QgsRasterCalcNode( const QString& rasterName )
     , mRight( 0 )
     , mNumber( 0 )
     , mRasterName( rasterName )
+    , mMatrix( 0 )
     , mOperator( opNONE )
 {
   if ( mRasterName.startsWith( '"' ) && mRasterName.endsWith( '"' ) )
@@ -81,8 +96,12 @@ bool QgsRasterCalcNode::calculate( QMap<QString, QgsRasterMatrix*>& rasterData,
 
     int nEntries = ( *it )->nColumns() * ( *it )->nRows();
     double* data = new double[nEntries];
-    memcpy( data, ( *it )->data(), nEntries * sizeof( double ) );
-    result.setData(( *it )->nColumns(), ( *it )->nRows(), data, ( *it )->nodataValue() );
+
+    for ( int i = 0; i < nEntries; ++i )
+    {
+      data[i] = ( *it )->data()[i] == ( *it )->nodataValue() ? result.nodataValue() : ( *it )->data()[i];
+    }
+    result.setData(( *it )->nColumns(), ( *it )->nRows(), data, result.nodataValue() );
     return true;
   }
   else if ( mType == tOperator )
@@ -186,6 +205,17 @@ bool QgsRasterCalcNode::calculate( QMap<QString, QgsRasterMatrix*>& rasterData,
     result.setData( 1, 1, data, result.nodataValue() );
     return true;
   }
+  else if ( mType == tMatrix )
+  {
+    int nEntries = mMatrix->nColumns() * mMatrix->nRows();
+    double* data = new double[nEntries];
+    for ( int i = 0; i < nEntries; ++i )
+    {
+      data[i] = mMatrix->data()[i] == mMatrix->nodataValue() ? result.nodataValue() : mMatrix->data()[i];
+    }
+    result.setData( mMatrix->nColumns(), mMatrix->nRows(), data, result.nodataValue() );
+    return true;
+  }
   return false;
 }
 
diff --git a/src/analysis/raster/qgsrastercalcnode.h b/src/analysis/raster/qgsrastercalcnode.h
index e46dc58422c..ad98f6bb8b8 100644
--- a/src/analysis/raster/qgsrastercalcnode.h
+++ b/src/analysis/raster/qgsrastercalcnode.h
@@ -31,7 +31,8 @@ class ANALYSIS_EXPORT QgsRasterCalcNode
     {
       tOperator = 1,
       tNumber,
-      tRasterRef
+      tRasterRef,
+      tMatrix
     };
 
     //! possible operators
@@ -65,6 +66,7 @@ class ANALYSIS_EXPORT QgsRasterCalcNode
 
     QgsRasterCalcNode();
     QgsRasterCalcNode( double number );
+    QgsRasterCalcNode( QgsRasterMatrix* matrix );
     QgsRasterCalcNode( Operator op, QgsRasterCalcNode* left, QgsRasterCalcNode* right );
     QgsRasterCalcNode( const QString& rasterName );
     ~QgsRasterCalcNode();
@@ -86,6 +88,7 @@ class ANALYSIS_EXPORT QgsRasterCalcNode
     QgsRasterCalcNode* mRight;
     double mNumber;
     QString mRasterName;
+    QgsRasterMatrix* mMatrix;
     Operator mOperator;
 };
 
diff --git a/tests/src/analysis/testqgsrastercalculator.cpp b/tests/src/analysis/testqgsrastercalculator.cpp
index 7c49afff407..028e4e9d00f 100644
--- a/tests/src/analysis/testqgsrastercalculator.cpp
+++ b/tests/src/analysis/testqgsrastercalculator.cpp
@@ -40,6 +40,14 @@ class TestQgsRasterCalculator : public QObject
     void singleOp_data();
     void singleOp(); //test operators which operate on a single value
 
+    void singleOpMatrices(); // test single op using matrix
+    void dualOpNumberMatrix(); // test dual op run on number and matrix
+    void dualOpMatrixNumber(); // test dual op run on matrix and number
+    void dualOpMatrixMatrix(); // test dual op run on matrix and matrix
+
+    void rasterRefOp();
+    void dualOpRasterRaster(); //test dual op on raster ref and raster ref
+
 };
 
 void  TestQgsRasterCalculator::initTestCase()
@@ -178,5 +186,209 @@ void TestQgsRasterCalculator::singleOp()
 
 }
 
+void TestQgsRasterCalculator::singleOpMatrices()
+{
+  // test single op run on matrix
+  double* d = new double[6];
+  d[0] = 1.0;
+  d[1] = 2.0;
+  d[2] = 3.0;
+  d[3] = 4.0;
+  d[4] = 5.0;
+  d[5] = -1.0;
+
+  QgsRasterMatrix m( 2, 3, d, -1.0 );
+
+  QgsRasterCalcNode node( QgsRasterCalcNode::opSIGN, new QgsRasterCalcNode( &m ) , 0 );
+
+  QgsRasterMatrix result;
+  result.setNodataValue( -9999 );
+  QMap<QString, QgsRasterMatrix*> rasterData;
+
+  QVERIFY( node.calculate( rasterData, result ) );
+
+  QCOMPARE( result.data()[0], -d[0] );
+  QCOMPARE( result.data()[1], -d[1] );
+  QCOMPARE( result.data()[2], -d[2] );
+  QCOMPARE( result.data()[3], -d[3] );
+  QCOMPARE( result.data()[4], -d[4] );
+  QCOMPARE( result.data()[5], -9999.0 );
+}
+
+void TestQgsRasterCalculator::dualOpNumberMatrix()
+{
+  // test dual op run on number and matrix
+  double* d = new double[6];
+  d[0] = 1.0;
+  d[1] = 2.0;
+  d[2] = 3.0;
+  d[3] = 4.0;
+  d[4] = 5.0;
+  d[5] = -1.0;
+
+  QgsRasterMatrix m( 2, 3, d, -1.0 );
+
+  QgsRasterCalcNode node( QgsRasterCalcNode::opPLUS, new QgsRasterCalcNode( 5.0 ), new QgsRasterCalcNode( &m ) );
+
+  QgsRasterMatrix result;
+  result.setNodataValue( -9999 );
+  QMap<QString, QgsRasterMatrix*> rasterData;
+
+  QVERIFY( node.calculate( rasterData, result ) );
+
+  QCOMPARE( result.data()[0], 6.0 );
+  QCOMPARE( result.data()[1], 7.0 );
+  QCOMPARE( result.data()[2], 8.0 );
+  QCOMPARE( result.data()[3], 9.0 );
+  QCOMPARE( result.data()[4], 10.0 );
+  QCOMPARE( result.data()[5], -9999.0 );
+
+  //also check adding no data number
+  QgsRasterCalcNode nodeNoData( QgsRasterCalcNode::opPLUS, new QgsRasterCalcNode( -9999 ), new QgsRasterCalcNode( &m ) );
+  QVERIFY( nodeNoData.calculate( rasterData, result ) );
+  QCOMPARE( result.data()[0], -9999.0 );
+}
+
+void TestQgsRasterCalculator::dualOpMatrixNumber()
+{
+  // test dual op run on matrix and number
+  double* d = new double[6];
+  d[0] = 1.0;
+  d[1] = 2.0;
+  d[2] = 3.0;
+  d[3] = 4.0;
+  d[4] = 5.0;
+  d[5] = -1.0;
+
+  QgsRasterMatrix m( 2, 3, d, -1.0 );
+
+  QgsRasterCalcNode node( QgsRasterCalcNode::opPLUS, new QgsRasterCalcNode( &m ), new QgsRasterCalcNode( 5.0 ) );
+
+  QgsRasterMatrix result;
+  result.setNodataValue( -9999 );
+  QMap<QString, QgsRasterMatrix*> rasterData;
+
+  QVERIFY( node.calculate( rasterData, result ) );
+
+  QCOMPARE( result.data()[0], 6.0 );
+  QCOMPARE( result.data()[1], 7.0 );
+  QCOMPARE( result.data()[2], 8.0 );
+  QCOMPARE( result.data()[3], 9.0 );
+  QCOMPARE( result.data()[4], 10.0 );
+  QCOMPARE( result.data()[5], -9999.0 );
+
+  //also check adding no data number
+  QgsRasterCalcNode nodeNoData( QgsRasterCalcNode::opPLUS, new QgsRasterCalcNode( &m ), new QgsRasterCalcNode( -9999 ) );
+  QVERIFY( nodeNoData.calculate( rasterData, result ) );
+  QCOMPARE( result.data()[0], -9999.0 );
+}
+
+void TestQgsRasterCalculator::dualOpMatrixMatrix()
+{
+  // test dual op run on matrix and matrix
+  double* d = new double[6];
+  d[0] = 1.0;
+  d[1] = 2.0;
+  d[2] = -2.0;
+  d[3] = -1.0; //nodata
+  d[4] = 5.0;
+  d[5] = -1.0; //nodata
+  QgsRasterMatrix m1( 2, 3, d, -1.0 );
+
+  double* d2 = new double[6];
+  d2[0] = -1.0;
+  d2[1] = -2.0; //nodata
+  d2[2] = 13.0;
+  d2[3] = -2.0; //nodata
+  d2[4] = 15.0;
+  d2[5] = -1.0;
+  QgsRasterMatrix m2( 2, 3, d2, -2.0 ); //different no data value
+
+  QgsRasterCalcNode node( QgsRasterCalcNode::opPLUS, new QgsRasterCalcNode( &m1 ), new QgsRasterCalcNode( &m2 ) );
+
+  QgsRasterMatrix result;
+  result.setNodataValue( -9999 );
+  QMap<QString, QgsRasterMatrix*> rasterData;
+
+  QVERIFY( node.calculate( rasterData, result ) );
+
+  QCOMPARE( result.data()[0], 0.0 );
+  QCOMPARE( result.data()[1], -9999.0 );
+  QCOMPARE( result.data()[2], 11.0 );
+  QCOMPARE( result.data()[3], -9999.0 );
+  QCOMPARE( result.data()[4], 20.0 );
+  QCOMPARE( result.data()[5], -9999.0 );
+}
+
+void TestQgsRasterCalculator::rasterRefOp()
+{
+  // test single op run on raster ref
+  QgsRasterCalcNode node( QgsRasterCalcNode::opSIGN, new QgsRasterCalcNode( "raster" ), 0 );
+
+  QgsRasterMatrix result;
+  result.setNodataValue( -9999 );
+  QMap<QString, QgsRasterMatrix*> rasterData;
+
+  //first test invalid raster ref
+  QVERIFY( !node.calculate( rasterData, result ) );
+
+  //now create raster ref
+  double* d = new double[6];
+  d[0] = 1.0;
+  d[1] = 2.0;
+  d[2] = 3.0;
+  d[3] = 4.0;
+  d[4] = 5.0;
+  d[5] = -1.0;
+  QgsRasterMatrix m( 2, 3, d, -1.0 );
+  rasterData.insert( "raster", &m );
+
+  QVERIFY( node.calculate( rasterData, result ) );
+  QCOMPARE( result.data()[0], -d[0] );
+  QCOMPARE( result.data()[1], -d[1] );
+  QCOMPARE( result.data()[2], -d[2] );
+  QCOMPARE( result.data()[3], -d[3] );
+  QCOMPARE( result.data()[4], -d[4] );
+  QCOMPARE( result.data()[5], -9999.0 );
+}
+
+void TestQgsRasterCalculator::dualOpRasterRaster()
+{
+  // test dual op run on matrix and matrix
+  double* d = new double[6];
+  d[0] = 1.0;
+  d[1] = 2.0;
+  d[2] = -2.0;
+  d[3] = -1.0; //nodata
+  d[4] = 5.0;
+  d[5] = -1.0; //nodata
+  QgsRasterMatrix m1( 2, 3, d, -1.0 );
+  QMap<QString, QgsRasterMatrix*> rasterData;
+  rasterData.insert( "raster1", &m1 );
+
+  double* d2 = new double[6];
+  d2[0] = -1.0;
+  d2[1] = -2.0; //nodata
+  d2[2] = 13.0;
+  d2[3] = -2.0; //nodata
+  d2[4] = 15.0;
+  d2[5] = -1.0;
+  QgsRasterMatrix m2( 2, 3, d2, -2.0 ); //different no data value
+  rasterData.insert( "raster2", &m2 );
+
+  QgsRasterCalcNode node( QgsRasterCalcNode::opPLUS, new QgsRasterCalcNode( "raster1" ), new QgsRasterCalcNode( "raster2" ) );
+
+  QgsRasterMatrix result;
+  result.setNodataValue( -9999 );
+
+  QVERIFY( node.calculate( rasterData, result ) );
+  QCOMPARE( result.data()[0], 0.0 );
+  QCOMPARE( result.data()[1], -9999.0 );
+  QCOMPARE( result.data()[2], 11.0 );
+  QCOMPARE( result.data()[3], -9999.0 );
+  QCOMPARE( result.data()[4], 20.0 );
+  QCOMPARE( result.data()[5], -9999.0 );
+}
+
 QTEST_MAIN( TestQgsRasterCalculator )
 #include "testqgsrastercalculator.moc"