Skip to content
Permalink
Browse files

[FEATURE]: Add rastercalculator

git-svn-id: http://svn.osgeo.org/qgis/trunk/qgis@14442 c8812cc2-4d05-0410-92ff-de0c093fc19c
  • Loading branch information
mhugent
mhugent committed Oct 28, 2010
1 parent 01705a4 commit 8a645170b08ef997ee1b729f31f435c013aacca0
@@ -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)
@@ -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)
@@ -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}
@@ -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);
}
@@ -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;
}


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

0 comments on commit 8a64517

Please sign in to comment.
You can’t perform that action at this time.