[FEATURE]: Add rastercalculator

git-svn-id: http://svn.osgeo.org/qgis/trunk@14442 c8812cc2-4d05-0410-92ff-de0c093fc19c
qgis · Oct 28, 2010 · 9e62825 · 9e62825
1 parent 9f9a9cf
commit 9e62825
Show file tree

Hide file tree

Showing 16 changed files with 2,640 additions and 428 deletions.
diff --git a/cmake/Bison.cmake b/cmake/Bison.cmake
@@ -65,3 +65,32 @@ MACRO(ADD_BISON_FILES _sources )
     SET(${_sources} ${${_sources}} ${_out} )
   ENDFOREACH (_current_FILE)
 ENDMACRO(ADD_BISON_FILES)
+
+MACRO(ADD_BISON_FILES_PREFIX _sources prefix)
+  FIND_BISON()
+
+  FOREACH (_current_FILE ${ARGN})
+    GET_FILENAME_COMPONENT(_in ${_current_FILE} ABSOLUTE)
+    GET_FILENAME_COMPONENT(_basename ${_current_FILE} NAME_WE)
+
+    SET(_out ${CMAKE_CURRENT_BINARY_DIR}/${_basename}.cpp)
+
+
+    # bison options:
+    # -t add debugging facilities
+    # -d produce additional header file (used in parser.l)
+    # -v produce additional *.output file with parser states
+
+    ADD_CUSTOM_COMMAND(
+      OUTPUT ${_out}
+      COMMAND ${BISON_EXECUTABLE}
+      ARGS
+      -p ${prefix}
+      -o${_out} -d -v -t
+      ${_in}
+      DEPENDS ${_in}
+      )
+
+    SET(${_sources} ${${_sources}} ${_out} )
+  ENDFOREACH (_current_FILE)
+ENDMACRO(ADD_BISON_FILES)
diff --git a/cmake/Flex.cmake b/cmake/Flex.cmake
@@ -47,3 +47,30 @@ MACRO(ADD_FLEX_FILES _sources )
       SET(${_sources} ${${_sources}} ${_out} )
    ENDFOREACH (_current_FILE)
 ENDMACRO(ADD_FLEX_FILES)
+
+
+MACRO(ADD_FLEX_FILES_PREFIX _sources prefix )
+    FIND_FLEX()
+
+    FOREACH (_current_FILE ${ARGN})
+      GET_FILENAME_COMPONENT(_in ${_current_FILE} ABSOLUTE)
+      GET_FILENAME_COMPONENT(_basename ${_current_FILE} NAME_WE)
+
+      SET(_out ${CMAKE_CURRENT_BINARY_DIR}/flex_${_basename}.cpp)
+
+
+      # -d option for flex means that it will produce output to stderr while analyzing 
+
+      ADD_CUSTOM_COMMAND(
+         OUTPUT ${_out}
+         COMMAND ${FLEX_EXECUTABLE}
+         ARGS
+	 -P ${prefix}
+         -o${_out} -d
+         ${_in}
+         DEPENDS ${_in}
+      )
+
+      SET(${_sources} ${${_sources}} ${_out} )
+   ENDFOREACH (_current_FILE)
+ENDMACRO(ADD_FLEX_FILES_PREFIX)
diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt
@@ -29,11 +29,19 @@ SET(QGIS_ANALYSIS_SRCS
   raster/qgsslopefilter.cpp
   raster/qgsaspectfilter.cpp
   raster/qgstotalcurvaturefilter.cpp
+  raster/qgsrastercalcnode.cpp
+  raster/qgsrastercalculator.cpp
+  raster/qgsrastermatrix.cpp
   vector/qgsgeometryanalyzer.cpp
   vector/qgszonalstatistics.cpp
   vector/qgsoverlayanalyzer.cpp
 )
 
+INCLUDE_DIRECTORIES(BEFORE raster)
+ADD_FLEX_FILES_PREFIX(QGIS_ANALYSIS_SRCS raster raster/qgsrastercalclexer.ll)
+
+ADD_BISON_FILES_PREFIX(QGIS_ANALYSIS_SRCS raster raster/qgsrastercalcparser.yy)
+
 SET(QGIS_ANALYSIS_MOC_HDRS
 )
 
