Skip to content

Commit 9e62825

Browse files
author
mhugent
committedOct 28, 2010
[FEATURE]: Add rastercalculator
git-svn-id: http://svn.osgeo.org/qgis/trunk@14442 c8812cc2-4d05-0410-92ff-de0c093fc19c
·
ltr-3_40Before-merge-new_vector_api
1 parent 9f9a9cf commit 9e62825

16 files changed

+2640
-428
lines changed
 

‎cmake/Bison.cmake

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -65,3 +65,32 @@ MACRO(ADD_BISON_FILES _sources )
6565
SET(${_sources} ${${_sources}} ${_out} )
6666
ENDFOREACH (_current_FILE)
6767
ENDMACRO(ADD_BISON_FILES)
68+
69+
MACRO(ADD_BISON_FILES_PREFIX _sources prefix)
70+
FIND_BISON()
71+
72+
FOREACH (_current_FILE ${ARGN})
73+
GET_FILENAME_COMPONENT(_in ${_current_FILE} ABSOLUTE)
74+
GET_FILENAME_COMPONENT(_basename ${_current_FILE} NAME_WE)
75+
76+
SET(_out ${CMAKE_CURRENT_BINARY_DIR}/${_basename}.cpp)
77+
78+
79+
# bison options:
80+
# -t add debugging facilities
81+
# -d produce additional header file (used in parser.l)
82+
# -v produce additional *.output file with parser states
83+
84+
ADD_CUSTOM_COMMAND(
85+
OUTPUT ${_out}
86+
COMMAND ${BISON_EXECUTABLE}
87+
ARGS
88+
-p ${prefix}
89+
-o${_out} -d -v -t
90+
${_in}
91+
DEPENDS ${_in}
92+
)
93+
94+
SET(${_sources} ${${_sources}} ${_out} )
95+
ENDFOREACH (_current_FILE)
96+
ENDMACRO(ADD_BISON_FILES)

‎cmake/Flex.cmake

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -47,3 +47,30 @@ MACRO(ADD_FLEX_FILES _sources )
4747
SET(${_sources} ${${_sources}} ${_out} )
4848
ENDFOREACH (_current_FILE)
4949
ENDMACRO(ADD_FLEX_FILES)
50+
51+
52+
MACRO(ADD_FLEX_FILES_PREFIX _sources prefix )
53+
FIND_FLEX()
54+
55+
FOREACH (_current_FILE ${ARGN})
56+
GET_FILENAME_COMPONENT(_in ${_current_FILE} ABSOLUTE)
57+
GET_FILENAME_COMPONENT(_basename ${_current_FILE} NAME_WE)
58+
59+
SET(_out ${CMAKE_CURRENT_BINARY_DIR}/flex_${_basename}.cpp)
60+
61+
62+
# -d option for flex means that it will produce output to stderr while analyzing
63+
64+
ADD_CUSTOM_COMMAND(
65+
OUTPUT ${_out}
66+
COMMAND ${FLEX_EXECUTABLE}
67+
ARGS
68+
-P ${prefix}
69+
-o${_out} -d
70+
${_in}
71+
DEPENDS ${_in}
72+
)
73+
74+
SET(${_sources} ${${_sources}} ${_out} )
75+
ENDFOREACH (_current_FILE)
76+
ENDMACRO(ADD_FLEX_FILES_PREFIX)

‎src/analysis/CMakeLists.txt

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,11 +29,19 @@ SET(QGIS_ANALYSIS_SRCS
2929
raster/qgsslopefilter.cpp
3030
raster/qgsaspectfilter.cpp
3131
raster/qgstotalcurvaturefilter.cpp
32+
raster/qgsrastercalcnode.cpp
33+
raster/qgsrastercalculator.cpp
34+
raster/qgsrastermatrix.cpp
3235
vector/qgsgeometryanalyzer.cpp
3336
vector/qgszonalstatistics.cpp
3437
vector/qgsoverlayanalyzer.cpp
3538
)
3639

40+
INCLUDE_DIRECTORIES(BEFORE raster)
41+
ADD_FLEX_FILES_PREFIX(QGIS_ANALYSIS_SRCS raster raster/qgsrastercalclexer.ll)
42+
43+
ADD_BISON_FILES_PREFIX(QGIS_ANALYSIS_SRCS raster raster/qgsrastercalcparser.yy)
44+
3745
SET(QGIS_ANALYSIS_MOC_HDRS
3846
)
3947

