[QGIS Commit] r14503 - in trunk/qgis/src: analysis/raster ui
svn_qgis at osgeo.org
svn_qgis at osgeo.org
Thu Nov 4 07:11:56 EDT 2010
Author: mhugent
Date: 2010-11-04 04:11:56 -0700 (Thu, 04 Nov 2010)
New Revision: 14503
Modified:
trunk/qgis/src/analysis/raster/qgsrastercalcnode.cpp
trunk/qgis/src/analysis/raster/qgsrastercalculator.cpp
trunk/qgis/src/analysis/raster/qgsrastermatrix.cpp
trunk/qgis/src/analysis/raster/qgsrastermatrix.h
trunk/qgis/src/ui/qgsrastercalcdialogbase.ui
Log:
Fix consideration of nodata values in raster calculator, code cleanups
Modified: trunk/qgis/src/analysis/raster/qgsrastercalcnode.cpp
===================================================================
--- trunk/qgis/src/analysis/raster/qgsrastercalcnode.cpp 2010-11-04 07:36:56 UTC (rev 14502)
+++ trunk/qgis/src/analysis/raster/qgsrastercalcnode.cpp 2010-11-04 11:11:56 UTC (rev 14503)
@@ -1,4 +1,5 @@
#include "qgsrastercalcnode.h"
+#include <cfloat>
QgsRasterCalcNode::QgsRasterCalcNode(): mLeft( 0 ), mRight( 0 ), mRasterMatrix( 0 ), mNumber( 0 )
{
@@ -18,11 +19,11 @@
QgsRasterCalcNode::~QgsRasterCalcNode()
{
- if ( mLeft )
+ if( mLeft )
{
delete mLeft;
}
- if ( mRight )
+ if( mRight )
{
delete mRight;
}
@@ -33,10 +34,10 @@
//if type is raster ref: return a copy of the corresponding matrix
//if type is operator, call the proper matrix operations
- if ( mType == tRasterRef )
+ if( mType == tRasterRef )
{
QMap<QString, QgsRasterMatrix*>::iterator it = rasterData.find( mRasterName );
- if ( it == rasterData.end() )
+ if( it == rasterData.end() )
{
return false;
}
@@ -44,23 +45,23 @@
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 );
+ result.setData(( *it )->nColumns(), ( *it )->nRows(), data, ( *it )->nodataValue() );
return true;
}
- else if ( mType == tOperator )
+ else if( mType == tOperator )
{
QgsRasterMatrix leftMatrix, rightMatrix;
QgsRasterMatrix resultMatrix;
- if ( !mLeft || !mLeft->calculate( rasterData, leftMatrix ) )
+ if( !mLeft || !mLeft->calculate( rasterData, leftMatrix ) )
{
return false;
}
- if ( mRight && !mRight->calculate( rasterData, rightMatrix ) )
+ if( mRight && !mRight->calculate( rasterData, rightMatrix ) )
{
return false;
}
- switch ( mOperator )
+ switch( mOperator )
{
case opPLUS:
leftMatrix.add( rightMatrix );
@@ -121,14 +122,14 @@
}
int newNColumns = leftMatrix.nColumns();
int newNRows = leftMatrix.nRows();
- result.setData( newNColumns, newNRows, leftMatrix.takeData() );
+ result.setData( newNColumns, newNRows, leftMatrix.takeData(), leftMatrix.nodataValue() );
return true;
}
- else if ( mType == tNumber )
+ else if( mType == tNumber )
{
float* data = new float[1];
data[0] = mNumber;
- result.setData( 1, 1, data );
+ result.setData( 1, 1, data, -FLT_MAX );
return true;
}
return false;
Modified: trunk/qgis/src/analysis/raster/qgsrastercalculator.cpp
===================================================================
--- trunk/qgis/src/analysis/raster/qgsrastercalculator.cpp 2010-11-04 07:36:56 UTC (rev 14502)
+++ trunk/qgis/src/analysis/raster/qgsrastercalculator.cpp 2010-11-04 11:11:56 UTC (rev 14503)
@@ -24,7 +24,7 @@
QgsRasterCalculator::QgsRasterCalculator( const QString& formulaString, const QString& outputFile, const QString& outputFormat,
const QgsRectangle& outputExtent, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry>& rasterEntries ): mFormulaString( formulaString ), mOutputFile( outputFile ), mOutputFormat( outputFormat ),
- mOutputRectangle( outputExtent ), mNumOutputColumns( nOutputColumns ), mNumOutputRows( nOutputRows ), mRasterEntries( rasterEntries )
+ mOutputRectangle( outputExtent ), mNumOutputColumns( nOutputColumns ), mNumOutputRows( nOutputRows ), mRasterEntries( rasterEntries )
{
}
@@ -37,7 +37,7 @@
//prepare search string / tree
QString errorString;
QgsRasterCalcNode* calcNode = QgsRasterCalcNode::parseRasterCalcString( mFormulaString, errorString );
- if ( !calcNode )
+ if( !calcNode )
{
//error
}
@@ -51,39 +51,47 @@
QVector< GDALDatasetH > mInputDatasets; //raster references and corresponding dataset
QVector<QgsRasterCalculatorEntry>::const_iterator it = mRasterEntries.constBegin();
- for ( ; it != mRasterEntries.constEnd(); ++it )
+ for( ; it != mRasterEntries.constEnd(); ++it )
{
- if ( !it->raster ) // no raster layer in entry
+ if( !it->raster ) // no raster layer in entry
{
return 2;
}
GDALDatasetH inputDataset = GDALOpen( it->raster->source().toLocal8Bit().data(), GA_ReadOnly );
- if ( inputDataset == NULL )
+ if( inputDataset == NULL )
{
return 2;
}
GDALRasterBandH inputRasterBand = GDALGetRasterBand( inputDataset, it->bandNumber );
- if ( inputRasterBand == NULL )
+ if( inputRasterBand == NULL )
{
return 2;
}
+ int nodataSuccess;
+ double nodataValue = GDALGetRasterNoDataValue( inputRasterBand, &nodataSuccess );
+
mInputDatasets.push_back( inputDataset );
mInputRasterBands.insert( it->ref, inputRasterBand );
- inputScanLineData.insert( it->ref, new QgsRasterMatrix( mNumOutputColumns, 1, new float[mNumOutputColumns] ) );
+ inputScanLineData.insert( it->ref, new QgsRasterMatrix( mNumOutputColumns, 1, new float[mNumOutputColumns], nodataValue ) );
}
//open output dataset for writing
GDALDriverH outputDriver = openOutputDriver();
- if ( outputDriver == NULL )
+ if( outputDriver == NULL )
{
return 1;
}
GDALDatasetH outputDataset = openOutputFile( outputDriver );
GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset, 1 );
+ int nodataSuccess;
+
+ float outputNodataValue = -FLT_MAX;
+ GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
+
float* resultScanLine = ( float * ) CPLMalloc( sizeof( float ) * mNumOutputColumns );
- if ( p )
+ if( p )
{
p->setMaximum( mNumOutputRows );
}
@@ -91,21 +99,21 @@
QgsRasterMatrix resultMatrix;
//read / write line by line
- for ( int i = 0; i < mNumOutputRows; ++i )
+ for( int i = 0; i < mNumOutputRows; ++i )
{
- if ( p )
+ if( p )
{
p->setValue( i );
}
- if ( p && p->wasCanceled() )
+ if( p && p->wasCanceled() )
{
break;
}
//fill buffers
QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin();
- for ( ; bufferIt != inputScanLineData.end(); ++bufferIt )
+ for( ; bufferIt != inputScanLineData.end(); ++bufferIt )
{
double sourceTransformation[6];
GDALRasterBandH sourceRasterBand = mInputRasterBands[bufferIt.key()];
@@ -114,18 +122,48 @@
readRasterPart( targetGeoTransform, 0, i, mNumOutputColumns, 1, sourceTransformation, sourceRasterBand, bufferIt.value()->data() );
}
- if ( calcNode->calculate( inputScanLineData, resultMatrix ) )
+ if( calcNode->calculate( inputScanLineData, resultMatrix ) )
{
+ bool resultIsNumber = resultMatrix.isNumber();
+ float* calcData;
+
+ if( resultIsNumber ) //scalar result. Insert number for every pixel
+ {
+ calcData = new float[mNumOutputColumns];
+ for( int j = 0; j < mNumOutputColumns; ++j )
+ {
+ calcData[j] = resultMatrix.number();
+ }
+ }
+ else //result is real matrix
+ {
+ calcData = resultMatrix.data();
+ }
+
+ //replace all matrix nodata values with output nodatas
+ for( int j = 0; j < mNumOutputColumns; ++j )
+ {
+ if( calcData[j] == resultMatrix.nodataValue() )
+ {
+ calcData[j] = outputNodataValue;
+ }
+ }
+
//write scanline to the dataset
- if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, resultMatrix.data(), mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
+ if( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
{
qWarning( "RasterIO error!" );
}
+
+ if( resultIsNumber )
+ {
+ delete[] calcData;
+ }
}
}
- if ( p )
+ if( p )
{
p->setValue( mNumOutputRows );
}
@@ -133,19 +171,19 @@
//close datasets and release memory
delete calcNode;
QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin();
- for ( ; bufferIt != inputScanLineData.end(); ++bufferIt )
+ for( ; bufferIt != inputScanLineData.end(); ++bufferIt )
{
delete bufferIt.value();
}
inputScanLineData.clear();
QVector< GDALDatasetH >::iterator datasetIt = mInputDatasets.begin();
- for ( ; datasetIt != mInputDatasets.end(); ++ datasetIt )
+ for( ; datasetIt != mInputDatasets.end(); ++ datasetIt )
{
GDALClose( *datasetIt );
}
- if ( p && p->wasCanceled() )
+ if( p && p->wasCanceled() )
{
//delete the dataset without closing (because it is faster)
GDALDeleteDataset( outputDriver, mOutputFile.toLocal8Bit().data() );
@@ -167,13 +205,13 @@
//open driver
GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
- if ( outputDriver == NULL )
+ if( outputDriver == NULL )
{
return outputDriver; //return NULL, driver does not exist
}
driverMetadata = GDALGetMetadata( outputDriver, NULL );
- if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
+ if( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
{
return NULL; //driver exist, but it does not support the create operation
}
@@ -186,7 +224,7 @@
//open output file
char **papszOptions = NULL;
GDALDatasetH outputDataset = GDALCreate( outputDriver, mOutputFile.toLocal8Bit().data(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions );
- if ( outputDataset == NULL )
+ if( outputDataset == NULL )
{
return outputDataset;
}
@@ -202,7 +240,7 @@
void QgsRasterCalculator::readRasterPart( double* targetGeotransform, int xOffset, int yOffset, int nCols, int nRows, double* sourceTransform, GDALRasterBandH sourceBand, float* rasterBuffer )
{
//If dataset transform is the same as the requested transform, do a normal GDAL raster io
- if ( transformationsEqual( targetGeotransform, sourceTransform ) )
+ if( transformationsEqual( targetGeotransform, sourceTransform ) )
{
GDALRasterIO( sourceBand, GF_Read, xOffset, yOffset, nCols, nRows, rasterBuffer, nCols, nRows, GDT_Float32, 0, 0 );
return;
@@ -213,14 +251,15 @@
double nodataValue = GDALGetRasterNoDataValue( sourceBand, &nodataSuccess );
QgsRectangle targetRect( targetGeotransform[0] + targetGeotransform[1] * xOffset, targetGeotransform[3] + yOffset * targetGeotransform[5] + nRows * targetGeotransform[5]
, targetGeotransform[0] + targetGeotransform[1] * xOffset + targetGeotransform[1] * nCols, targetGeotransform[3] + yOffset * targetGeotransform[5] );
- QgsRectangle sourceRect( sourceTransform[0], sourceTransform[3] + GDALGetRasterBandXSize( sourceBand ) * sourceTransform[5], sourceTransform[0] + GDALGetRasterBandXSize( sourceBand )* sourceTransform[1], sourceTransform[3] );
+ QgsRectangle sourceRect( sourceTransform[0], sourceTransform[3] + GDALGetRasterBandYSize( sourceBand ) * sourceTransform[5],
+ sourceTransform[0] + GDALGetRasterBandXSize( sourceBand )* sourceTransform[1], sourceTransform[3] );
QgsRectangle intersection = targetRect.intersect( &sourceRect );
//no intersection, fill all the pixels with nodata values
- if ( intersection.isEmpty() )
+ if( intersection.isEmpty() )
{
int nPixels = nCols * nRows;
- for ( int i = 0; i < nPixels; ++i )
+ for( int i = 0; i < nPixels; ++i )
{
rasterBuffer[i] = nodataValue;
}
@@ -246,17 +285,17 @@
double targetPixelY = targetGeotransform[3] + targetGeotransform[5] * yOffset + targetGeotransform[5] / 2.0; //coordinates of current target pixel
int sourceIndexX, sourceIndexY; //current raster index in source pixels
double sx, sy;
- for ( int i = 0; i < nRows; ++i )
+ for( int i = 0; i < nRows; ++i )
{
targetPixelX = targetPixelXMin;
- for ( int j = 0; j < nCols; ++j )
+ for( int j = 0; j < nCols; ++j )
{
sx = ( targetPixelX - sourceRasterXMin ) / sourceTransform[1];
sourceIndexX = sx > 0 ? sx : floor( sx );
sy = ( targetPixelY - sourceRasterYMax ) / sourceTransform[5];
sourceIndexY = sy > 0 ? sy : floor( sy );
- if ( sourceIndexX >= 0 && sourceIndexX < nSourcePixelsX
- && sourceIndexY >= 0 && sourceIndexY < nSourcePixelsY )
+ if( sourceIndexX >= 0 && sourceIndexX < nSourcePixelsX
+ && sourceIndexY >= 0 && sourceIndexY < nSourcePixelsY )
{
rasterBuffer[j + i*nRows] = sourceRaster[ sourceIndexX + nSourcePixelsX * sourceIndexY ];
}
@@ -271,86 +310,13 @@
CPLFree( sourceRaster );
return;
-
-#if 0
- //If dataset transform is the same as the requested transform, do a normal GDAL raster io
- if ( transformationsEqual( targetGeotransform, sourceTransform ) )
- {
- GDALRasterIO( sourceBand, GF_Read, xOffset, yOffset, nCols, nRows, rasterBuffer, nCols, nRows, GDT_Float32, 0, 0 );
- return;
- }
-
- //pixel calculation needed because of different raster position / resolution
-
- //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)
- //calculate offset
-
- //QgsRectangle targetRect( targetGeotransform[0], targetGeotransform[3] + mNumOutputColumns * targetGeotransform[5], targetGeotransform[0] + mNumOutputRows * targetGeotransform[1], targetGeotransform[3]);
- QgsRectangle targetRect( targetGeotransform[0] + targetGeotransform[1] * xOffset, targetGeotransform[3] + yOffset * targetGeotransform[5] + nRows * targetGeotransform[5]
- , targetGeotransform[0] + targetGeotransform[1] * xOffset + targetGeotransform[1] * nCols, targetGeotransform[3] + yOffset * targetGeotransform[5] );
- QgsRectangle sourceRect( sourceTransform[0], sourceTransform[3] + GDALGetRasterBandXSize( sourceBand ) * sourceTransform[5], sourceTransform[0] + GDALGetRasterBandYSize( sourceBand )* sourceTransform[1], sourceTransform[3] );
- QgsRectangle intersection = targetRect.intersect( &sourceRect );
-
- //no intersection, fill all the pixels with nodata values
- if ( intersection.isEmpty() )
- {
- int nPixels = nCols * nRows;
- for ( int i = 0; i < nPixels; ++i )
- {
- rasterBuffer[i] = 0;
- }
- return;
- }
-
- //do raster io in source resolution
- int sourcePixelOffsetX = ( intersection.xMinimum() - sourceTransform[0] ) / sourceTransform[1];
- int sourcePixelOffsetY = ( intersection.yMaximum() - sourceTransform[3] ) / sourceTransform[5];
- int nSourcePixelsX = ceil( fabs( intersection.width() / sourceTransform[1] ) );
- int nSourcePixelsY = ceil( fabs( intersection.height() / sourceTransform[5] ) );
- float* sourceRaster = ( float * ) CPLMalloc( sizeof( float ) * nSourcePixelsX * nSourcePixelsY );
-
- double sourceRasterXMin = sourceRect.xMinimum() + sourcePixelOffsetX * sourceTransform[1];
- double sourceRasterYMax = sourceRect.yMaximum() + sourcePixelOffsetY * sourceTransform[5];
-
- GDALRasterIO( sourceBand, GF_Read, sourcePixelOffsetX, sourcePixelOffsetY, nSourcePixelsX, nSourcePixelsY,
- sourceRaster, nSourcePixelsX, nSourcePixelsY, GDT_Float32, 0, 0 );
-
-
- double targetPixelX;
- double targetPixelXMin = targetGeotransform[0] + targetGeotransform[1] * xOffset + targetGeotransform[1] / 2.0;
- double targetPixelY = targetGeotransform[3] + targetGeotransform[5] * yOffset + targetGeotransform[5] / 2.0; //coordinates of current target pixel
- int sourceIndexX, sourceIndexY; //current raster index in source pixels
- for ( int i = 0; i < nRows; ++i )
- {
- targetPixelX = targetPixelXMin;
- for ( int j = 0; j < nCols; ++j )
- {
- sourceIndexX = ( targetPixelX - sourceRasterXMin ) / sourceTransform[1];
- sourceIndexY = ( targetPixelY - sourceRasterYMax ) / sourceTransform[5];
- if ( sourceIndexX >= 0 && sourceIndexX < nSourcePixelsX
- && sourceIndexY >= 0 && sourceIndexY < nSourcePixelsY )
- {
- rasterBuffer[j + i*nCols] = sourceRaster[ sourceIndexX + nSourcePixelsX * sourceIndexY ];
- }
- else
- {
- rasterBuffer[j + i*nCols] = 0; //todo: insert null value of rasterband here
- }
- targetPixelX += targetGeotransform[1];
- }
- targetPixelY += targetGeotransform[5];
- }
-
- CPLFree( sourceRaster );
- return;
-#endif //0
}
bool QgsRasterCalculator::transformationsEqual( double* t1, double* t2 ) const
{
- for ( int i = 0; i < 6; ++i )
+ for( int i = 0; i < 6; ++i )
{
- if ( !doubleNear( t1[i], t2[i] ) )
+ if( !doubleNear( t1[i], t2[i], 0.00001 ) )
{
return false;
}
Modified: trunk/qgis/src/analysis/raster/qgsrastermatrix.cpp
===================================================================
--- trunk/qgis/src/analysis/raster/qgsrastermatrix.cpp 2010-11-04 07:36:56 UTC (rev 14502)
+++ trunk/qgis/src/analysis/raster/qgsrastermatrix.cpp 2010-11-04 11:11:56 UTC (rev 14503)
@@ -28,7 +28,8 @@
{
}
-QgsRasterMatrix::QgsRasterMatrix( int nCols, int nRows, float* data ): mColumns( nCols ), mRows( nRows ), mData( data )
+QgsRasterMatrix::QgsRasterMatrix( int nCols, int nRows, float* data, double nodataValue ):
+ mColumns( nCols ), mRows( nRows ), mData( data ), mNodataValue( nodataValue )
{
}
@@ -50,15 +51,17 @@
int nEntries = mColumns * mRows;
mData = new float[nEntries];
memcpy( mData, m.mData, sizeof( float ) * nEntries );
+ mNodataValue = m.nodataValue();
return *this;
}
-void QgsRasterMatrix::setData( int cols, int rows, float* data )
+void QgsRasterMatrix::setData( int cols, int rows, float* data, double nodataValue )
{
delete[] mData;
mColumns = cols;
mRows = rows;
mData = data;
+ mNodataValue = nodataValue;
}
float* QgsRasterMatrix::takeData()
@@ -125,123 +128,101 @@
bool QgsRasterMatrix::squareRoot()
{
- if ( !mData )
- {
- return false;
- }
-
- int nEntries = mColumns * mRows;
- for ( int i = 0; i < nEntries; ++i )
- {
- double value = mData[i];
- if ( value >= 0 )
- {
- mData[i] = sqrt( value );
- }
- else
- {
- mData[i] = -10000; //no complex numbers
- }
- }
- return true;
+ return oneArgumentOperation( opSQRT );
}
bool QgsRasterMatrix::sinus()
{
- if ( !mData )
- {
- return false;
- }
-
- int nEntries = mColumns * mRows;
- for ( int i = 0; i < nEntries; ++i )
- {
- mData[i] = sin( mData[i] );
- }
- return true;
+ return oneArgumentOperation( opSIN );
}
bool QgsRasterMatrix::asinus()
{
- if ( !mData )
- {
- return false;
- }
-
- int nEntries = mColumns * mRows;
- for ( int i = 0; i < nEntries; ++i )
- {
- mData[i] = asin( mData[i] );
- }
- return true;
+ return oneArgumentOperation( opASIN );
}
bool QgsRasterMatrix::cosinus()
{
- if ( !mData )
- {
- return false;
- }
-
- int nEntries = mColumns * mRows;
- for ( int i = 0; i < nEntries; ++i )
- {
- mData[i] = cos( mData[i] );
- }
- return true;
+ return oneArgumentOperation( opCOS );
}
bool QgsRasterMatrix::acosinus()
{
- if ( !mData )
- {
- return false;
- }
-
- int nEntries = mColumns * mRows;
- for ( int i = 0; i < nEntries; ++i )
- {
- mData[i] = acos( mData[i] );
- }
- return true;
+ return oneArgumentOperation( opACOS );
}
bool QgsRasterMatrix::tangens()
{
- if ( !mData )
- {
- return false;
- }
-
- int nEntries = mColumns * mRows;
- for ( int i = 0; i < nEntries; ++i )
- {
- mData[i] = tan( mData[i] );
- }
- return true;
+ return oneArgumentOperation( opTAN );
}
bool QgsRasterMatrix::atangens()
{
- if ( !mData )
+ return oneArgumentOperation( opATAN );
+}
+
+bool QgsRasterMatrix::oneArgumentOperation( OneArgOperator op )
+{
+ if( !mData )
{
return false;
}
int nEntries = mColumns * mRows;
- for ( int i = 0; i < nEntries; ++i )
+ double value;
+ for( int i = 0; i < nEntries; ++i )
{
- mData[i] = atan( mData[i] );
+ value = mData[i];
+ if( value != mNodataValue )
+ {
+ switch( op )
+ {
+ case opSQRT:
+ if( value < 0 ) //no complex numbers
+ {
+ mData[i] = mNodataValue;
+ }
+ else
+ {
+ mData[i] = sqrt( value );
+ }
+ break;
+ case opSIN:
+ mData[i] = sin( value );
+ break;
+ case opCOS:
+ mData[i] = cos( value );
+ break;
+ case opTAN:
+ mData[i] = tan( value );
+ break;
+ case opASIN:
+ mData[i] = asin( value );
+ break;
+ case opACOS:
+ mData[i] = acos( value );
+ break;
+ case opATAN:
+ mData[i] = atan( value );
+ break;
+ }
+ }
}
return true;
}
bool QgsRasterMatrix::twoArgumentOperation( TwoArgOperator op, const QgsRasterMatrix& other )
{
- if ( isNumber() && other.isNumber() )
+ if( isNumber() && other.isNumber() ) //operation on two 1x1 matrices
{
- switch ( op )
+ //operations with nodata values always generate nodata
+ if( mData[0] == mNodataValue || other.number() == other.nodataValue() )
{
+ mData[0] = mNodataValue;
+ return true;
+ }
+ switch( op )
+ {
case opPLUS:
mData[0] = number() + other.number();
break;
@@ -252,9 +233,9 @@
mData[0] = number() * other.number();
break;
case opDIV:
- if ( other.number() == 0 )
+ if( other.number() == 0 )
{
- mData[0] = -10000;
+ mData[0] = mNodataValue;
}
else
{
@@ -262,9 +243,9 @@
}
break;
case opPOW:
- if ( !testPowerValidity( mData[0], ( float ) other.number() ) )
+ if( !testPowerValidity( mData[0], ( float ) other.number() ) )
{
- mData[0] = -10000;
+ mData[0] = mNodataValue;
}
else
{
@@ -292,286 +273,237 @@
}
return true;
}
+
//two matrices
- if ( !isNumber() && !other.isNumber() )
+ if( !isNumber() && !other.isNumber() )
{
float* matrix = other.mData;
int nEntries = mColumns * mRows;
- switch ( op )
+ double value1, value2;
+
+ for( int i = 0; i < nEntries; ++i )
{
- case opPLUS:
- for ( int i = 0; i < nEntries; ++i )
+ value1 = mData[i]; value2 = matrix[i];
+ if( value1 == mNodataValue || value2 == other.mNodataValue )
+ {
+ mData[i] = mNodataValue;
+ }
+ else
+ {
+ switch( op )
{
- mData[i] = mData[i] + matrix[i];
+ case opPLUS:
+ mData[i] = value1 + value2;
+ break;
+ case opMINUS:
+ mData[i] = value1 - value2;
+ break;
+ case opMUL:
+ mData[i] = value1 * value2;
+ break;
+ case opDIV:
+ if( value2 == 0 )
+ {
+ mData[i] = mNodataValue;
+ }
+ else
+ {
+ mData[i] = value1 / value2;
+ }
+ break;
+ case opPOW:
+ if( !testPowerValidity( value1, value2 ) )
+ {
+ mData[i] = mNodataValue;
+ }
+ else
+ {
+ mData[i] = pow( value1, value2 );
+ }
+ break;
+ case opEQ:
+ mData[i] = ( value1 == value2 ? 1 : 0 );
+ break;
+ case opNE:
+ mData[i] = ( value1 == value2 ? 0 : 1 );
+ break;
+ case opGT:
+ mData[i] = ( value1 > value2 ? 1 : 0 );
+ break;
+ case opLT:
+ mData[i] = ( value1 < value2 ? 1 : 0 );
+ break;
+ case opGE:
+ mData[i] = ( value1 >= value2 ? 1 : 0 );
+ break;
+ case opLE:
+ mData[i] = ( value1 <= value2 ? 1 : 0 );
+ break;
}
- break;
- case opMINUS:
- for ( int i = 0; i < nEntries; ++i )
- {
- mData[i] = mData[i] - matrix[i];
- }
- break;
- case opMUL:
- for ( int i = 0; i < nEntries; ++i )
- {
- mData[i] = mData[i] * matrix[i];
- }
- break;
- case opDIV:
- for ( int i = 0; i < nEntries; ++i )
- {
- if ( matrix[i] == 0 )
- {
- mData[i] = -10000;
- }
- else
- {
- mData[i] = mData[i] / matrix[i];
- }
- }
- break;
- case opPOW:
- for ( int i = 0; i < nEntries; ++i )
- {
- if ( !testPowerValidity( mData[i], matrix[i] ) )
- {
- mData[i] = -10000;
- }
- else
- {
- mData[i] = pow( mData[i], matrix[i] );
- }
- }
- break;
- case opEQ:
- for ( int i = 0; i < nEntries; ++i )
- {
- mData[i] = ( mData[i] == matrix[i] ? 1 : 0 );
- }
- break;
- case opNE:
- for ( int i = 0; i < nEntries; ++i )
- {
- mData[i] = ( mData[i] == matrix[i] ? 0 : 1 );
- }
- break;
- case opGT:
- for ( int i = 0; i < nEntries; ++i )
- {
- mData[i] = ( mData[i] > matrix[i] ? 1 : 0 );
- }
- break;
- case opLT:
- for ( int i = 0; i < nEntries; ++i )
- {
- mData[i] = ( mData[i] < matrix[i] ? 1 : 0 );
- }
- break;
- case opGE:
- for ( int i = 0; i < nEntries; ++i )
- {
- mData[i] = ( mData[i] >= matrix[i] ? 1 : 0 );
- }
- break;
- case opLE:
- for ( int i = 0; i < nEntries; ++i )
- {
- mData[i] = ( mData[i] <= matrix[i] ? 1 : 0 );
- }
- break;
+ }
}
return true;
}
- double value = 0;
- if ( isNumber() )
+
+ //this matrix is a single number and the other one a real matrix
+ if( isNumber() )
{
float* matrix = other.mData;
int nEntries = other.nColumns() * other.nRows();
- value = mData[0];
+ double value = mData[0];
delete[] mData;
mData = new float[nEntries]; mColumns = other.nColumns(); mRows = other.nRows();
+ mNodataValue = other.nodataValue();
- switch ( op )
+ if( value == mNodataValue )
{
- case opPLUS:
- for ( int i = 0; i < nEntries; ++i )
- {
+ for( int i = 0; i < nEntries; ++i )
+ {
+ mData[i] = mNodataValue;
+ }
+ return true;
+ }
+
+ for( int i = 0; i < nEntries; ++i )
+ {
+ if( matrix[i] == other.mNodataValue )
+ {
+ mData[i] = mNodataValue;
+ continue;
+ }
+
+ switch( op )
+ {
+ case opPLUS:
mData[i] = value + matrix[i];
- }
- break;
- case opMINUS:
- for ( int i = 0; i < nEntries; ++i )
- {
+ break;
+ case opMINUS:
mData[i] = value - matrix[i];
- }
- break;
- case opMUL:
- for ( int i = 0; i < nEntries; ++i )
- {
+ break;
+ case opMUL:
mData[i] = value * matrix[i];
- }
- break;
- case opDIV:
- for ( int i = 0; i < nEntries; ++i )
- {
- if ( matrix[i] == 0 )
+ break;
+ case opDIV:
+ if( matrix[i] == 0 )
{
- mData[i] = -10000;
+ mData[i] = mNodataValue;
}
else
{
mData[i] = value / matrix[i];
}
- }
- break;
- case opPOW:
- for ( int i = 0; i < nEntries; ++i )
- {
- if ( !testPowerValidity( value, matrix[i] ) )
+ break;
+ case opPOW:
+ if( !testPowerValidity( value, matrix[i] ) )
{
- mData[i] = -10000;
+ mData[i] = mNodataValue;
}
else
{
mData[i] = pow(( float ) value, matrix[i] );
}
- }
- break;
- case opEQ:
- for ( int i = 0; i < nEntries; ++i )
- {
+ break;
+ case opEQ:
mData[i] = ( value == matrix[i] ? 1 : 0 );
- }
- break;
- case opNE:
- for ( int i = 0; i < nEntries; ++i )
- {
+ break;
+ case opNE:
mData[i] = ( value == matrix[i] ? 0 : 1 );
- }
- break;
- case opGT:
- for ( int i = 0; i < nEntries; ++i )
- {
+ break;
+ case opGT:
mData[i] = ( value > matrix[i] ? 1 : 0 );
- }
- break;
- case opLT:
- for ( int i = 0; i < nEntries; ++i )
- {
+ break;
+ case opLT:
mData[i] = ( value < matrix[i] ? 1 : 0 );
- }
- break;
- case opGE:
- for ( int i = 0; i < nEntries; ++i )
- {
+ break;
+ case opGE:
mData[i] = ( value >= matrix[i] ? 1 : 0 );
- }
- break;
- case opLE:
- for ( int i = 0; i < nEntries; ++i )
- {
+ break;
+ case opLE:
mData[i] = ( value <= matrix[i] ? 1 : 0 );
- }
- break;
+ break;
+ }
}
+ return true;
}
else //this matrix is a real matrix and the other a number
{
- value = other.number();
+ double value = other.number();
int nEntries = mColumns * mRows;
- switch ( op )
+
+ if( other.number() == other.mNodataValue )
{
- case opPLUS:
- for ( int i = 0; i < nEntries; ++i )
- {
+ for( int i = 0; i < nEntries; ++i )
+ {
+ mData[i] = mNodataValue;
+ }
+ return true;
+ }
+
+ for( int i = 0; i < nEntries; ++i )
+ {
+ if( mData[i] == mNodataValue )
+ {
+ continue;
+ }
+
+ switch( op )
+ {
+ case opPLUS:
mData[i] = mData[i] + value;
- }
- break;
- case opMINUS:
- for ( int i = 0; i < nEntries; ++i )
- {
+ break;
+ case opMINUS:
mData[i] = mData[i] - value;
- }
- break;
- case opMUL:
- for ( int i = 0; i < nEntries; ++i )
- {
+ break;
+ case opMUL:
mData[i] = mData[i] * value;
- }
- break;
- case opDIV:
- if ( value == 0 )
- {
- for ( int i = 0; i < nEntries; ++i )
+ break;
+ case opDIV:
+ if( value == 0 )
{
- mData[i] = -10000;
+ mData[i] = mNodataValue;
}
- }
- else
- {
- for ( int i = 0; i < nEntries; ++i )
+ else
{
mData[i] = mData[i] / value;
}
- }
- break;
- case opPOW:
- for ( int i = 0; i < nEntries; ++i )
- {
- if ( !testPowerValidity( mData[i], value ) )
+ break;
+ case opPOW:
+ if( !testPowerValidity( mData[i], value ) )
{
- mData[i] = -10000;
+ mData[i] = mNodataValue;
}
else
{
mData[i] = pow( mData[i], ( float ) value );
}
- }
- break;
- case opEQ:
- for ( int i = 0; i < nEntries; ++i )
- {
+ break;
+ case opEQ:
mData[i] = ( mData[i] == value ? 1 : 0 );
- }
- break;
- case opNE:
- for ( int i = 0; i < nEntries; ++i )
- {
+ break;
+ case opNE:
mData[i] = ( mData[i] == value ? 0 : 1 );
- }
- break;
- case opGT:
- for ( int i = 0; i < nEntries; ++i )
- {
+ break;
+ case opGT:
mData[i] = ( mData[i] > value ? 1 : 0 );
- }
- break;
- case opLT:
- for ( int i = 0; i < nEntries; ++i )
- {
+ break;
+ case opLT:
mData[i] = ( mData[i] < value ? 1 : 0 );
- }
- break;
- case opGE:
- for ( int i = 0; i < nEntries; ++i )
- {
+ break;
+ case opGE:
mData[i] = ( mData[i] >= value ? 1 : 0 );
- }
- break;
- case opLE:
- for ( int i = 0; i < nEntries; ++i )
- {
+ break;
+ case opLE:
mData[i] = ( mData[i] <= value ? 1 : 0 );
- }
- break;
+ break;
+ }
}
+ return true;
}
- return true;
}
bool QgsRasterMatrix::testPowerValidity( double base, double power )
{
- if (( base == 0 && power < 0 ) || ( power < 0 && ( power - floor( power ) ) > 0 ) )
+ if(( base == 0 && power < 0 ) || ( power < 0 && ( power - floor( power ) ) > 0 ) )
{
return false;
}
Modified: trunk/qgis/src/analysis/raster/qgsrastermatrix.h
===================================================================
--- trunk/qgis/src/analysis/raster/qgsrastermatrix.h 2010-11-04 07:36:56 UTC (rev 14502)
+++ trunk/qgis/src/analysis/raster/qgsrastermatrix.h 2010-11-04 11:11:56 UTC (rev 14503)
@@ -50,7 +50,7 @@
/**Takes ownership of data array*/
QgsRasterMatrix();
- QgsRasterMatrix( int nCols, int nRows, float* data );
+ QgsRasterMatrix( int nCols, int nRows, float* data, double nodataValue );
QgsRasterMatrix( const QgsRasterMatrix& m );
~QgsRasterMatrix();
@@ -63,11 +63,14 @@
/**Returns data and ownership. Sets data and nrows, ncols of this matrix to 0*/
float* takeData();
- void setData( int cols, int rows, float* data );
+ void setData( int cols, int rows, float* data, double nodataValue );
int nColumns() const { return mColumns; }
int nRows() const { return mRows; }
+ double nodataValue() const { return mNodataValue; }
+ void setNodataValue( double d ) { mNodataValue = d; }
+
QgsRasterMatrix& operator=( const QgsRasterMatrix& m );
/**Adds another matrix to this one*/
bool add( const QgsRasterMatrix& other );
@@ -95,9 +98,12 @@
int mColumns;
int mRows;
float* mData;
+ double mNodataValue;
- /**+,-,*,/,^*/
+ /**+,-,*,/,^,<,>,<=,>=,=,!=*/
bool twoArgumentOperation( TwoArgOperator op, const QgsRasterMatrix& other );
+ /*sqrt, sin, cos, tan, asin, acos, atan*/
+ bool oneArgumentOperation( OneArgOperator op );
bool testPowerValidity( double base, double power );
};
Modified: trunk/qgis/src/ui/qgsrastercalcdialogbase.ui
===================================================================
--- trunk/qgis/src/ui/qgsrastercalcdialogbase.ui 2010-11-04 07:36:56 UTC (rev 14502)
+++ trunk/qgis/src/ui/qgsrastercalcdialogbase.ui 2010-11-04 11:11:56 UTC (rev 14503)
@@ -95,6 +95,9 @@
</item>
<item row="0" column="1">
<widget class="QDoubleSpinBox" name="mXMinSpinBox">
+ <property name="decimals">
+ <number>5</number>
+ </property>
<property name="minimum">
<double>-999999999.000000000000000</double>
</property>
@@ -112,6 +115,9 @@
</item>
<item row="0" column="3">
<widget class="QDoubleSpinBox" name="mXMaxSpinBox">
+ <property name="decimals">
+ <number>5</number>
+ </property>
<property name="minimum">
<double>-999999999.000000000000000</double>
</property>
@@ -129,6 +135,9 @@
</item>
<item row="1" column="1">
<widget class="QDoubleSpinBox" name="mYMinSpinBox">
+ <property name="decimals">
+ <number>5</number>
+ </property>
<property name="minimum">
<double>-999999999.000000000000000</double>
</property>
@@ -146,6 +155,9 @@
</item>
<item row="1" column="3">
<widget class="QDoubleSpinBox" name="mYMaxSpinBox">
+ <property name="decimals">
+ <number>5</number>
+ </property>
<property name="minimum">
<double>-999999999.000000000000000</double>
</property>
More information about the QGIS-commit
mailing list