23 #include <QProgressDialog>
26 #include <cpl_string.h>
27 #include <gdalwarper.h>
29 #if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800
30 #define TO8(x) (x).toUtf8().constData()
31 #define TO8F(x) (x).toUtf8().constData()
33 #define TO8(x) (x).toLocal8Bit().constData()
34 #define TO8F(x) QFile::encodeName( x ).constData()
39 : mFormulaString( formulaString )
40 , mOutputFile( outputFile )
41 , mOutputFormat( outputFormat )
42 , mOutputRectangle( outputExtent )
43 , mNumOutputColumns( nOutputColumns )
44 , mNumOutputRows( nOutputRows )
45 , mRasterEntries( rasterEntries )
48 mOutputCrs = mRasterEntries.
at( 0 ).raster->crs();
53 : mFormulaString( formulaString )
54 , mOutputFile( outputFile )
55 , mOutputFormat( outputFormat )
56 , mOutputRectangle( outputExtent )
57 , mOutputCrs( outputCrs )
58 , mNumOutputColumns( nOutputColumns )
59 , mNumOutputRows( nOutputRows )
60 , mRasterEntries( rasterEntries )
82 for ( ; it != mRasterEntries.
constEnd(); ++it )
91 if ( it->raster->crs() != mOutputCrs )
94 proj.
setCRS( it->raster->crs(), mOutputCrs );
95 proj.
setInput( it->raster->dataProvider() );
98 block = proj.
block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows );
102 block = it->raster->dataProvider()->block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows );
104 inputBlocks.
insert( it->ref, block );
108 GDALDriverH outputDriver = openOutputDriver();
109 if ( outputDriver == NULL )
114 GDALDatasetH outputDataset = openOutputFile( outputDriver );
116 GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset, 1 );
118 float outputNodataValue = -FLT_MAX;
119 GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
130 for (
int i = 0; i < mNumOutputRows; ++i )
142 if ( calcNode->
calculate( inputBlocks, resultMatrix, i ) )
144 bool resultIsNumber = resultMatrix.
isNumber();
145 float* calcData =
new float[mNumOutputColumns];
147 for (
int j = 0; j < mNumOutputColumns; ++j )
149 calcData[j] = ( float )( resultIsNumber ? resultMatrix.
number() : resultMatrix.
data()[j] );
153 if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
155 qWarning(
"RasterIO error!" );
170 qDeleteAll( inputBlocks );
176 GDALDeleteDataset( outputDriver,
TO8F( mOutputFile ) );
179 GDALClose( outputDataset );
184 QgsRasterCalculator::QgsRasterCalculator()
185 : mNumOutputColumns( 0 )
186 , mNumOutputRows( 0 )
190 GDALDriverH QgsRasterCalculator::openOutputDriver()
192 char **driverMetadata;
195 GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.
toLocal8Bit().
data() );
197 if ( outputDriver == NULL )
202 driverMetadata = GDALGetMetadata( outputDriver, NULL );
203 if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE,
false ) )
211 GDALDatasetH QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver )
214 char **papszOptions = NULL;
215 GDALDatasetH outputDataset = GDALCreate( outputDriver,
TO8F( mOutputFile ), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions );
216 if ( outputDataset == NULL )
218 return outputDataset;
222 double geotransform[6];
223 outputGeoTransform( geotransform );
224 GDALSetGeoTransform( outputDataset, geotransform );
226 return outputDataset;
229 void QgsRasterCalculator::outputGeoTransform(
double* transform )
const
231 transform[0] = mOutputRectangle.
xMinimum();
232 transform[1] = mOutputRectangle.
width() / mNumOutputColumns;
234 transform[3] = mOutputRectangle.
yMaximum();
236 transform[5] = -mOutputRectangle.
height() / mNumOutputRows;
A rectangle specified with double values.
bool calculate(QMap< QString, QgsRasterBlock * > &rasterData, QgsRasterMatrix &result, int row=-1) const
Calculates result of raster calculation (might be real matrix or single number).
void setCRS(const QgsCoordinateReferenceSystem &theSrcCRS, const QgsCoordinateReferenceSystem &theDestCRS, int srcDatumTransform=-1, int destDatumTransform=-1)
set source and destination CRS
void setMaximum(int maximum)
double yMaximum() const
Get the y maximum value (top side of rectangle)
void setNodataValue(double d)
const_iterator constEnd() const
QString toWkt() const
A helper to get an wkt representation of this srs.
int processCalculation(QProgressDialog *p=0)
Starts the calculation and writes new raster.
QgsRasterCalculator(const QString &formulaString, const QString &outputFile, const QString &outputFormat, const QgsRectangle &outputExtent, int nOutputColumns, int nOutputRows, const QVector< QgsRasterCalculatorEntry > &rasterEntries)
QgsRasterCalculator constructor.
void setValue(int progress)
QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height) override
Read block of data using given extent and size.
QByteArray toLocal8Bit() const
bool isNumber() const
Returns true if matrix is 1x1 (=scalar number)
void setPrecision(Precision precision)
virtual bool setInput(QgsRasterInterface *input)
Set input.
const T & at(int i) const
const_iterator constBegin() const
Class for storing a coordinate reference system (CRS)
double * data()
Returns data array (but not ownership)
static QgsRasterCalcNode * parseRasterCalcString(const QString &str, QString &parserErrorMsg)
iterator insert(const Key &key, const T &value)
double width() const
Width of the rectangle.
double xMinimum() const
Get the x minimum value (left side of rectangle)
double height() const
Height of the rectangle.