@@ -45,6 +53,7 @@ INCLUDE_DIRECTORIES(
   ${CMAKE_CURRENT_SOURCE_DIR}/../core/
   ${CMAKE_CURRENT_SOURCE_DIR}/../core/renderer
   ${CMAKE_CURRENT_SOURCE_DIR}/../core/spatialindex
+  ${CMAKE_CURRENT_SOURCE_DIR}/../core/raster
   interpolation
   ${PROJ_INCLUDE_DIR}
   ${GEOS_INCLUDE_DIR}

diff --git a/src/analysis/raster/qgsrastercalclexer.ll b/src/analysis/raster/qgsrastercalclexer.ll
@@ -0,0 +1,75 @@
+/***************************************************************************
+                          qgsrastercalclexer.ll
+          Rules for lexical analysis of raster calc strings done by Flex
+                          --------------------
+    begin                : 2010-10-01
+    copyright            : (C) 2010 by Marco Hugentobler
+    email                : marco dot hugentobler at sourcepole dot ch
+***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+%option noyywrap
+%option case-insensitive
+%option never-interactive
+
+ // ensure that lexer will be 8-bit (and not just 7-bit)
+%option 8bit
+
+%{
+  //directly included in the output program
+  #include "qgsrastercalcnode.h"
+  #include "qgsrastercalcparser.hpp"
+
+  // if not defined, searches for isatty()
+  // which doesn't in MSVC compiler
+  #define YY_NEVER_INTERACTIVE 1
+
+  #ifdef _MSC_VER
+  #define YY_NO_UNISTD_H
+  #endif
+%}
+
+white       [ \t\r\n]+
+
+dig     [0-9]
+num1    {dig}+\.?([eE][-+]?{dig}+)?
+num2    {dig}*\.{dig}+([eE][-+]?{dig}+)?
+number  {num1}|{num2}
+
+non_ascii    [\x80-\xFF]
+raster_ref_char  [A-Za-z0-9_]|{non_ascii}
+raster_band_ref ({raster_ref_char}+)@{dig}
+
+%%
+
+"sqrt" { rasterlval.op = QgsRasterCalcNode::opSQRT; return FUNCTION;}
+"sin"  { rasterlval.op = QgsRasterCalcNode::opSIN; return FUNCTION;}
+"cos"  { rasterlval.op = QgsRasterCalcNode::opCOS; return FUNCTION;}
+"tan"  { rasterlval.op = QgsRasterCalcNode::opTAN; return FUNCTION;}
+"asin" { rasterlval.op = QgsRasterCalcNode::opASIN; return FUNCTION;}
+"acos" { rasterlval.op = QgsRasterCalcNode::opACOS; return FUNCTION;}
+"atan" { rasterlval.op = QgsRasterCalcNode::opATAN; return FUNCTION;}
+
+[+-/*^] { return yytext[0]; }
+
+[()] { return yytext[0]; }
+
+{number} { rasterlval.number  = atof(rastertext); return NUMBER; }
+
+{raster_band_ref} { return RASTER_BAND_REF; }
+
+{white}    /* skip blanks and tabs */
+%%
+
+void set_raster_input_buffer(const char* buffer)
+{
+  raster_scan_string(buffer);
+}
diff --git a/src/analysis/raster/qgsrastercalcnode.cpp b/src/analysis/raster/qgsrastercalcnode.cpp
@@ -0,0 +1,119 @@
+#include "qgsrastercalcnode.h"
+
+QgsRasterCalcNode::QgsRasterCalcNode(): mLeft( 0 ), mRight( 0 ), mRasterMatrix( 0 ), mNumber( 0 )
+{
+}
+
+QgsRasterCalcNode::QgsRasterCalcNode( double number ): mType( tNumber ), mLeft( 0 ), mRight( 0 ), mRasterMatrix( 0 ), mNumber( number )
+{
+}
+
+QgsRasterCalcNode::QgsRasterCalcNode( Operator op, QgsRasterCalcNode* left, QgsRasterCalcNode* right ): mType( tOperator ), mLeft( left ), mRight( right ), mRasterMatrix( 0 ), mNumber( 0 ), mOperator( op )
+{
+}
+
+QgsRasterCalcNode::QgsRasterCalcNode( const QString& rasterName ): mType( tRasterRef ), mLeft( 0 ), mRight( 0 ), mRasterMatrix( 0 ), mNumber( 0 ), mRasterName( rasterName )
+{
+}
+
+QgsRasterCalcNode::~QgsRasterCalcNode()
+{
+  if( mLeft )
+  {
+    delete mLeft;
+  }
+  if( mRight )
+  {
+    delete mRight;
+  }
+}
+
+bool QgsRasterCalcNode::calculate( QMap<QString, QgsRasterMatrix*>& rasterData, QgsRasterMatrix& result ) 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 );
+    if( it == rasterData.end() )
+    {
+      return false;
+    }
+
+    int nEntries = ( *it )->nColumns() * ( *it )->nRows();
+    float* data = new float[nEntries];
+    memcpy( data, ( *it )->data(), nEntries * sizeof( float ) );
+    result.setData(( *it )->nColumns(), ( *it )->nRows(), data );
+    return true;
+  }
+  else if( mType == tOperator )
+  {
+    QgsRasterMatrix leftMatrix, rightMatrix;
+    QgsRasterMatrix resultMatrix;
+    if( !mLeft || !mLeft->calculate( rasterData, leftMatrix ) )
+    {
+      return false;
+    }
+    if( mRight && !mRight->calculate( rasterData, rightMatrix ) )
+    {
+      return false;
+    }
+
+    switch( mOperator )
+    {
+      case opPLUS:
+        leftMatrix.add( rightMatrix );
+        break;
+      case opMINUS:
+        leftMatrix.subtract( rightMatrix );
+        break;
+      case opMUL:
+        leftMatrix.multiply( rightMatrix );
+        break;
+      case opDIV:
+        leftMatrix.divide( rightMatrix );
+        break;
+      case opPOW:
+        leftMatrix.power( rightMatrix );
+        break;
+      case opSQRT:
+        leftMatrix.squareRoot();
+        break;
+      case opSIN:
+        leftMatrix.sinus();
+        break;
+      case opCOS:
+        leftMatrix.cosinus();
+        break;
+      case opTAN:
+        leftMatrix.tangens();
+        break;
+      case opASIN:
+        leftMatrix.asinus();
+        break;
+      case opACOS:
+        leftMatrix.acosinus();
+        break;
+      case opATAN:
+        leftMatrix.atangens();
+        break;
+      default:
+        return false;
+    }
+    int newNColumns = leftMatrix.nColumns();
+    int newNRows = leftMatrix.nRows();
+    result.setData( newNColumns, newNRows, leftMatrix.takeData() );
+    return true;
+  }
+  else if( mType == tNumber )
+  {
+    float* data = new float[1];
+    data[0] = mNumber;
+    result.setData( 1, 1, data );
+    return true;
+  }
+  return false;
+}
+
+
diff --git a/src/analysis/raster/qgsrastercalcnode.h b/src/analysis/raster/qgsrastercalcnode.h
@@ -0,0 +1,79 @@
+/***************************************************************************
+                          qgsrastercalcnode.h
+            Node for raster calculator tree
+                          --------------------
+    begin                : 2010-10-23
+    copyright            : (C) 20010 by Marco Hugentobler
+    email                : marco dot hugentobler at sourcepole dot ch
+***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef QGSRASTERCALCNODE_H
+#define QGSRASTERCALCNODE_H
+
+#include "qgsrastermatrix.h"
+#include <QMap>
+#include <QString>
+
+class ANALYSIS_EXPORT QgsRasterCalcNode
+{
+  public:
+    //! defines possible types of node
+    enum Type
+    {
+      tOperator = 1,
+      tNumber,
+      tRasterRef
+    };
+
+    //! possible operators
+    enum Operator
+    {
+      opPLUS,
+      opMINUS,
+      opMUL,
+      opDIV,
+      opPOW,
+      opSQRT,
+      opSIN,
+      opCOS,
+      opTAN,
+      opASIN,
+      opACOS,
+      opATAN
+    };
+
+    QgsRasterCalcNode();
+    QgsRasterCalcNode( double number );
+    QgsRasterCalcNode( Operator op, QgsRasterCalcNode* left, QgsRasterCalcNode* right );
+    QgsRasterCalcNode( const QString& rasterName );
+    ~QgsRasterCalcNode();
+
+    Type type() const { return mType; }
+
+    //set left node
+    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;
+
+  private:
+    Type mType;
+    QgsRasterCalcNode* mLeft;
+    QgsRasterCalcNode* mRight;
+    QgsRasterMatrix* mRasterMatrix;
+    QString mRasterName;
+    double mNumber;
+    Operator mOperator;
+};
+
+#endif // QGSRASTERCALCNODE_H