Quantum GIS API Documentation
1.7.4
|
00001 /*************************************************************************** 00002 qgsrastercalculator.cpp - description 00003 ----------------------- 00004 begin : September 28th, 2010 00005 copyright : (C) 2010 by Marco Hugentobler 00006 email : marco dot hugentobler at sourcepole dot ch 00007 ***************************************************************************/ 00008 00009 /*************************************************************************** 00010 * * 00011 * This program is free software; you can redistribute it and/or modify * 00012 * it under the terms of the GNU General Public License as published by * 00013 * the Free Software Foundation; either version 2 of the License, or * 00014 * (at your option) any later version. * 00015 * * 00016 ***************************************************************************/ 00017 00018 #include "qgsrastercalculator.h" 00019 #include "qgsrastercalcnode.h" 00020 #include "qgsrasterlayer.h" 00021 #include "qgsrastermatrix.h" 00022 #include "cpl_string.h" 00023 #include <QProgressDialog> 00024 00025 #include "gdalwarper.h" 00026 #include <ogr_srs_api.h> 00027 00028 #if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800 00029 #define TO8(x) (x).toUtf8().constData() 00030 #else 00031 #define TO8(x) (x).toLocal8Bit().constData() 00032 #endif 00033 00034 QgsRasterCalculator::QgsRasterCalculator( const QString& formulaString, const QString& outputFile, const QString& outputFormat, 00035 const QgsRectangle& outputExtent, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry>& rasterEntries ): mFormulaString( formulaString ), mOutputFile( outputFile ), mOutputFormat( outputFormat ), 00036 mOutputRectangle( outputExtent ), mNumOutputColumns( nOutputColumns ), mNumOutputRows( nOutputRows ), mRasterEntries( rasterEntries ) 00037 { 00038 } 00039 00040 QgsRasterCalculator::~QgsRasterCalculator() 00041 { 00042 } 00043 00044 int QgsRasterCalculator::processCalculation( QProgressDialog* p ) 00045 { 00046 //prepare search string / tree 00047 QString errorString; 00048 QgsRasterCalcNode* calcNode = QgsRasterCalcNode::parseRasterCalcString( mFormulaString, errorString ); 00049 if ( !calcNode ) 00050 { 00051 //error 00052 } 00053 00054 double targetGeoTransform[6]; 00055 outputGeoTransform( targetGeoTransform ); 00056 00057 //open all input rasters for reading 00058 QMap< QString, GDALRasterBandH > mInputRasterBands; //raster references and corresponding scanline data 00059 QMap< QString, QgsRasterMatrix* > inputScanLineData; //stores raster references and corresponding scanline data 00060 QVector< GDALDatasetH > mInputDatasets; //raster references and corresponding dataset 00061 00062 QVector<QgsRasterCalculatorEntry>::const_iterator it = mRasterEntries.constBegin(); 00063 for ( ; it != mRasterEntries.constEnd(); ++it ) 00064 { 00065 if ( !it->raster ) // no raster layer in entry 00066 { 00067 return 2; 00068 } 00069 GDALDatasetH inputDataset = GDALOpen( TO8( it->raster->source() ), GA_ReadOnly ); 00070 if ( inputDataset == NULL ) 00071 { 00072 return 2; 00073 } 00074 00075 //check if the input dataset is south up or rotated. If yes, use GDALAutoCreateWarpedVRT to create a north up raster 00076 double inputGeoTransform[6]; 00077 if ( GDALGetGeoTransform( inputDataset, inputGeoTransform ) == CE_None 00078 && ( inputGeoTransform[1] < 0.0 00079 || inputGeoTransform[2] != 0.0 00080 || inputGeoTransform[4] != 0.0 00081 || inputGeoTransform[5] > 0.0 ) ) 00082 { 00083 GDALDatasetH vDataset = GDALAutoCreateWarpedVRT( inputDataset, NULL, NULL, GRA_NearestNeighbour, 0.2, NULL ); 00084 mInputDatasets.push_back( vDataset ); 00085 mInputDatasets.push_back( inputDataset ); 00086 inputDataset = vDataset; 00087 } 00088 else 00089 { 00090 mInputDatasets.push_back( inputDataset ); 00091 } 00092 00093 00094 GDALRasterBandH inputRasterBand = GDALGetRasterBand( inputDataset, it->bandNumber ); 00095 if ( inputRasterBand == NULL ) 00096 { 00097 return 2; 00098 } 00099 00100 int nodataSuccess; 00101 double nodataValue = GDALGetRasterNoDataValue( inputRasterBand, &nodataSuccess ); 00102 00103 mInputRasterBands.insert( it->ref, inputRasterBand ); 00104 inputScanLineData.insert( it->ref, new QgsRasterMatrix( mNumOutputColumns, 1, new float[mNumOutputColumns], nodataValue ) ); 00105 } 00106 00107 //open output dataset for writing 00108 GDALDriverH outputDriver = openOutputDriver(); 00109 if ( outputDriver == NULL ) 00110 { 00111 return 1; 00112 } 00113 GDALDatasetH outputDataset = openOutputFile( outputDriver ); 00114 00115 //copy the projection info from the first input raster 00116 if ( mRasterEntries.size() > 0 ) 00117 { 00118 QgsRasterLayer* rl = mRasterEntries.at( 0 ).raster; 00119 if ( rl ) 00120 { 00121 char* crsWKT = 0; 00122 OGRSpatialReferenceH ogrSRS = OSRNewSpatialReference( NULL ); 00123 if ( OSRSetFromUserInput( ogrSRS, rl->crs().authid().toUtf8().constData() ) == OGRERR_NONE ) 00124 { 00125 OSRExportToWkt( ogrSRS, &crsWKT ); 00126 GDALSetProjection( outputDataset, crsWKT ); 00127 } 00128 else 00129 { 00130 GDALSetProjection( outputDataset, TO8( rl->crs().toWkt() ) ); 00131 } 00132 OSRDestroySpatialReference( ogrSRS ); 00133 CPLFree( crsWKT ); 00134 } 00135 } 00136 00137 00138 GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset, 1 ); 00139 00140 float outputNodataValue = -FLT_MAX; 00141 GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue ); 00142 00143 float* resultScanLine = ( float * ) CPLMalloc( sizeof( float ) * mNumOutputColumns ); 00144 00145 if ( p ) 00146 { 00147 p->setMaximum( mNumOutputRows ); 00148 } 00149 00150 QgsRasterMatrix resultMatrix; 00151 00152 //read / write line by line 00153 for ( int i = 0; i < mNumOutputRows; ++i ) 00154 { 00155 if ( p ) 00156 { 00157 p->setValue( i ); 00158 } 00159 00160 if ( p && p->wasCanceled() ) 00161 { 00162 break; 00163 } 00164 00165 //fill buffers 00166 QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin(); 00167 for ( ; bufferIt != inputScanLineData.end(); ++bufferIt ) 00168 { 00169 double sourceTransformation[6]; 00170 GDALRasterBandH sourceRasterBand = mInputRasterBands[bufferIt.key()]; 00171 GDALGetGeoTransform( GDALGetBandDataset( sourceRasterBand ), sourceTransformation ); 00172 //the function readRasterPart calls GDALRasterIO (and ev. does some conversion if raster transformations are not the same) 00173 readRasterPart( targetGeoTransform, 0, i, mNumOutputColumns, 1, sourceTransformation, sourceRasterBand, bufferIt.value()->data() ); 00174 } 00175 00176 if ( calcNode->calculate( inputScanLineData, resultMatrix ) ) 00177 { 00178 bool resultIsNumber = resultMatrix.isNumber(); 00179 float* calcData; 00180 00181 if ( resultIsNumber ) //scalar result. Insert number for every pixel 00182 { 00183 calcData = new float[mNumOutputColumns]; 00184 for ( int j = 0; j < mNumOutputColumns; ++j ) 00185 { 00186 calcData[j] = resultMatrix.number(); 00187 } 00188 } 00189 else //result is real matrix 00190 { 00191 calcData = resultMatrix.data(); 00192 } 00193 00194 //replace all matrix nodata values with output nodatas 00195 for ( int j = 0; j < mNumOutputColumns; ++j ) 00196 { 00197 if ( calcData[j] == resultMatrix.nodataValue() ) 00198 { 00199 calcData[j] = outputNodataValue; 00200 } 00201 } 00202 00203 //write scanline to the dataset 00204 if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None ) 00205 { 00206 qWarning( "RasterIO error!" ); 00207 } 00208 00209 if ( resultIsNumber ) 00210 { 00211 delete[] calcData; 00212 } 00213 } 00214 00215 } 00216 00217 if ( p ) 00218 { 00219 p->setValue( mNumOutputRows ); 00220 } 00221 00222 //close datasets and release memory 00223 delete calcNode; 00224 QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin(); 00225 for ( ; bufferIt != inputScanLineData.end(); ++bufferIt ) 00226 { 00227 delete bufferIt.value(); 00228 } 00229 inputScanLineData.clear(); 00230 00231 QVector< GDALDatasetH >::iterator datasetIt = mInputDatasets.begin(); 00232 for ( ; datasetIt != mInputDatasets.end(); ++ datasetIt ) 00233 { 00234 GDALClose( *datasetIt ); 00235 } 00236 00237 if ( p && p->wasCanceled() ) 00238 { 00239 //delete the dataset without closing (because it is faster) 00240 GDALDeleteDataset( outputDriver, mOutputFile.toLocal8Bit().data() ); 00241 return 3; 00242 } 00243 GDALClose( outputDataset ); 00244 CPLFree( resultScanLine ); 00245 return 0; 00246 } 00247 00248 QgsRasterCalculator::QgsRasterCalculator() 00249 { 00250 } 00251 00252 GDALDriverH QgsRasterCalculator::openOutputDriver() 00253 { 00254 char **driverMetadata; 00255 00256 //open driver 00257 GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() ); 00258 00259 if ( outputDriver == NULL ) 00260 { 00261 return outputDriver; //return NULL, driver does not exist 00262 } 00263 00264 driverMetadata = GDALGetMetadata( outputDriver, NULL ); 00265 if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) ) 00266 { 00267 return NULL; //driver exist, but it does not support the create operation 00268 } 00269 00270 return outputDriver; 00271 } 00272 00273 GDALDatasetH QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver ) 00274 { 00275 //open output file 00276 char **papszOptions = NULL; 00277 GDALDatasetH outputDataset = GDALCreate( outputDriver, mOutputFile.toLocal8Bit().data(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions ); 00278 if ( outputDataset == NULL ) 00279 { 00280 return outputDataset; 00281 } 00282 00283 //assign georef information 00284 double geotransform[6]; 00285 outputGeoTransform( geotransform ); 00286 GDALSetGeoTransform( outputDataset, geotransform ); 00287 00288 return outputDataset; 00289 } 00290 00291 void QgsRasterCalculator::readRasterPart( double* targetGeotransform, int xOffset, int yOffset, int nCols, int nRows, double* sourceTransform, GDALRasterBandH sourceBand, float* rasterBuffer ) 00292 { 00293 //If dataset transform is the same as the requested transform, do a normal GDAL raster io 00294 if ( transformationsEqual( targetGeotransform, sourceTransform ) ) 00295 { 00296 GDALRasterIO( sourceBand, GF_Read, xOffset, yOffset, nCols, nRows, rasterBuffer, nCols, nRows, GDT_Float32, 0, 0 ); 00297 return; 00298 } 00299 00300 //pixel calculation needed because of different raster position / resolution 00301 int nodataSuccess; 00302 double nodataValue = GDALGetRasterNoDataValue( sourceBand, &nodataSuccess ); 00303 QgsRectangle targetRect( targetGeotransform[0] + targetGeotransform[1] * xOffset, targetGeotransform[3] + yOffset * targetGeotransform[5] + nRows * targetGeotransform[5] 00304 , targetGeotransform[0] + targetGeotransform[1] * xOffset + targetGeotransform[1] * nCols, targetGeotransform[3] + yOffset * targetGeotransform[5] ); 00305 QgsRectangle sourceRect( sourceTransform[0], sourceTransform[3] + GDALGetRasterBandYSize( sourceBand ) * sourceTransform[5], 00306 sourceTransform[0] + GDALGetRasterBandXSize( sourceBand )* sourceTransform[1], sourceTransform[3] ); 00307 QgsRectangle intersection = targetRect.intersect( &sourceRect ); 00308 00309 //no intersection, fill all the pixels with nodata values 00310 if ( intersection.isEmpty() ) 00311 { 00312 int nPixels = nCols * nRows; 00313 for ( int i = 0; i < nPixels; ++i ) 00314 { 00315 rasterBuffer[i] = nodataValue; 00316 } 00317 return; 00318 } 00319 00320 //do raster io in source resolution 00321 int sourcePixelOffsetXMin = floor(( intersection.xMinimum() - sourceTransform[0] ) / sourceTransform[1] ); 00322 int sourcePixelOffsetXMax = ceil(( intersection.xMaximum() - sourceTransform[0] ) / sourceTransform[1] ); 00323 int nSourcePixelsX = sourcePixelOffsetXMax - sourcePixelOffsetXMin; 00324 int sourcePixelOffsetYMax = floor(( intersection.yMaximum() - sourceTransform[3] ) / sourceTransform[5] ); 00325 int sourcePixelOffsetYMin = ceil(( intersection.yMinimum() - sourceTransform[3] ) / sourceTransform[5] ); 00326 int nSourcePixelsY = sourcePixelOffsetYMin - sourcePixelOffsetYMax; 00327 float* sourceRaster = ( float * ) CPLMalloc( sizeof( float ) * nSourcePixelsX * nSourcePixelsY ); 00328 double sourceRasterXMin = sourceRect.xMinimum() + sourcePixelOffsetXMin * sourceTransform[1]; 00329 double sourceRasterYMax = sourceRect.yMaximum() + sourcePixelOffsetYMax * sourceTransform[5]; 00330 GDALRasterIO( sourceBand, GF_Read, sourcePixelOffsetXMin, sourcePixelOffsetYMax, nSourcePixelsX, nSourcePixelsY, 00331 sourceRaster, nSourcePixelsX, nSourcePixelsY, GDT_Float32, 0, 0 ); 00332 00333 00334 double targetPixelX; 00335 double targetPixelXMin = targetGeotransform[0] + targetGeotransform[1] * xOffset + targetGeotransform[1] / 2.0; 00336 double targetPixelY = targetGeotransform[3] + targetGeotransform[5] * yOffset + targetGeotransform[5] / 2.0; //coordinates of current target pixel 00337 int sourceIndexX, sourceIndexY; //current raster index in source pixels 00338 double sx, sy; 00339 for ( int i = 0; i < nRows; ++i ) 00340 { 00341 targetPixelX = targetPixelXMin; 00342 for ( int j = 0; j < nCols; ++j ) 00343 { 00344 sx = ( targetPixelX - sourceRasterXMin ) / sourceTransform[1]; 00345 sourceIndexX = sx > 0 ? sx : floor( sx ); 00346 sy = ( targetPixelY - sourceRasterYMax ) / sourceTransform[5]; 00347 sourceIndexY = sy > 0 ? sy : floor( sy ); 00348 if ( sourceIndexX >= 0 && sourceIndexX < nSourcePixelsX 00349 && sourceIndexY >= 0 && sourceIndexY < nSourcePixelsY ) 00350 { 00351 rasterBuffer[j + i*nRows] = sourceRaster[ sourceIndexX + nSourcePixelsX * sourceIndexY ]; 00352 } 00353 else 00354 { 00355 rasterBuffer[j + i*j] = nodataValue; 00356 } 00357 targetPixelX += targetGeotransform[1]; 00358 } 00359 targetPixelY += targetGeotransform[5]; 00360 } 00361 00362 CPLFree( sourceRaster ); 00363 return; 00364 } 00365 00366 bool QgsRasterCalculator::transformationsEqual( double* t1, double* t2 ) const 00367 { 00368 for ( int i = 0; i < 6; ++i ) 00369 { 00370 if ( !doubleNear( t1[i], t2[i], 0.00001 ) ) 00371 { 00372 return false; 00373 } 00374 } 00375 return true; 00376 } 00377 00378 void QgsRasterCalculator::outputGeoTransform( double* transform ) const 00379 { 00380 transform[0] = mOutputRectangle.xMinimum(); 00381 transform[1] = mOutputRectangle.width() / mNumOutputColumns; 00382 transform[2] = 0; 00383 transform[3] = mOutputRectangle.yMaximum(); 00384 transform[4] = 0; 00385 transform[5] = -mOutputRectangle.height() / mNumOutputRows; 00386 } 00387 00388