@@ -45,6 +53,7 @@ INCLUDE_DIRECTORIES(
4553
${CMAKE_CURRENT_SOURCE_DIR}/../core/
4654
${CMAKE_CURRENT_SOURCE_DIR}/../core/renderer
4755
${CMAKE_CURRENT_SOURCE_DIR}/../core/spatialindex
56+
${CMAKE_CURRENT_SOURCE_DIR}/../core/raster
4857
interpolation
4958
${PROJ_INCLUDE_DIR}
5059
${GEOS_INCLUDE_DIR}

‎src/analysis/raster/qgsrastercalclexer.ll

Lines changed: 75 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,75 @@
1+
/***************************************************************************
2+
qgsrastercalclexer.ll
3+
Rules for lexical analysis of raster calc strings done by Flex
4+
--------------------
5+
begin : 2010-10-01
6+
copyright : (C) 2010 by Marco Hugentobler
7+
email : marco dot hugentobler at sourcepole dot ch
8+
***************************************************************************/
9+
10+
/***************************************************************************
11+
* *
12+
* This program is free software; you can redistribute it and/or modify *
13+
* it under the terms of the GNU General Public License as published by *
14+
* the Free Software Foundation; either version 2 of the License, or *
15+
* (at your option) any later version. *
16+
* *
17+
***************************************************************************/
18+
19+
%option noyywrap
20+
%option case-insensitive
21+
%option never-interactive
22+
23+
// ensure that lexer will be 8-bit (and not just 7-bit)
24+
%option 8bit
25+
26+
%{
27+
//directly included in the output program
28+
#include "qgsrastercalcnode.h"
29+
#include "qgsrastercalcparser.hpp"
30+
31+
// if not defined, searches for isatty()
32+
// which doesn't in MSVC compiler
33+
#define YY_NEVER_INTERACTIVE 1
34+
35+
#ifdef _MSC_VER
36+
#define YY_NO_UNISTD_H
37+
#endif
38+
%}
39+
40+
white [ \t\r\n]+
41+
42+
dig [0-9]
43+
num1 {dig}+\.?([eE][-+]?{dig}+)?
44+
num2 {dig}*\.{dig}+([eE][-+]?{dig}+)?
45+
number {num1}|{num2}
46+
47+
non_ascii [\x80-\xFF]
48+
raster_ref_char [A-Za-z0-9_]|{non_ascii}
49+
raster_band_ref ({raster_ref_char}+)@{dig}
50+
51+
%%
52+
53+
"sqrt" { rasterlval.op = QgsRasterCalcNode::opSQRT; return FUNCTION;}
54+
"sin" { rasterlval.op = QgsRasterCalcNode::opSIN; return FUNCTION;}
55+
"cos" { rasterlval.op = QgsRasterCalcNode::opCOS; return FUNCTION;}
56+
"tan" { rasterlval.op = QgsRasterCalcNode::opTAN; return FUNCTION;}
57+
"asin" { rasterlval.op = QgsRasterCalcNode::opASIN; return FUNCTION;}
58+
"acos" { rasterlval.op = QgsRasterCalcNode::opACOS; return FUNCTION;}
59+
"atan" { rasterlval.op = QgsRasterCalcNode::opATAN; return FUNCTION;}
60+
61+
[+-/*^] { return yytext[0]; }
62+
63+
[()] { return yytext[0]; }
64+
65+
{number} { rasterlval.number = atof(rastertext); return NUMBER; }
66+
67+
{raster_band_ref} { return RASTER_BAND_REF; }
68+
69+
{white} /* skip blanks and tabs */
70+
%%
71+
72+
void set_raster_input_buffer(const char* buffer)
73+
{
74+
raster_scan_string(buffer);
75+
}

‎src/analysis/raster/qgsrastercalcnode.cpp

Lines changed: 119 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,119 @@
1+
#include "qgsrastercalcnode.h"
2+
3+
QgsRasterCalcNode::QgsRasterCalcNode(): mLeft( 0 ), mRight( 0 ), mRasterMatrix( 0 ), mNumber( 0 )
4+
{
5+
}
6+
7+
QgsRasterCalcNode::QgsRasterCalcNode( double number ): mType( tNumber ), mLeft( 0 ), mRight( 0 ), mRasterMatrix( 0 ), mNumber( number )
8+
{
9+
}
10+
11+
QgsRasterCalcNode::QgsRasterCalcNode( Operator op, QgsRasterCalcNode* left, QgsRasterCalcNode* right ): mType( tOperator ), mLeft( left ), mRight( right ), mRasterMatrix( 0 ), mNumber( 0 ), mOperator( op )
12+
{
13+
}
14+
15+
QgsRasterCalcNode::QgsRasterCalcNode( const QString& rasterName ): mType( tRasterRef ), mLeft( 0 ), mRight( 0 ), mRasterMatrix( 0 ), mNumber( 0 ), mRasterName( rasterName )
16+
{
17+
}
18+
19+
QgsRasterCalcNode::~QgsRasterCalcNode()
20+
{
21+
if( mLeft )
22+
{
23+
delete mLeft;
24+
}
25+
if( mRight )
26+
{
27+
delete mRight;
28+
}
29+
}
30+
31+
bool QgsRasterCalcNode::calculate( QMap<QString, QgsRasterMatrix*>& rasterData, QgsRasterMatrix& result ) const
32+
{
33+
//if type is raster ref: return a copy of the corresponding matrix
34+
35+
//if type is operator, call the proper matrix operations
36+
if( mType == tRasterRef )
37+
{
38+
QMap<QString, QgsRasterMatrix*>::iterator it = rasterData.find( mRasterName );
39+
if( it == rasterData.end() )
40+
{
41+
return false;
42+
}
43+
44+
int nEntries = ( *it )->nColumns() * ( *it )->nRows();
45+
float* data = new float[nEntries];
46+
memcpy( data, ( *it )->data(), nEntries * sizeof( float ) );
47+
result.setData(( *it )->nColumns(), ( *it )->nRows(), data );
48+
return true;
49+
}
50+
else if( mType == tOperator )
51+
{
52+
QgsRasterMatrix leftMatrix, rightMatrix;
53+
QgsRasterMatrix resultMatrix;
54+
if( !mLeft || !mLeft->calculate( rasterData, leftMatrix ) )
55+
{
56+
return false;
57+
}
58+
if( mRight && !mRight->calculate( rasterData, rightMatrix ) )
59+
{
60+
return false;
61+
}
62+
63+
switch( mOperator )
64+
{
65+
case opPLUS:
66+
leftMatrix.add( rightMatrix );
67+
break;
68+
case opMINUS:
69+
leftMatrix.subtract( rightMatrix );
70+
break;
71+
case opMUL:
72+
leftMatrix.multiply( rightMatrix );
73+
break;
74+
case opDIV:
75+
leftMatrix.divide( rightMatrix );
76+
break;
77+
case opPOW:
78+
leftMatrix.power( rightMatrix );
79+
break;
80+
case opSQRT:
81+
leftMatrix.squareRoot();
82+
break;
83+
case opSIN:
84+
leftMatrix.sinus();
85+
break;
86+
case opCOS:
87+
leftMatrix.cosinus();
88+
break;
89+
case opTAN:
90+
leftMatrix.tangens();
91+
break;
92+
case opASIN:
93+
leftMatrix.asinus();
94+
break;
95+
case opACOS:
96+
leftMatrix.acosinus();
97+
break;
98+
case opATAN:
99+
leftMatrix.atangens();
100+
break;
101+
default:
102+
return false;
103+
}
104+
int newNColumns = leftMatrix.nColumns();
105+
int newNRows = leftMatrix.nRows();
106+
result.setData( newNColumns, newNRows, leftMatrix.takeData() );
107+
return true;
108+
}
109+
else if( mType == tNumber )
110+
{
111+
float* data = new float[1];
112+
data[0] = mNumber;
113+
result.setData( 1, 1, data );
114+
return true;
115+
}
116+
return false;
117+
}
118+
119+

‎src/analysis/raster/qgsrastercalcnode.h

Lines changed: 79 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
/***************************************************************************
2+
qgsrastercalcnode.h
3+
Node for raster calculator tree
4+
--------------------
5+
begin : 2010-10-23
6+
copyright : (C) 20010 by Marco Hugentobler
7+
email : marco dot hugentobler at sourcepole dot ch
8+
***************************************************************************/
9+
10+
/***************************************************************************
11+
* *
12+
* This program is free software; you can redistribute it and/or modify *
13+
* it under the terms of the GNU General Public License as published by *
14+
* the Free Software Foundation; either version 2 of the License, or *
15+
* (at your option) any later version. *
16+
* *
17+
***************************************************************************/
18+
19+
#ifndef QGSRASTERCALCNODE_H
20+
#define QGSRASTERCALCNODE_H
21+
22+
#include "qgsrastermatrix.h"
23+
#include <QMap>
24+
#include <QString>
25+
26+
class ANALYSIS_EXPORT QgsRasterCalcNode
27+
{
28+
public:
29+
//! defines possible types of node
30+
enum Type
31+
{
32+
tOperator = 1,
33+
tNumber,
34+
tRasterRef
35+
};
36+
37+
//! possible operators
38+
enum Operator
39+
{
40+
opPLUS,
41+
opMINUS,
42+
opMUL,
43+
opDIV,
44+
opPOW,
45+
opSQRT,
46+
opSIN,
47+
opCOS,
48+
opTAN,
49+
opASIN,
50+
opACOS,
51+
opATAN
52+
};
53+
54+
QgsRasterCalcNode();
55+
QgsRasterCalcNode( double number );
56+
QgsRasterCalcNode( Operator op, QgsRasterCalcNode* left, QgsRasterCalcNode* right );
57+
QgsRasterCalcNode( const QString& rasterName );
58+
~QgsRasterCalcNode();
59+
60+
Type type() const { return mType; }
61+
62+
//set left node
63+
void setLeft( QgsRasterCalcNode* left ) { delete mLeft; mLeft = left; }
64+
void setRight( QgsRasterCalcNode* right ) { delete mRight; mRight = right; }
65+
66+
/**Calculates result (might be real matrix or single number)*/
67+
bool calculate( QMap<QString, QgsRasterMatrix*>& rasterData, QgsRasterMatrix& result ) const;
68+
69+
private:
70+
Type mType;
71+
QgsRasterCalcNode* mLeft;
72+
QgsRasterCalcNode* mRight;
73+
QgsRasterMatrix* mRasterMatrix;
74+
QString mRasterName;
75+
double mNumber;
76+
Operator mOperator;
77+
};
78+
79+
#endif // QGSRASTERCALCNODE_H

‎src/analysis/raster/qgsrastercalculator.cpp

Lines changed: 373 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,373 @@
1+
/***************************************************************************
2+
qgsrastercalculator.cpp - description
3+
-----------------------
4+
begin : September 28th, 2010
5+
copyright : (C) 2010 by Marco Hugentobler
6+
email : marco dot hugentobler at sourcepole dot ch
7+
***************************************************************************/
8+
9+
/***************************************************************************
10+
* *
11+
* This program is free software; you can redistribute it and/or modify *
12+
* it under the terms of the GNU General Public License as published by *
13+
* the Free Software Foundation; either version 2 of the License, or *
14+
* (at your option) any later version. *
15+
* *
16+
***************************************************************************/
17+
18+
#include "qgsrastercalculator.h"
19+
#include "qgsrastercalcnode.h"
20+
#include "qgsrasterlayer.h"
21+
#include "qgsrastermatrix.h"
22+
#include "cpl_string.h"
23+
#include <QProgressDialog>
24+
25+
extern QgsRasterCalcNode* parseRasterCalcString( const QString& str, QString& parserErrorMsg );
26+
27+
QgsRasterCalculator::QgsRasterCalculator( const QString& formulaString, const QString& outputFile, const QString& outputFormat,
28+
const QgsRectangle& outputExtent, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry>& rasterEntries ): mFormulaString( formulaString ), mOutputFile( outputFile ), mOutputFormat( outputFormat ),
29+
mOutputRectangle( outputExtent ), mNumOutputColumns( nOutputColumns ), mNumOutputRows( nOutputRows ), mRasterEntries( rasterEntries )
30+
{
31+
}
32+
33+
QgsRasterCalculator::~QgsRasterCalculator()
34+
{
35+
}
36+
37+
int QgsRasterCalculator::processCalculation( QProgressDialog* p )
38+
{
39+
//prepare search string / tree
40+
QString errorString;
41+
QgsRasterCalcNode* calcNode = parseRasterCalcString( mFormulaString, errorString );
42+
if( !calcNode )
43+
{
44+
//error
45+
}
46+
47+
double targetGeoTransform[6];
48+
outputGeoTransform( targetGeoTransform );
49+
50+
//open all input rasters for reading
51+
QMap< QString, GDALRasterBandH > mInputRasterBands; //raster references and corresponding scanline data
52+
QMap< QString, QgsRasterMatrix* > inputScanLineData; //stores raster references and corresponding scanline data
53+
QVector< GDALDatasetH > mInputDatasets; //raster references and corresponding dataset
54+
55+
QVector<QgsRasterCalculatorEntry>::const_iterator it = mRasterEntries.constBegin();
56+
for( ; it != mRasterEntries.constEnd(); ++it )
57+
{
58+
if( !it->raster ) // no raster layer in entry
59+
{
60+
return 2;
61+
}
62+
GDALDatasetH inputDataset = GDALOpen( it->raster->source().toLocal8Bit().data(), GA_ReadOnly );
63+
if( inputDataset == NULL )
64+
{
65+
return 2;
66+
}
67+
GDALRasterBandH inputRasterBand = GDALGetRasterBand( inputDataset, it->bandNumber );
68+
if( inputRasterBand == NULL )
69+
{
70+
return 2;
71+
}
72+
73+
mInputDatasets.push_back( inputDataset );
74+
mInputRasterBands.insert( it->ref, inputRasterBand );
75+
inputScanLineData.insert( it->ref, new QgsRasterMatrix( mNumOutputColumns, 1, new float[mNumOutputColumns] ) );
76+
}
77+
78+
//open output dataset for writing
79+
GDALDriverH outputDriver = openOutputDriver();
80+
if( outputDriver == NULL )
81+
{
82+
return 1;
83+
}
84+
GDALDatasetH outputDataset = openOutputFile( outputDriver );
85+
GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset, 1 );
86+
float* resultScanLine = ( float * ) CPLMalloc( sizeof( float ) * mNumOutputColumns );
87+
88+
if( p )
89+
{
90+
p->setMaximum( mNumOutputRows );
91+
}
92+
93+
QgsRasterMatrix resultMatrix;
94+
95+
//read / write line by line
96+
for( int i = 0; i < mNumOutputRows; ++i )
97+
{
98+
if( p )
99+
{
100+
p->setValue( i );
101+
}
102+
103+
if( p && p->wasCanceled() )
104+
{
105+
break;
106+
}
107+
108+
//fill buffers
109+
QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin();
110+
for( ; bufferIt != inputScanLineData.end(); ++bufferIt )
111+
{
112+
double sourceTransformation[6];
113+
GDALRasterBandH sourceRasterBand = mInputRasterBands[bufferIt.key()];
114+
GDALGetGeoTransform( GDALGetBandDataset( sourceRasterBand ), sourceTransformation );
115+
//the function readRasterPart calls GDALRasterIO (and ev. does some conversion if raster transformations are not the same)
116+
readRasterPart( targetGeoTransform, 0, i, mNumOutputColumns, 1, sourceTransformation, sourceRasterBand, bufferIt.value()->data() );
117+
}
118+
119+
if( calcNode->calculate( inputScanLineData, resultMatrix ) )
120+
{
121+
//write scanline to the dataset
122+
if( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, resultMatrix.data(), mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
123+
{
124+
qWarning( "RasterIO error!" );
125+
}
126+
}
127+
128+
}
129+
130+
if( p )
131+
{
132+
p->setValue( mNumOutputRows );
133+
}
134+
135+
//close datasets and release memory
136+
delete calcNode;
137+
QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin();
138+
for( ; bufferIt != inputScanLineData.end(); ++bufferIt )
139+
{
140+
delete bufferIt.value();
141+
}
142+
inputScanLineData.clear();
143+
144+
QVector< GDALDatasetH >::iterator datasetIt = mInputDatasets.begin();
145+
for( ; datasetIt != mInputDatasets.end(); ++ datasetIt )
146+
{
147+
GDALClose( *datasetIt );
148+
}
149+
150+
if( p && p->wasCanceled() )
151+
{
152+
//delete the dataset without closing (because it is faster)
153+
GDALDeleteDataset( outputDriver, mOutputFile.toLocal8Bit().data() );
154+
return 3;
155+
}
156+
GDALClose( outputDataset );
157+
CPLFree( resultScanLine );
158+
return 0;
159+
}
160+
161+
QgsRasterCalculator::QgsRasterCalculator()
162+
{
163+
}
164+
165+
GDALDriverH QgsRasterCalculator::openOutputDriver()
166+
{
167+
char **driverMetadata;
168+
169+
//open driver
170+
GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
171+
172+
if( outputDriver == NULL )
173+
{
174+
return outputDriver; //return NULL, driver does not exist
175+
}
176+
177+
driverMetadata = GDALGetMetadata( outputDriver, NULL );
178+
if( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
179+
{
180+
return NULL; //driver exist, but it does not support the create operation
181+
}
182+
183+
return outputDriver;
184+
}
185+
186+
GDALDatasetH QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver )
187+
{
188+
//open output file
189+
char **papszOptions = NULL;
190+
GDALDatasetH outputDataset = GDALCreate( outputDriver, mOutputFile.toLocal8Bit().data(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions );
191+
if( outputDataset == NULL )
192+
{
193+
return outputDataset;
194+
}
195+
196+
//assign georef information
197+
double geotransform[6];
198+
outputGeoTransform( geotransform );
199+
GDALSetGeoTransform( outputDataset, geotransform );
200+
201+
return outputDataset;
202+
}
203+
204+
void QgsRasterCalculator::readRasterPart( double* targetGeotransform, int xOffset, int yOffset, int nCols, int nRows, double* sourceTransform, GDALRasterBandH sourceBand, float* rasterBuffer )
205+
{
206+
//If dataset transform is the same as the requested transform, do a normal GDAL raster io
207+
if( transformationsEqual( targetGeotransform, sourceTransform ) )
208+
{
209+
GDALRasterIO( sourceBand, GF_Read, xOffset, yOffset, nCols, nRows, rasterBuffer, nCols, nRows, GDT_Float32, 0, 0 );
210+
return;
211+
}
212+
213+
//pixel calculation needed because of different raster position / resolution
214+
int nodataSuccess;
215+
double nodataValue = GDALGetRasterNoDataValue( sourceBand, &nodataSuccess );
216+
QgsRectangle targetRect( targetGeotransform[0] + targetGeotransform[1] * xOffset, targetGeotransform[3] + yOffset * targetGeotransform[5] + nRows * targetGeotransform[5]
217+
, targetGeotransform[0] + targetGeotransform[1] * xOffset + targetGeotransform[1] * nCols, targetGeotransform[3] + yOffset * targetGeotransform[5] );
218+
QgsRectangle sourceRect( sourceTransform[0], sourceTransform[3] + GDALGetRasterBandXSize( sourceBand ) * sourceTransform[5], sourceTransform[0] + GDALGetRasterBandXSize( sourceBand )* sourceTransform[1], sourceTransform[3] );
219+
QgsRectangle intersection = targetRect.intersect( &sourceRect );
220+
221+
//no intersection, fill all the pixels with nodata values
222+
if( intersection.isEmpty() )
223+
{
224+
int nPixels = nCols * nRows;
225+
for( int i = 0; i < nPixels; ++i )
226+
{
227+
rasterBuffer[i] = nodataValue;
228+
}
229+
return;
230+
}
231+
232+
//do raster io in source resolution
233+
int sourcePixelOffsetXMin = floor(( intersection.xMinimum() - sourceTransform[0] ) / sourceTransform[1] );
234+
int sourcePixelOffsetXMax = ceil(( intersection.xMaximum() - sourceTransform[0] ) / sourceTransform[1] );
235+
int nSourcePixelsX = sourcePixelOffsetXMax - sourcePixelOffsetXMin;
236+
int sourcePixelOffsetYMax = floor(( intersection.yMaximum() - sourceTransform[3] ) / sourceTransform[5] );
237+
int sourcePixelOffsetYMin = ceil(( intersection.yMinimum() - sourceTransform[3] ) / sourceTransform[5] );
238+
int nSourcePixelsY = sourcePixelOffsetYMin - sourcePixelOffsetYMax;
239+
float* sourceRaster = ( float * ) CPLMalloc( sizeof( float ) * nSourcePixelsX * nSourcePixelsY );
240+
double sourceRasterXMin = sourceRect.xMinimum() + sourcePixelOffsetXMin * sourceTransform[1];
241+
double sourceRasterYMax = sourceRect.yMaximum() + sourcePixelOffsetYMax * sourceTransform[5];
242+
GDALRasterIO( sourceBand, GF_Read, sourcePixelOffsetXMin, sourcePixelOffsetYMax, nSourcePixelsX, nSourcePixelsY,
243+
sourceRaster, nSourcePixelsX, nSourcePixelsY, GDT_Float32, 0, 0 );
244+
245+
246+
double targetPixelX;
247+
double targetPixelXMin = targetGeotransform[0] + targetGeotransform[1] * xOffset + targetGeotransform[1] / 2.0;
248+
double targetPixelY = targetGeotransform[3] + targetGeotransform[5] * yOffset + targetGeotransform[5] / 2.0; //coordinates of current target pixel
249+
int sourceIndexX, sourceIndexY; //current raster index in source pixels
250+
double sx, sy;
251+
for( int i = 0; i < nRows; ++i )
252+
{
253+
targetPixelX = targetPixelXMin;
254+
for( int j = 0; j < nCols; ++j )
255+
{
256+
sx = ( targetPixelX - sourceRasterXMin ) / sourceTransform[1];
257+
sourceIndexX = sx > 0 ? sx : floor( sx );
258+
sy = ( targetPixelY - sourceRasterYMax ) / sourceTransform[5];
259+
sourceIndexY = sy > 0 ? sy : floor( sy );
260+
if( sourceIndexX >= 0 && sourceIndexX < nSourcePixelsX
261+
&& sourceIndexY >= 0 && sourceIndexY < nSourcePixelsY )
262+
{
263+
rasterBuffer[j + i*nRows] = sourceRaster[ sourceIndexX + nSourcePixelsX * sourceIndexY ];
264+
}
265+
else
266+
{
267+
rasterBuffer[j + i*j] = nodataValue;
268+
}
269+
targetPixelX += targetGeotransform[1];
270+
}
271+
targetPixelY += targetGeotransform[5];
272+
}
273+
274+
CPLFree( sourceRaster );
275+
return;
276+
277+
#if 0
278+
//If dataset transform is the same as the requested transform, do a normal GDAL raster io
279+
if( transformationsEqual( targetGeotransform, sourceTransform ) )
280+
{
281+
GDALRasterIO( sourceBand, GF_Read, xOffset, yOffset, nCols, nRows, rasterBuffer, nCols, nRows, GDT_Float32, 0, 0 );
282+
return;
283+
}
284+
285+
//pixel calculation needed because of different raster position / resolution
286+
287+
//calculate raster span in source layer to fetch (which is a bit larger than the target area to fetch if the two resolutions don't match)
288+
//calculate offset
289+
290+
//QgsRectangle targetRect( targetGeotransform[0], targetGeotransform[3] + mNumOutputColumns * targetGeotransform[5], targetGeotransform[0] + mNumOutputRows * targetGeotransform[1], targetGeotransform[3]);
291+
QgsRectangle targetRect( targetGeotransform[0] + targetGeotransform[1] * xOffset, targetGeotransform[3] + yOffset * targetGeotransform[5] + nRows * targetGeotransform[5]
292+
, targetGeotransform[0] + targetGeotransform[1] * xOffset + targetGeotransform[1] * nCols, targetGeotransform[3] + yOffset * targetGeotransform[5] );
293+
QgsRectangle sourceRect( sourceTransform[0], sourceTransform[3] + GDALGetRasterBandXSize( sourceBand ) * sourceTransform[5], sourceTransform[0] + GDALGetRasterBandYSize( sourceBand )* sourceTransform[1], sourceTransform[3] );
294+
QgsRectangle intersection = targetRect.intersect( &sourceRect );
295+
296+
//no intersection, fill all the pixels with nodata values
297+
if( intersection.isEmpty() )
298+
{
299+
int nPixels = nCols * nRows;
300+
for( int i = 0; i < nPixels; ++i )
301+
{
302+
rasterBuffer[i] = 0;
303+
}
304+
return;
305+
}
306+
307+
//do raster io in source resolution
308+
int sourcePixelOffsetX = ( intersection.xMinimum() - sourceTransform[0] ) / sourceTransform[1];
309+
int sourcePixelOffsetY = ( intersection.yMaximum() - sourceTransform[3] ) / sourceTransform[5];
310+
int nSourcePixelsX = ceil( fabs( intersection.width() / sourceTransform[1] ) );
311+
int nSourcePixelsY = ceil( fabs( intersection.height() / sourceTransform[5] ) );
312+
float* sourceRaster = ( float * ) CPLMalloc( sizeof( float ) * nSourcePixelsX * nSourcePixelsY );
313+
314+
double sourceRasterXMin = sourceRect.xMinimum() + sourcePixelOffsetX * sourceTransform[1];
315+
double sourceRasterYMax = sourceRect.yMaximum() + sourcePixelOffsetY * sourceTransform[5];
316+
317+
GDALRasterIO( sourceBand, GF_Read, sourcePixelOffsetX, sourcePixelOffsetY, nSourcePixelsX, nSourcePixelsY,
318+
sourceRaster, nSourcePixelsX, nSourcePixelsY, GDT_Float32, 0, 0 );
319+
320+
321+
double targetPixelX;
322+
double targetPixelXMin = targetGeotransform[0] + targetGeotransform[1] * xOffset + targetGeotransform[1] / 2.0;
323+
double targetPixelY = targetGeotransform[3] + targetGeotransform[5] * yOffset + targetGeotransform[5] / 2.0; //coordinates of current target pixel
324+
int sourceIndexX, sourceIndexY; //current raster index in source pixels
325+
for( int i = 0; i < nRows; ++i )
326+
{
327+
targetPixelX = targetPixelXMin;
328+
for( int j = 0; j < nCols; ++j )
329+
{
330+
sourceIndexX = ( targetPixelX - sourceRasterXMin ) / sourceTransform[1];
331+
sourceIndexY = ( targetPixelY - sourceRasterYMax ) / sourceTransform[5];
332+
if( sourceIndexX >= 0 && sourceIndexX < nSourcePixelsX
333+
&& sourceIndexY >= 0 && sourceIndexY < nSourcePixelsY )
334+
{
335+
rasterBuffer[j + i*nCols] = sourceRaster[ sourceIndexX + nSourcePixelsX * sourceIndexY ];
336+
}
337+
else
338+
{
339+
rasterBuffer[j + i*nCols] = 0; //todo: insert null value of rasterband here
340+
}
341+
targetPixelX += targetGeotransform[1];
342+
}
343+
targetPixelY += targetGeotransform[5];
344+
}
345+
346+
CPLFree( sourceRaster );
347+
return;
348+
#endif //0
349+
}
350+
351+
bool QgsRasterCalculator::transformationsEqual( double* t1, double* t2 ) const
352+
{
353+
for( int i = 0; i < 6; ++i )
354+
{
355+
if( !doubleNear( t1[i], t2[i] ) )
356+
{
357+
return false;
358+
}
359+
}
360+
return true;
361+
}
362+
363+
void QgsRasterCalculator::outputGeoTransform( double* transform ) const
364+
{
365+
transform[0] = mOutputRectangle.xMinimum();
366+
transform[1] = mOutputRectangle.width() / mNumOutputColumns;
367+
transform[2] = 0;
368+
transform[3] = mOutputRectangle.yMaximum();
369+
transform[4] = 0;
370+
transform[5] = -mOutputRectangle.height() / mNumOutputRows;
371+
}
372+
373+

‎src/analysis/raster/qgsrastercalculator.h

Lines changed: 89 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,89 @@
1+
/***************************************************************************
2+
qgsrastercalculator.h - description
3+
---------------------
4+
begin : September 28th, 2010
5+
copyright : (C) 2010 by Marco Hugentobler
6+
email : marco dot hugentobler at sourcepole dot ch
7+
***************************************************************************/
8+
9+
/***************************************************************************
10+
* *
11+
* This program is free software; you can redistribute it and/or modify *
12+
* it under the terms of the GNU General Public License as published by *
13+
* the Free Software Foundation; either version 2 of the License, or *
14+
* (at your option) any later version. *
15+
* *
16+
***************************************************************************/
17+
18+
#ifndef QGSRASTERCALCULATOR_H
19+
#define QGSRASTERCALCULATOR_H
20+
21+
#include "qgsfield.h"
22+
#include "qgsrectangle.h"
23+
#include <QString>
24+
#include <QVector>
25+
#include "gdal.h"
26+
27+
class QgsRasterLayer;
28+
class QProgressDialog;
29+
30+
31+
struct ANALYSIS_EXPORT QgsRasterCalculatorEntry
32+
{
33+
QString ref; //name
34+
QgsRasterLayer* raster; //pointer to rasterlayer
35+
int bandNumber; //raster band number
36+
};
37+
38+
/**Raster calculator class*/
39+
class ANALYSIS_EXPORT QgsRasterCalculator
40+
{
41+
public:
42+
QgsRasterCalculator( const QString& formulaString, const QString& outputFile, const QString& outputFormat,
43+
const QgsRectangle& outputExtent, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry>& rasterEntries );
44+
~QgsRasterCalculator();
45+
46+
/**Starts the calculation and writes new raster
47+
@param p progress bar (or 0 if called from non-gui code)
48+
@return 0 in case of success*/
49+
int processCalculation( QProgressDialog* p = 0 );
50+
51+
private:
52+
//default constructor forbidden. We need formula, output file, output format and output raster resolution obligatory
53+
QgsRasterCalculator();
54+
55+
/**Opens the output driver and tests if it supports the creation of a new dataset
56+
@return NULL on error and the driver handle on success*/
57+
GDALDriverH openOutputDriver();
58+
59+
/**Opens the output file and sets the same geotransform and CRS as the input data
60+
@return the output dataset or NULL in case of error*/
61+
GDALDatasetH openOutputFile( GDALDriverH outputDriver );
62+
63+
/**Reads raster pixels from a dataset/band
64+
@param targetGeotransformation transformation parameters of the requested raster array (not necessarily the same as the transform of the source dataset */
65+
void readRasterPart( double* targetGeotransform, int xOffset, int yOffset, int nCols, int nRows, double* sourceTransform, GDALRasterBandH sourceBand, float* rasterBuffer );
66+
67+
/**Compares two geotransformations (six parameter double arrays*/
68+
bool transformationsEqual( double* t1, double* t2 ) const;
69+
70+
/**Sets gdal 6 parameters array from mOutputRectangle, mNumOutputColumns, mNumOutputRows
71+
@param transform double[6] array that receives the GDAL parameters*/
72+
void outputGeoTransform( double* transform ) const;
73+
74+
QString mFormulaString;
75+
QString mOutputFile;
76+
QString mOutputFormat;
77+
78+
/**Output raster extent*/
79+
QgsRectangle mOutputRectangle;
80+
/**Number of output columns*/
81+
int mNumOutputColumns;
82+
/**Number of output rows*/
83+
int mNumOutputRows;
84+
85+
/***/
86+
QVector<QgsRasterCalculatorEntry> mRasterEntries;
87+
};
88+
89+
#endif // QGSRASTERCALCULATOR_H

‎src/analysis/raster/qgsrastermatrix.cpp

Lines changed: 423 additions & 0 deletions
Large diffs are not rendered by default.

‎src/analysis/raster/qgsrastermatrix.h

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
/***************************************************************************
2+
qgsrastermatrix.h
3+
-----------------
4+
begin : 2010-10-23
5+
copyright : (C) 20010 by Marco Hugentobler
6+
email : marco dot hugentobler at sourcepole dot ch
7+
***************************************************************************/
8+
9+
/***************************************************************************
10+
* *
11+
* This program is free software; you can redistribute it and/or modify *
12+
* it under the terms of the GNU General Public License as published by *
13+
* the Free Software Foundation; either version 2 of the License, or *
14+
* (at your option) any later version. *
15+
* *
16+
***************************************************************************/
17+
18+
#ifndef QGSRASTERMATRIX_H
19+
#define QGSRASTERMATRIX_H
20+
21+
class ANALYSIS_EXPORT QgsRasterMatrix
22+
{
23+
public:
24+
25+
enum TwoArgOperator
26+
{
27+
opPLUS,
28+
opMINUS,
29+
opMUL,
30+
opDIV,
31+
opPOW,
32+
};
33+
34+
enum OneArgOperator
35+
{
36+
opSQRT,
37+
opSIN,
38+
opCOS,
39+
opTAN,
40+
opASIN,
41+
opACOS,
42+
opATAN
43+
};
44+
45+
/**Takes ownership of data array*/
46+
QgsRasterMatrix();
47+
QgsRasterMatrix( int nCols, int nRows, float* data );
48+
QgsRasterMatrix( const QgsRasterMatrix& m );
49+
~QgsRasterMatrix();
50+
51+
/**Returns true if matrix is 1x1 (=scalar number)*/
52+
bool isNumber() const { return ( mColumns == 1 && mRows == 1 ); }
53+
double number() const { return mData[0]; }
54+
55+
/**Returns data array (but not ownership)*/
56+
float* data() { return mData; }
57+
/**Returns data and ownership. Sets data and nrows, ncols of this matrix to 0*/
58+
float* takeData();
59+
60+
void setData( int cols, int rows, float* data );
61+
62+
int nColumns() const { return mColumns; }
63+
int nRows() const { return mRows; }
64+
65+
QgsRasterMatrix& operator=( const QgsRasterMatrix& m );
66+
/**Adds another matrix to this one*/
67+
bool add( const QgsRasterMatrix& other );
68+
/**Subtracts another matrix from this one*/
69+
bool subtract( const QgsRasterMatrix& other );
70+
bool multiply( const QgsRasterMatrix& other );
71+
bool divide( const QgsRasterMatrix& other );
72+
bool power( const QgsRasterMatrix& other );
73+
74+
bool squareRoot();
75+
bool sinus();
76+
bool asinus();
77+
bool cosinus();
78+
bool acosinus();
79+
bool tangens();
80+
bool atangens();
81+
82+
private:
83+
int mColumns;
84+
int mRows;
85+
float* mData;
86+
87+
/**+,-,*,/,^*/
88+
bool twoArgumentOperation( TwoArgOperator op, const QgsRasterMatrix& other );
89+
bool testPowerValidity( double base, double power );
90+
};
91+
92+
#endif // QGSRASTERMATRIX_H

‎src/app/CMakeLists.txt

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -72,6 +72,7 @@ SET(QGIS_APP_SRCS
7272
qgspluginmetadata.cpp
7373
qgspluginregistry.cpp
7474
qgsprojectproperties.cpp
75+
qgsrastercalcdialog.cpp
7576
qgsrasterlayerproperties.cpp
7677
qgssearchquerybuilder.cpp
7778
qgstextannotationdialog.cpp
@@ -204,6 +205,7 @@ SET (QGIS_APP_MOC_HDRS
204205
qgspastetransformations.h
205206
qgspluginmanager.h
206207
qgsprojectproperties.h
208+
qgsrastercalcdialog.h
207209
qgsrasterlayerproperties.h
208210
qgssearchquerybuilder.h
209211
qgstextannotationdialog.h
@@ -325,6 +327,7 @@ INCLUDE_DIRECTORIES(
325327
${CMAKE_CURRENT_BINARY_DIR}/../ui
326328
${QWT_INCLUDE_DIR}
327329
${QT_QTUITOOLS_INCLUDE_DIR}
330+
../analysis/raster
328331
../core
329332
../core/gps ../core/gps/qextserialport
330333
../core/composer ../core/raster ../core/renderer ../core/symbology ../core/symbology-ng
@@ -376,6 +379,7 @@ TARGET_LINK_LIBRARIES(qgis
376379
${QT_QTMAIN_LIBRARY}
377380
qgis_core
378381
qgis_gui
382+
qgis_analysis
379383
)
380384

381385
IF( WIN32 )

‎src/app/qgisapp.cpp

Lines changed: 456 additions & 428 deletions
Large diffs are not rendered by default.

‎src/app/qgisapp.h

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -520,6 +520,9 @@ class QgisApp : public QMainWindow
520520
void fileNew();
521521
//! As above but allows forcing without prompt
522522
void fileNew( bool thePromptToSaveFlag );
523+
//! Calculate new rasters from existing ones
524+
void showRasterCalculator();
525+
523526
//! Create a new empty vector layer
524527
void newVectorLayer();
525528
#ifdef HAVE_SPATIALITE
@@ -927,6 +930,7 @@ class QgisApp : public QMainWindow
927930

928931
QAction *mActionNewVectorLayer;
929932
QAction *mActionNewSpatialiteLayer;
933+
QAction *mActionShowRasterCalculator;
930934
QAction *mActionAddOgrLayer;
931935
QAction *mActionAddRasterLayer;
932936
QAction *mActionAddPgLayer;

‎src/app/qgsrastercalcdialog.cpp

Lines changed: 391 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,391 @@
1+
/***************************************************************************
2+
qgsrastercalcdialog.h - description
3+
---------------------
4+
begin : September 28th, 2010
5+
copyright : (C) 2010 by Marco Hugentobler
6+
email : marco dot hugentobler at sourcepole dot ch
7+
***************************************************************************/
8+
9+
/***************************************************************************
10+
* *
11+
* This program is free software; you can redistribute it and/or modify *
12+
* it under the terms of the GNU General Public License as published by *
13+
* the Free Software Foundation; either version 2 of the License, or *
14+
* (at your option) any later version. *
15+
* *
16+
***************************************************************************/
17+
18+
#include "qgsrastercalcdialog.h"
19+
#include "qgsmaplayerregistry.h"
20+
#include "qgsrastercalcnode.h"
21+
#include "qgsrasterlayer.h"
22+
#include "cpl_string.h"
23+
#include "gdal.h"
24+
25+
#include <QFileDialog>
26+
#include <QSettings>
27+
28+
extern QgsRasterCalcNode* parseRasterCalcString( const QString& str, QString& parserErrorMsg );
29+
30+
QgsRasterCalcDialog::QgsRasterCalcDialog( QWidget * parent, Qt::WindowFlags f ): QDialog( parent, f )
31+
{
32+
setupUi( this );
33+
34+
//add supported output formats
35+
insertAvailableOutputFormats();
36+
insertAvailableRasterBands();
37+
}
38+
39+
QgsRasterCalcDialog::~QgsRasterCalcDialog()
40+
{
41+
}
42+
43+
QString QgsRasterCalcDialog::formulaString() const
44+
{
45+
return mExpressionTextEdit->toPlainText();
46+
}
47+
48+
QString QgsRasterCalcDialog::outputFile() const
49+
{
50+
QString outputFileName = mOutputLayerLineEdit->text();
51+
QFileInfo fileInfo( outputFileName );
52+
QString suffix = fileInfo.suffix();
53+
if( !suffix.isEmpty() )
54+
{
55+
return outputFileName;
56+
}
57+
58+
//add the file format extension if the user did not specify it
59+
int index = mOutputFormatComboBox->currentIndex();
60+
if( index == -1 )
61+
{
62+
return outputFileName;
63+
}
64+
65+
QString driverShortName = mOutputFormatComboBox->itemData( index ).toString();
66+
QMap<QString, QString>::const_iterator it = mDriverExtensionMap.find( driverShortName );
67+
if( it == mDriverExtensionMap.constEnd() )
68+
{
69+
return outputFileName;
70+
}
71+
72+
return ( outputFileName + "." + it.value() );
73+
}
74+
75+
QString QgsRasterCalcDialog::outputFormat() const
76+
{
77+
int index = mOutputFormatComboBox->currentIndex();
78+
if( index == -1 )
79+
{
80+
return "";
81+
}
82+
return mOutputFormatComboBox->itemData( index ).toString();
83+
}
84+
85+
bool QgsRasterCalcDialog::addLayerToProject() const
86+
{
87+
return mAddResultToProjectCheckBox->isChecked();
88+
}
89+
90+
QVector<QgsRasterCalculatorEntry> QgsRasterCalcDialog::rasterEntries() const
91+
{
92+
QVector<QgsRasterCalculatorEntry> entries;
93+
QString expressionString = mExpressionTextEdit->toPlainText();
94+
95+
QList<QgsRasterCalculatorEntry>::const_iterator bandIt = mAvailableRasterBands.constBegin();
96+
for( ; bandIt != mAvailableRasterBands.constEnd(); ++bandIt )
97+
{
98+
if( expressionString.contains( bandIt->ref ) )
99+
{
100+
entries.push_back( *bandIt );
101+
}
102+
}
103+
104+
return entries;
105+
}
106+
107+
void QgsRasterCalcDialog::insertAvailableRasterBands()
108+
{
109+
QMap<QString, QgsMapLayer*>& layers = QgsMapLayerRegistry::instance()->mapLayers();
110+
QMap<QString, QgsMapLayer*>::const_iterator layerIt = layers.constBegin();
111+
112+
bool firstLayer = true;
113+
for( ; layerIt != layers.constEnd(); ++layerIt )
114+
{
115+
QgsRasterLayer* rlayer = dynamic_cast<QgsRasterLayer*>( layerIt.value() );
116+
if( rlayer && !rlayer->usesProvider() )
117+
{
118+
if( firstLayer ) //set bounding box / resolution of output to the values of the first possible input layer
119+
{
120+
mNColumnsSpinBox->setValue( rlayer->width() );
121+
mNRowsSpinBox->setValue( rlayer->height() );
122+
QgsRectangle bbox = rlayer->extent();
123+
mXMinSpinBox->setValue( bbox.xMinimum() );
124+
mXMaxSpinBox->setValue( bbox.xMaximum() );
125+
mYMinSpinBox->setValue( bbox.yMinimum() );
126+
mYMaxSpinBox->setValue( bbox.yMaximum() );
127+
firstLayer = false;
128+
}
129+
//get number of bands
130+
for( int i = 0; i < rlayer->bandCount(); ++i )
131+
{
132+
QgsRasterCalculatorEntry entry;
133+
entry.raster = rlayer;
134+
entry.bandNumber = i + 1;
135+
entry.ref = rlayer->name() + "@" + QString::number( i + 1 );
136+
mAvailableRasterBands.push_back( entry );
137+
mRasterBandsListWidget->addItem( entry.ref );
138+
}
139+
}
140+
}
141+
}
142+
143+
void QgsRasterCalcDialog::insertAvailableOutputFormats()
144+
{
145+
GDALAllRegister();
146+
147+
int nDrivers = GDALGetDriverCount();
148+
for( int i = 0; i < nDrivers; ++i )
149+
{
150+
GDALDriverH driver = GDALGetDriver( i );
151+
if( driver != NULL )
152+
{
153+
char** driverMetadata = GDALGetMetadata( driver, NULL );
154+
if( CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
155+
{
156+
mOutputFormatComboBox->addItem( GDALGetDriverLongName( driver ), QVariant( GDALGetDriverShortName( driver ) ) );
157+
158+
//store the driver shortnames and the corresponding extensions
159+
//(just in case the user does not give an extension for the output file name)
160+
int index = 0;
161+
while(( driverMetadata ) && driverMetadata[index] != 0 )
162+
{
163+
QStringList metadataTokens = QString( driverMetadata[index] ).split( "=", QString::SkipEmptyParts );
164+
if( metadataTokens.size() < 1 )
165+
{
166+
break;
167+
}
168+
169+
if( metadataTokens[0] == "DMD_EXTENSION" )
170+
{
171+
if( metadataTokens.size() < 2 )
172+
{
173+
++index;
174+
continue;
175+
}
176+
mDriverExtensionMap.insert( QString( GDALGetDriverShortName( driver ) ), metadataTokens[1] );
177+
break;
178+
}
179+
++index;
180+
}
181+
}
182+
}
183+
}
184+
185+
//and set last used driver in combo box
186+
QSettings s;
187+
QString lastUsedDriver = s.value( "/RasterCalculator/lastOutputFormat", "GeoTIFF" ).toString();
188+
int lastDriverIndex = mOutputFormatComboBox->findText( lastUsedDriver );
189+
if( lastDriverIndex != -1 )
190+
{
191+
mOutputFormatComboBox->setCurrentIndex( lastDriverIndex );
192+
}
193+
}
194+
195+
QgsRectangle QgsRasterCalcDialog::outputRectangle() const
196+
{
197+
return QgsRectangle( mXMinSpinBox->value(), mYMinSpinBox->value(), mXMaxSpinBox->value(), mYMaxSpinBox->value() );
198+
}
199+
200+
int QgsRasterCalcDialog::numberOfColumns() const
201+
{
202+
return mNColumnsSpinBox->value();
203+
}
204+
205+
int QgsRasterCalcDialog::numberOfRows() const
206+
{
207+
return mNRowsSpinBox->value();
208+
}
209+
210+
//slots
211+
212+
void QgsRasterCalcDialog::on_mButtonBox_accepted()
213+
{
214+
//save last output format
215+
QSettings s;
216+
s.setValue( "/RasterCalculator/lastOutputFormat", QVariant( mOutputFormatComboBox->currentText() ) );
217+
}
218+
219+
void QgsRasterCalcDialog::on_mOutputLayerPushButton_clicked()
220+
{
221+
QString saveFileName = QFileDialog::getSaveFileName( 0, tr( "Enter result file" ) );
222+
if( !saveFileName.isNull() )
223+
{
224+
mOutputLayerLineEdit->setText( saveFileName );
225+
}
226+
}
227+
228+
void QgsRasterCalcDialog::on_mCurrentLayerExtentButton_clicked()
229+
{
230+
QListWidgetItem* currentLayerItem = mRasterBandsListWidget->currentItem();
231+
if( currentLayerItem )
232+
{
233+
QgsRasterLayer* rlayer = 0;
234+
QList<QgsRasterCalculatorEntry>::const_iterator rasterIt = mAvailableRasterBands.constBegin();
235+
for( ; rasterIt != mAvailableRasterBands.constEnd(); ++rasterIt )
236+
{
237+
if( rasterIt->ref == currentLayerItem->text() )
238+
{
239+
rlayer = rasterIt->raster;
240+
}
241+
}
242+
243+
if( !rlayer )
244+
{
245+
return;
246+
}
247+
248+
QgsRectangle layerExtent = rlayer->extent();
249+
mXMinSpinBox->setValue( layerExtent.xMinimum() );
250+
mXMaxSpinBox->setValue( layerExtent.xMaximum() );
251+
mYMinSpinBox->setValue( layerExtent.yMinimum() );
252+
mYMaxSpinBox->setValue( layerExtent.yMaximum() );
253+
mNColumnsSpinBox->setValue( rlayer->width() );
254+
mNRowsSpinBox->setValue( rlayer->height() );
255+
}
256+
}
257+
258+
void QgsRasterCalcDialog::on_mExpressionTextEdit_textChanged()
259+
{
260+
if( expressionValid() )
261+
{
262+
mExpressionValidLabel->setText( tr( "Expression valid" ) );
263+
if( filePathValid() )
264+
{
265+
mButtonBox->button( QDialogButtonBox::Ok )->setEnabled( true );
266+
return;
267+
}
268+
}
269+
else
270+
{
271+
mExpressionValidLabel->setText( tr( "Expression invalid" ) );
272+
}
273+
mButtonBox->button( QDialogButtonBox::Ok )->setEnabled( false );
274+
}
275+
276+
void QgsRasterCalcDialog::on_mOutputLayerLineEdit_textChanged( const QString& text )
277+
{
278+
setAcceptButtonState();
279+
}
280+
281+
void QgsRasterCalcDialog::setAcceptButtonState()
282+
{
283+
if( expressionValid() && filePathValid() )
284+
{
285+
mButtonBox->button( QDialogButtonBox::Ok )->setEnabled( true );
286+
}
287+
else
288+
{
289+
mButtonBox->button( QDialogButtonBox::Ok )->setEnabled( false );
290+
}
291+
}
292+
293+
bool QgsRasterCalcDialog::expressionValid() const
294+
{
295+
QString errorString;
296+
QgsRasterCalcNode* testNode = parseRasterCalcString( mExpressionTextEdit->toPlainText(), errorString );
297+
if( testNode )
298+
{
299+
delete testNode;
300+
return true;
301+
}
302+
return false;
303+
}
304+
305+
bool QgsRasterCalcDialog::filePathValid() const
306+
{
307+
QString outputPath = QFileInfo( mOutputLayerLineEdit->text() ).absolutePath();
308+
if( QFileInfo( outputPath ).isWritable() )
309+
{
310+
return true;
311+
}
312+
else
313+
{
314+
return false;
315+
}
316+
}
317+
318+
void QgsRasterCalcDialog::on_mRasterBandsListWidget_itemDoubleClicked( QListWidgetItem* item )
319+
{
320+
mExpressionTextEdit->insertPlainText( item->text() );
321+
}
322+
323+
void QgsRasterCalcDialog::on_mPlusPushButton_clicked()
324+
{
325+
mExpressionTextEdit->insertPlainText( " + " );
326+
}
327+
328+
void QgsRasterCalcDialog::on_mMinusPushButton_clicked()
329+
{
330+
mExpressionTextEdit->insertPlainText( " - " );
331+
}
332+
333+
void QgsRasterCalcDialog::on_mMultiplyPushButton_clicked()
334+
{
335+
mExpressionTextEdit->insertPlainText( " * " );
336+
}
337+
338+
void QgsRasterCalcDialog::on_mDividePushButton_clicked()
339+
{
340+
mExpressionTextEdit->insertPlainText( " / " );
341+
}
342+
343+
void QgsRasterCalcDialog::on_mSqrtButton_clicked()
344+
{
345+
mExpressionTextEdit->insertPlainText( " sqrt ( " );
346+
}
347+
348+
void QgsRasterCalcDialog::on_mCosButton_clicked()
349+
{
350+
mExpressionTextEdit->insertPlainText( " cos ( " );
351+
}
352+
353+
void QgsRasterCalcDialog::on_mSinButton_clicked()
354+
{
355+
mExpressionTextEdit->insertPlainText( " sin ( " );
356+
}
357+
358+
void QgsRasterCalcDialog::on_mASinButton_clicked()
359+
{
360+
mExpressionTextEdit->insertPlainText( " asin ( " );
361+
}
362+
363+
void QgsRasterCalcDialog::on_mExpButton_clicked()
364+
{
365+
mExpressionTextEdit->insertPlainText( " ^ " );
366+
}
367+
368+
void QgsRasterCalcDialog::on_mTanButton_clicked()
369+
{
370+
mExpressionTextEdit->insertPlainText( " tan ( " );
371+
}
372+
373+
void QgsRasterCalcDialog::on_mACosButton_clicked()
374+
{
375+
mExpressionTextEdit->insertPlainText( " acos ( " );
376+
}
377+
378+
void QgsRasterCalcDialog::on_mATanButton_clicked()
379+
{
380+
mExpressionTextEdit->insertPlainText( " atan ( " );
381+
}
382+
383+
void QgsRasterCalcDialog::on_mOpenBracketPushButton_clicked()
384+
{
385+
mExpressionTextEdit->insertPlainText( " ( " );
386+
}
387+
388+
void QgsRasterCalcDialog::on_mCloseBracketPushButton_clicked()
389+
{
390+
mExpressionTextEdit->insertPlainText( " ) " );
391+
}

‎src/app/qgsrastercalcdialog.h

Lines changed: 89 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,89 @@
1+
/***************************************************************************
2+
qgsrastercalcdialog.h - description
3+
---------------------
4+
begin : September 28th, 2010
5+
copyright : (C) 2010 by Marco Hugentobler
6+
email : marco dot hugentobler at sourcepole dot ch
7+
***************************************************************************/
8+
9+
/***************************************************************************
10+
* *
11+
* This program is free software; you can redistribute it and/or modify *
12+
* it under the terms of the GNU General Public License as published by *
13+
* the Free Software Foundation; either version 2 of the License, or *
14+
* (at your option) any later version. *
15+
* *
16+
***************************************************************************/
17+
18+
#ifndef QGSRASTERCALCDIALOG_H
19+
#define QGSRASTERCALCDIALOG_H
20+
21+
#include "ui_qgsrastercalcdialogbase.h"
22+
#include "qgsrastercalculator.h"
23+
24+
/**A dialog to enter a raster calculation expression*/
25+
class QgsRasterCalcDialog: public QDialog, private Ui::QgsRasterCalcDialogBase
26+
{
27+
Q_OBJECT
28+
public:
29+
QgsRasterCalcDialog( QWidget * parent = 0, Qt::WindowFlags f = 0 );
30+
~QgsRasterCalcDialog();
31+
32+
QString formulaString() const;
33+
QString outputFile() const;
34+
QString outputFormat() const;
35+
bool addLayerToProject() const;
36+
37+
/**Bounding box for output raster*/
38+
QgsRectangle outputRectangle() const;
39+
/**Number of pixels in x-direction*/
40+
int numberOfColumns() const;
41+
/**Number of pixels in y-direction*/
42+
int numberOfRows() const;
43+
44+
QVector<QgsRasterCalculatorEntry> rasterEntries() const;
45+
46+
private slots:
47+
void on_mOutputLayerPushButton_clicked();
48+
void on_mRasterBandsListWidget_itemDoubleClicked( QListWidgetItem* item );
49+
void on_mButtonBox_accepted();
50+
void on_mCurrentLayerExtentButton_clicked();
51+
void on_mExpressionTextEdit_textChanged();
52+
void on_mOutputLayerLineEdit_textChanged( const QString& text );
53+
/**Enables ok button if calculator expression is valid and output file path exists*/
54+
void setAcceptButtonState();
55+
56+
//calculator buttons
57+
void on_mPlusPushButton_clicked();
58+
void on_mMinusPushButton_clicked();
59+
void on_mMultiplyPushButton_clicked();
60+
void on_mDividePushButton_clicked();
61+
void on_mSqrtButton_clicked();
62+
void on_mCosButton_clicked();
63+
void on_mSinButton_clicked();
64+
void on_mASinButton_clicked();
65+
void on_mExpButton_clicked();
66+
void on_mTanButton_clicked();
67+
void on_mACosButton_clicked();
68+
void on_mATanButton_clicked();
69+
void on_mOpenBracketPushButton_clicked();
70+
void on_mCloseBracketPushButton_clicked();
71+
72+
private:
73+
//insert available GDAL drivers that support the create() option
74+
void insertAvailableOutputFormats();
75+
/**Accesses the available raster layers/bands from the layer registry*/
76+
void insertAvailableRasterBands();
77+
78+
/**Returns true if raster calculator expression has valid syntax*/
79+
bool expressionValid() const;
80+
/**Returns true if output file directory exists*/
81+
bool filePathValid() const;
82+
83+
/**Stores relation between driver name and extension*/
84+
QMap<QString, QString> mDriverExtensionMap;
85+
86+
QList<QgsRasterCalculatorEntry> mAvailableRasterBands;
87+
};
88+
89+
#endif // QGSRASTERCALCDIALOG_H

‎src/ui/qgsrastercalcdialogbase.ui

Lines changed: 381 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,381 @@
1+
<?xml version="1.0" encoding="UTF-8"?>
2+
<ui version="4.0">
3+
<class>QgsRasterCalcDialogBase</class>
4+
<widget class="QDialog" name="QgsRasterCalcDialogBase">
5+
<property name="geometry">
6+
<rect>
7+
<x>0</x>
8+
<y>0</y>
9+
<width>651</width>
10+
<height>512</height>
11+
</rect>
12+
</property>
13+
<property name="windowTitle">
14+
<string>Dialog</string>
15+
</property>
16+
<layout class="QGridLayout" name="gridLayout_5">
17+
<item row="0" column="0" colspan="2">
18+
<widget class="QGroupBox" name="mRasterBandsGroupBox">
19+
<property name="title">
20+
<string>Raster bands</string>
21+
</property>
22+
<layout class="QGridLayout" name="gridLayout_2">
23+
<item row="0" column="0">
24+
<widget class="QListWidget" name="mRasterBandsListWidget"/>
25+
</item>
26+
</layout>
27+
</widget>
28+
</item>
29+
<item row="0" column="2">
30+
<widget class="QGroupBox" name="mResultGroupBox">
31+
<property name="title">
32+
<string>Result layer</string>
33+
</property>
34+
<layout class="QGridLayout" name="gridLayout_4">
35+
<item row="0" column="0">
36+
<widget class="QLabel" name="mOutputLayerLabel">
37+
<property name="text">
38+
<string>Output layer</string>
39+
</property>
40+
<property name="buddy">
41+
<cstring>mOutputLayerPushButton</cstring>
42+
</property>
43+
</widget>
44+
</item>
45+
<item row="0" column="2">
46+
<widget class="QLineEdit" name="mOutputLayerLineEdit"/>
47+
</item>
48+
<item row="0" column="3">
49+
<widget class="QPushButton" name="mOutputLayerPushButton">
50+
<property name="sizePolicy">
51+
<sizepolicy hsizetype="Preferred" vsizetype="Fixed">
52+
<horstretch>0</horstretch>
53+
<verstretch>0</verstretch>
54+
</sizepolicy>
55+
</property>
56+
<property name="minimumSize">
57+
<size>
58+
<width>20</width>
59+
<height>0</height>
60+
</size>
61+
</property>
62+
<property name="text">
63+
<string>...</string>
64+
</property>
65+
</widget>
66+
</item>
67+
<item row="1" column="0">
68+
<widget class="QPushButton" name="mCurrentLayerExtentButton">
69+
<property name="text">
70+
<string>Current layer extent</string>
71+
</property>
72+
</widget>
73+
</item>
74+
<item row="1" column="1" colspan="3">
75+
<spacer name="horizontalSpacer">
76+
<property name="orientation">
77+
<enum>Qt::Horizontal</enum>
78+
</property>
79+
<property name="sizeHint" stdset="0">
80+
<size>
81+
<width>200</width>
82+
<height>20</height>
83+
</size>
84+
</property>
85+
</spacer>
86+
</item>
87+
<item row="2" column="0" colspan="4">
88+
<layout class="QGridLayout" name="gridLayout_3">
89+
<item row="0" column="0">
90+
<widget class="QLabel" name="mXMinLabel">
91+
<property name="text">
92+
<string>X min</string>
93+
</property>
94+
</widget>
95+
</item>
96+
<item row="0" column="1">
97+
<widget class="QDoubleSpinBox" name="mXMinSpinBox">
98+
<property name="minimum">
99+
<double>-999999999.000000000000000</double>
100+
</property>
101+
<property name="maximum">
102+
<double>999999999.000000000000000</double>
103+
</property>
104+
</widget>
105+
</item>
106+
<item row="0" column="2">
107+
<widget class="QLabel" name="mXMaxLabel">
108+
<property name="text">
109+
<string>XMax</string>
110+
</property>
111+
</widget>
112+
</item>
113+
<item row="0" column="3">
114+
<widget class="QDoubleSpinBox" name="mXMaxSpinBox">
115+
<property name="minimum">
116+
<double>-999999999.000000000000000</double>
117+
</property>
118+
<property name="maximum">
119+
<double>999999999.000000000000000</double>
120+
</property>
121+
</widget>
122+
</item>
123+
<item row="1" column="0">
124+
<widget class="QLabel" name="mYMinLabel">
125+
<property name="text">
126+
<string>Y min</string>
127+
</property>
128+
</widget>
129+
</item>
130+
<item row="1" column="1">
131+
<widget class="QDoubleSpinBox" name="mYMinSpinBox">
132+
<property name="minimum">
133+
<double>-999999999.000000000000000</double>
134+
</property>
135+
<property name="maximum">
136+
<double>999999999.000000000000000</double>
137+
</property>
138+
</widget>
139+
</item>
140+
<item row="1" column="2">
141+
<widget class="QLabel" name="mYMaxLabel">
142+
<property name="text">
143+
<string>Y max</string>
144+
</property>
145+
</widget>
146+
</item>
147+
<item row="1" column="3">
148+
<widget class="QDoubleSpinBox" name="mYMaxSpinBox">
149+
<property name="minimum">
150+
<double>-999999999.000000000000000</double>
151+
</property>
152+
<property name="maximum">
153+
<double>999999999.000000000000000</double>
154+
</property>
155+
</widget>
156+
</item>
157+
<item row="2" column="0">
158+
<widget class="QLabel" name="mColumnsLabel">
159+
<property name="text">
160+
<string>Columns</string>
161+
</property>
162+
</widget>
163+
</item>
164+
<item row="2" column="1">
165+
<widget class="QSpinBox" name="mNColumnsSpinBox">
166+
<property name="maximum">
167+
<number>999999999</number>
168+
</property>
169+
</widget>
170+
</item>
171+
<item row="2" column="2">
172+
<widget class="QLabel" name="mRowsLabel">
173+
<property name="text">
174+
<string>Rows</string>
175+
</property>
176+
</widget>
177+
</item>
178+
<item row="2" column="3">
179+
<widget class="QSpinBox" name="mNRowsSpinBox">
180+
<property name="maximum">
181+
<number>999999999</number>
182+
</property>
183+
</widget>
184+
</item>
185+
</layout>
186+
</item>
187+
<item row="3" column="0">
188+
<widget class="QLabel" name="mOutputFormatLabel">
189+
<property name="text">
190+
<string>Output format</string>
191+
</property>
192+
</widget>
193+
</item>
194+
<item row="3" column="2">
195+
<widget class="QComboBox" name="mOutputFormatComboBox"/>
196+
</item>
197+
<item row="4" column="0" colspan="2">
198+
<widget class="QCheckBox" name="mAddResultToProjectCheckBox">
199+
<property name="text">
200+
<string>Add result to project</string>
201+
</property>
202+
</widget>
203+
</item>
204+
</layout>
205+
</widget>
206+
</item>
207+
<item row="1" column="0" colspan="3">
208+
<widget class="QGroupBox" name="mOperatorsGroupBox">
209+
<property name="title">
210+
<string>Operators</string>
211+
</property>
212+
<layout class="QGridLayout" name="gridLayout">
213+
<item row="0" column="0">
214+
<widget class="QPushButton" name="mPlusPushButton">
215+
<property name="text">
216+
<string>+</string>
217+
</property>
218+
</widget>
219+
</item>
220+
<item row="0" column="1">
221+
<widget class="QPushButton" name="mMultiplyPushButton">
222+
<property name="text">
223+
<string>*</string>
224+
</property>
225+
</widget>
226+
</item>
227+
<item row="0" column="2">
228+
<widget class="QPushButton" name="mSqrtButton">
229+
<property name="text">
230+
<string>sqrt</string>
231+
</property>
232+
</widget>
233+
</item>
234+
<item row="0" column="3">
235+
<widget class="QPushButton" name="mSinButton">
236+
<property name="text">
237+
<string>sin</string>
238+
</property>
239+
</widget>
240+
</item>
241+
<item row="0" column="4">
242+
<widget class="QPushButton" name="mExpButton">
243+
<property name="text">
244+
<string>^</string>
245+
</property>
246+
</widget>
247+
</item>
248+
<item row="0" column="5">
249+
<widget class="QPushButton" name="mACosButton">
250+
<property name="text">
251+
<string>acos</string>
252+
</property>
253+
</widget>
254+
</item>
255+
<item row="0" column="6">
256+
<widget class="QPushButton" name="mOpenBracketPushButton">
257+
<property name="text">
258+
<string>(</string>
259+
</property>
260+
</widget>
261+
</item>
262+
<item row="1" column="0">
263+
<widget class="QPushButton" name="mMinusPushButton">
264+
<property name="text">
265+
<string>-</string>
266+
</property>
267+
</widget>
268+
</item>
269+
<item row="1" column="1">
270+
<widget class="QPushButton" name="mDividePushButton">
271+
<property name="text">
272+
<string>/</string>
273+
</property>
274+
</widget>
275+
</item>
276+
<item row="1" column="2">
277+
<widget class="QPushButton" name="mCosButton">
278+
<property name="text">
279+
<string>cos</string>
280+
</property>
281+
</widget>
282+
</item>
283+
<item row="1" column="3">
284+
<widget class="QPushButton" name="mASinButton">
285+
<property name="text">
286+
<string>asin</string>
287+
</property>
288+
</widget>
289+
</item>
290+
<item row="1" column="4">
291+
<widget class="QPushButton" name="mTanButton">
292+
<property name="text">
293+
<string>tan</string>
294+
</property>
295+
</widget>
296+
</item>
297+
<item row="1" column="5">
298+
<widget class="QPushButton" name="mATanButton">
299+
<property name="text">
300+
<string>atan</string>
301+
</property>
302+
</widget>
303+
</item>
304+
<item row="1" column="6">
305+
<widget class="QPushButton" name="mCloseBracketPushButton">
306+
<property name="text">
307+
<string>)</string>
308+
</property>
309+
</widget>
310+
</item>
311+
</layout>
312+
</widget>
313+
</item>
314+
<item row="2" column="0">
315+
<widget class="QLabel" name="mRasterCalculatorExpressionLabel">
316+
<property name="text">
317+
<string>Raster calculator expression</string>
318+
</property>
319+
<property name="buddy">
320+
<cstring>mExpressionTextEdit</cstring>
321+
</property>
322+
</widget>
323+
</item>
324+
<item row="3" column="0" colspan="3">
325+
<widget class="QTextEdit" name="mExpressionTextEdit"/>
326+
</item>
327+
<item row="4" column="1" colspan="2">
328+
<widget class="QDialogButtonBox" name="mButtonBox">
329+
<property name="orientation">
330+
<enum>Qt::Horizontal</enum>
331+
</property>
332+
<property name="standardButtons">
333+
<set>QDialogButtonBox::Cancel|QDialogButtonBox::Ok</set>
334+
</property>
335+
</widget>
336+
</item>
337+
<item row="4" column="0">
338+
<widget class="QLabel" name="mExpressionValidLabel">
339+
<property name="text">
340+
<string/>
341+
</property>
342+
</widget>
343+
</item>
344+
</layout>
345+
</widget>
346+
<resources/>
347+
<connections>
348+
<connection>
349+
<sender>mButtonBox</sender>
350+
<signal>accepted()</signal>
351+
<receiver>QgsRasterCalcDialogBase</receiver>
352+
<slot>accept()</slot>
353+
<hints>
354+
<hint type="sourcelabel">
355+
<x>248</x>
356+
<y>254</y>
357+
</hint>
358+
<hint type="destinationlabel">
359+
<x>157</x>
360+
<y>274</y>
361+
</hint>
362+
</hints>
363+
</connection>
364+
<connection>
365+
<sender>mButtonBox</sender>
366+
<signal>rejected()</signal>
367+
<receiver>QgsRasterCalcDialogBase</receiver>
368+
<slot>reject()</slot>
369+
<hints>
370+
<hint type="sourcelabel">
371+
<x>316</x>
372+
<y>260</y>
373+
</hint>
374+
<hint type="destinationlabel">
375+
<x>286</x>
376+
<y>274</y>
377+
</hint>
378+
</hints>
379+
</connection>
380+
</connections>
381+
</ui>

0 commit comments

Comments
 (0)
Please sign in to comment.