Quantum GIS API Documentation
1.8
|
00001 /*************************************************************************** 00002 qgsrelief.cpp - description 00003 --------------------------- 00004 begin : November 2011 00005 copyright : (C) 2011 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 "qgsrelief.h" 00019 #include "qgsaspectfilter.h" 00020 #include "qgshillshadefilter.h" 00021 #include "qgsslopefilter.h" 00022 #include "qgis.h" 00023 #include "cpl_string.h" 00024 #include <QProgressDialog> 00025 #include <cfloat> 00026 00027 #include <QFile> 00028 #include <QTextStream> 00029 00030 #if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800 00031 #define TO8(x) (x).toUtf8().constData() 00032 #else 00033 #define TO8(x) (x).toLocal8Bit().constData() 00034 #endif 00035 00036 QgsRelief::QgsRelief( const QString& inputFile, const QString& outputFile, const QString& outputFormat ): \ 00037 mInputFile( inputFile ), mOutputFile( outputFile ), mOutputFormat( outputFormat ), mZFactor( 1.0 ) 00038 { 00039 mSlopeFilter = new QgsSlopeFilter( inputFile, outputFile, outputFormat ); 00040 mAspectFilter = new QgsAspectFilter( inputFile, outputFile, outputFormat ); 00041 mHillshadeFilter285 = new QgsHillshadeFilter( inputFile, outputFile, outputFormat, 285, 30 ); 00042 mHillshadeFilter300 = new QgsHillshadeFilter( inputFile, outputFile, outputFormat, 300, 30 ); 00043 mHillshadeFilter315 = new QgsHillshadeFilter( inputFile, outputFile, outputFormat, 315, 30 ); 00044 00045 /*mReliefColors = calculateOptimizedReliefClasses(); 00046 setDefaultReliefColors();*/ 00047 } 00048 00049 QgsRelief::~QgsRelief() 00050 { 00051 delete mSlopeFilter; 00052 delete mAspectFilter; 00053 delete mHillshadeFilter285; 00054 delete mHillshadeFilter300; 00055 delete mHillshadeFilter315; 00056 } 00057 00058 void QgsRelief::clearReliefColors() 00059 { 00060 mReliefColors.clear(); 00061 } 00062 00063 void QgsRelief::addReliefColorClass( const ReliefColor& color ) 00064 { 00065 mReliefColors.push_back( color ); 00066 } 00067 00068 void QgsRelief::setDefaultReliefColors() 00069 { 00070 clearReliefColors(); 00071 addReliefColorClass( ReliefColor( QColor( 9, 176, 76 ), 0, 200 ) ); 00072 addReliefColorClass( ReliefColor( QColor( 20, 228, 128 ), 200, 500 ) ); 00073 addReliefColorClass( ReliefColor( QColor( 167, 239, 153 ), 500, 1000 ) ); 00074 addReliefColorClass( ReliefColor( QColor( 218, 188, 143 ), 1000, 2000 ) ); 00075 addReliefColorClass( ReliefColor( QColor( 233, 158, 91 ), 2000, 4000 ) ); 00076 addReliefColorClass( ReliefColor( QColor( 255, 255, 255 ), 4000, 9000 ) ); 00077 } 00078 00079 int QgsRelief::processRaster( QProgressDialog* p ) 00080 { 00081 //open input file 00082 int xSize, ySize; 00083 GDALDatasetH inputDataset = openInputFile( xSize, ySize ); 00084 if ( inputDataset == NULL ) 00085 { 00086 return 1; //opening of input file failed 00087 } 00088 00089 //output driver 00090 GDALDriverH outputDriver = openOutputDriver(); 00091 if ( outputDriver == 0 ) 00092 { 00093 return 2; 00094 } 00095 00096 GDALDatasetH outputDataset = openOutputFile( inputDataset, outputDriver ); 00097 if ( outputDataset == NULL ) 00098 { 00099 return 3; //create operation on output file failed 00100 } 00101 00102 //initialize dependency filters with cell sizes 00103 mHillshadeFilter285->setCellSizeX( mCellSizeX ); 00104 mHillshadeFilter285->setCellSizeY( mCellSizeY ); 00105 mHillshadeFilter285->setZFactor( mZFactor ); 00106 mHillshadeFilter300->setCellSizeX( mCellSizeX ); 00107 mHillshadeFilter300->setCellSizeY( mCellSizeY ); 00108 mHillshadeFilter300->setZFactor( mZFactor ); 00109 mHillshadeFilter315->setCellSizeX( mCellSizeX ); 00110 mHillshadeFilter315->setCellSizeY( mCellSizeY ); 00111 mHillshadeFilter315->setZFactor( mZFactor ); 00112 mSlopeFilter->setCellSizeX( mCellSizeX ); 00113 mSlopeFilter->setCellSizeY( mCellSizeY ); 00114 mSlopeFilter->setZFactor( mZFactor ); 00115 mAspectFilter->setCellSizeX( mCellSizeX ); 00116 mAspectFilter->setCellSizeY( mCellSizeY ); 00117 mAspectFilter->setZFactor( mZFactor ); 00118 00119 //open first raster band for reading (operation is only for single band raster) 00120 GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, 1 ); 00121 if ( rasterBand == NULL ) 00122 { 00123 GDALClose( inputDataset ); 00124 GDALClose( outputDataset ); 00125 return 4; 00126 } 00127 mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, NULL ); 00128 mSlopeFilter->setInputNodataValue( mInputNodataValue ); 00129 mAspectFilter->setInputNodataValue( mInputNodataValue ); 00130 mHillshadeFilter285->setInputNodataValue( mInputNodataValue ); 00131 mHillshadeFilter300->setInputNodataValue( mInputNodataValue ); 00132 mHillshadeFilter315->setInputNodataValue( mInputNodataValue ); 00133 00134 GDALRasterBandH outputRedBand = GDALGetRasterBand( outputDataset, 1 ); 00135 GDALRasterBandH outputGreenBand = GDALGetRasterBand( outputDataset, 2 ); 00136 GDALRasterBandH outputBlueBand = GDALGetRasterBand( outputDataset, 3 ); 00137 00138 if ( outputRedBand == NULL || outputGreenBand == NULL || outputBlueBand == NULL ) 00139 { 00140 GDALClose( inputDataset ); 00141 GDALClose( outputDataset ); 00142 return 5; 00143 } 00144 //try to set -9999 as nodata value 00145 GDALSetRasterNoDataValue( outputRedBand, -9999 ); 00146 GDALSetRasterNoDataValue( outputGreenBand, -9999 ); 00147 GDALSetRasterNoDataValue( outputBlueBand, -9999 ); 00148 mOutputNodataValue = GDALGetRasterNoDataValue( outputRedBand, NULL ); 00149 mSlopeFilter->setOutputNodataValue( mOutputNodataValue ); 00150 mAspectFilter->setOutputNodataValue( mOutputNodataValue ); 00151 mHillshadeFilter285->setOutputNodataValue( mOutputNodataValue ); 00152 mHillshadeFilter300->setOutputNodataValue( mOutputNodataValue ); 00153 mHillshadeFilter315->setOutputNodataValue( mOutputNodataValue ); 00154 00155 if ( ySize < 3 ) //we require at least three rows (should be true for most datasets) 00156 { 00157 GDALClose( inputDataset ); 00158 GDALClose( outputDataset ); 00159 return 6; 00160 } 00161 00162 //keep only three scanlines in memory at a time 00163 float* scanLine1 = ( float * ) CPLMalloc( sizeof( float ) * xSize ); 00164 float* scanLine2 = ( float * ) CPLMalloc( sizeof( float ) * xSize ); 00165 float* scanLine3 = ( float * ) CPLMalloc( sizeof( float ) * xSize ); 00166 00167 unsigned char* resultRedLine = ( unsigned char * ) CPLMalloc( sizeof( unsigned char ) * xSize ); 00168 unsigned char* resultGreenLine = ( unsigned char * ) CPLMalloc( sizeof( unsigned char ) * xSize ); 00169 unsigned char* resultBlueLine = ( unsigned char * ) CPLMalloc( sizeof( unsigned char ) * xSize ); 00170 00171 if ( p ) 00172 { 00173 p->setMaximum( ySize ); 00174 } 00175 00176 bool resultOk; 00177 00178 //values outside the layer extent (if the 3x3 window is on the border) are sent to the processing method as (input) nodata values 00179 for ( int i = 0; i < ySize; ++i ) 00180 { 00181 if ( p ) 00182 { 00183 p->setValue( i ); 00184 } 00185 00186 if ( p && p->wasCanceled() ) 00187 { 00188 break; 00189 } 00190 00191 if ( i == 0 ) 00192 { 00193 //fill scanline 1 with (input) nodata for the values above the first row and feed scanline2 with the first row 00194 for ( int a = 0; a < xSize; ++a ) 00195 { 00196 scanLine1[a] = mInputNodataValue; 00197 } 00198 GDALRasterIO( rasterBand, GF_Read, 0, 0, xSize, 1, scanLine2, xSize, 1, GDT_Float32, 0, 0 ); 00199 } 00200 else 00201 { 00202 //normally fetch only scanLine3 and release scanline 1 if we move forward one row 00203 CPLFree( scanLine1 ); 00204 scanLine1 = scanLine2; 00205 scanLine2 = scanLine3; 00206 scanLine3 = ( float * ) CPLMalloc( sizeof( float ) * xSize ); 00207 } 00208 00209 if ( i == ySize - 1 ) //fill the row below the bottom with nodata values 00210 { 00211 for ( int a = 0; a < xSize; ++a ) 00212 { 00213 scanLine3[a] = mInputNodataValue; 00214 } 00215 } 00216 else 00217 { 00218 GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, scanLine3, xSize, 1, GDT_Float32, 0, 0 ); 00219 } 00220 00221 for ( int j = 0; j < xSize; ++j ) 00222 { 00223 if ( j == 0 ) 00224 { 00225 resultOk = processNineCellWindow( &mInputNodataValue, &scanLine1[j], &scanLine1[j+1], &mInputNodataValue, &scanLine2[j], \ 00226 &scanLine2[j+1], &mInputNodataValue, &scanLine3[j], &scanLine3[j+1], \ 00227 &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] ); 00228 } 00229 else if ( j == xSize - 1 ) 00230 { 00231 resultOk = processNineCellWindow( &scanLine1[j-1], &scanLine1[j], &mInputNodataValue, &scanLine2[j-1], &scanLine2[j], \ 00232 &mInputNodataValue, &scanLine3[j-1], &scanLine3[j], &mInputNodataValue, \ 00233 &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] ); 00234 } 00235 else 00236 { 00237 resultOk = processNineCellWindow( &scanLine1[j-1], &scanLine1[j], &scanLine1[j+1], &scanLine2[j-1], &scanLine2[j], \ 00238 &scanLine2[j+1], &scanLine3[j-1], &scanLine3[j], &scanLine3[j+1], \ 00239 &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] ); 00240 } 00241 00242 if ( !resultOk ) 00243 { 00244 resultRedLine[j] = mOutputNodataValue; 00245 resultGreenLine[j] = mOutputNodataValue; 00246 resultBlueLine[j] = mOutputNodataValue; 00247 } 00248 } 00249 00250 GDALRasterIO( outputRedBand, GF_Write, 0, i, xSize, 1, resultRedLine, xSize, 1, GDT_Byte, 0, 0 ); 00251 GDALRasterIO( outputGreenBand, GF_Write, 0, i, xSize, 1, resultGreenLine, xSize, 1, GDT_Byte, 0, 0 ); 00252 GDALRasterIO( outputBlueBand, GF_Write, 0, i, xSize, 1, resultBlueLine, xSize, 1, GDT_Byte, 0, 0 ); 00253 } 00254 00255 if ( p ) 00256 { 00257 p->setValue( ySize ); 00258 } 00259 00260 CPLFree( resultRedLine ); 00261 CPLFree( resultBlueLine ); 00262 CPLFree( resultGreenLine ); 00263 CPLFree( scanLine1 ); 00264 CPLFree( scanLine2 ); 00265 CPLFree( scanLine3 ); 00266 00267 GDALClose( inputDataset ); 00268 00269 if ( p && p->wasCanceled() ) 00270 { 00271 //delete the dataset without closing (because it is faster) 00272 GDALDeleteDataset( outputDriver, mOutputFile.toLocal8Bit().data() ); 00273 return 7; 00274 } 00275 GDALClose( outputDataset ); 00276 00277 return 0; 00278 } 00279 00280 bool QgsRelief::processNineCellWindow( float* x1, float* x2, float* x3, float* x4, float* x5, float* x6, float* x7, float* x8, float* x9, 00281 unsigned char* red, unsigned char* green, unsigned char* blue ) 00282 { 00283 //1. component: color and hillshade from 300 degrees 00284 int r = 0; 00285 int g = 0; 00286 int b = 0; 00287 00288 float hillShadeValue300 = mHillshadeFilter300->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 ); 00289 if ( hillShadeValue300 != mOutputNodataValue ) 00290 { 00291 if ( !setElevationColor( *x5, &r, &g, &b ) ) 00292 { 00293 r = hillShadeValue300; 00294 g = hillShadeValue300; 00295 b = hillShadeValue300; 00296 } 00297 else 00298 { 00299 r = r / 2.0 + hillShadeValue300 / 2.0; 00300 g = g / 2.0 + hillShadeValue300 / 2.0; 00301 b = b / 2.0 + hillShadeValue300 / 2.0; 00302 } 00303 } 00304 00305 //2. component: hillshade and slope 00306 float hillShadeValue315 = mHillshadeFilter315->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 ); 00307 float slope = mSlopeFilter->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 ); 00308 if ( hillShadeValue315 != mOutputNodataValue && slope != mOutputNodataValue ) 00309 { 00310 int r2, g2, b2; 00311 if ( slope > 15 ) 00312 { 00313 r2 = 0 / 2.0 + hillShadeValue315 / 2.0; 00314 g2 = 0 / 2.0 + hillShadeValue315 / 2.0; 00315 b2 = 0 / 2.0 + hillShadeValue315 / 2.0; 00316 } 00317 else if ( slope >= 1 ) 00318 { 00319 int slopeValue = 255 - ( slope / 15.0 * 255.0 ); 00320 r2 = slopeValue / 2.0 + hillShadeValue315 / 2.0; 00321 g2 = slopeValue / 2.0 + hillShadeValue315 / 2.0; 00322 b2 = slopeValue / 2.0 + hillShadeValue315 / 2.0; 00323 } 00324 else 00325 { 00326 r2 = hillShadeValue315; g2 = hillShadeValue315; b2 = hillShadeValue315; 00327 } 00328 00329 //combine with r,g,b with 70 percentage coverage 00330 r = r * 0.7 + r2 * 0.3; 00331 g = g * 0.7 + g2 * 0.3; 00332 b = b * 0.7 + b2 * 0.3; 00333 } 00334 00335 //3. combine yellow aspect with 10% transparency, illumination from 285 degrees 00336 float hillShadeValue285 = mHillshadeFilter285->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 ); 00337 float aspect = mAspectFilter->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 ); 00338 if ( hillShadeValue285 != mOutputNodataValue && aspect != mOutputNodataValue ) 00339 { 00340 double angle_diff = qAbs( 285 - aspect ); 00341 if ( angle_diff > 180 ) 00342 { 00343 angle_diff -= 180; 00344 } 00345 00346 int r3, g3, b3; 00347 if ( angle_diff < 90 ) 00348 { 00349 int aspectVal = ( 1 - cos( angle_diff * M_PI / 180 ) ) * 255; 00350 r3 = 0.5 * 255 + hillShadeValue315 * 0.5; 00351 g3 = 0.5 * 255 + hillShadeValue315 * 0.5; 00352 b3 = 0.5 * aspectVal + hillShadeValue315 * 0.5; 00353 } 00354 else //white 00355 { 00356 r3 = 0.5 * 255 + hillShadeValue315 * 0.5; 00357 g3 = 0.5 * 255 + hillShadeValue315 * 0.5; 00358 b3 = 0.5 * 255 + hillShadeValue315 * 0.5; 00359 } 00360 00361 r = r3 * 0.1 + r * 0.9; 00362 g = g3 * 0.1 + g * 0.9; 00363 b = b3 * 0.1 + b * 0.9; 00364 } 00365 00366 *red = ( unsigned char )r; 00367 *green = ( unsigned char )g; 00368 *blue = ( unsigned char )b; 00369 return true; 00370 } 00371 00372 bool QgsRelief::setElevationColor( double elevation, int* red, int* green, int* blue ) 00373 { 00374 QList< ReliefColor >::const_iterator reliefColorIt = mReliefColors.constBegin(); 00375 for ( ; reliefColorIt != mReliefColors.constEnd(); ++reliefColorIt ) 00376 { 00377 if ( elevation >= reliefColorIt->minElevation && elevation <= reliefColorIt->maxElevation ) 00378 { 00379 const QColor& c = reliefColorIt->color; 00380 *red = c.red(); 00381 *green = c.green(); 00382 *blue = c.blue(); 00383 00384 return true; 00385 } 00386 } 00387 return false; 00388 } 00389 00390 //duplicated from QgsNineCellFilter. Todo: make common base class 00391 GDALDatasetH QgsRelief::openInputFile( int& nCellsX, int& nCellsY ) 00392 { 00393 GDALDatasetH inputDataset = GDALOpen( TO8( mInputFile ), GA_ReadOnly ); 00394 if ( inputDataset != NULL ) 00395 { 00396 nCellsX = GDALGetRasterXSize( inputDataset ); 00397 nCellsY = GDALGetRasterYSize( inputDataset ); 00398 00399 //we need at least one band 00400 if ( GDALGetRasterCount( inputDataset ) < 1 ) 00401 { 00402 GDALClose( inputDataset ); 00403 return NULL; 00404 } 00405 } 00406 return inputDataset; 00407 } 00408 00409 GDALDriverH QgsRelief::openOutputDriver() 00410 { 00411 char **driverMetadata; 00412 00413 //open driver 00414 GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() ); 00415 00416 if ( outputDriver == NULL ) 00417 { 00418 return outputDriver; //return NULL, driver does not exist 00419 } 00420 00421 driverMetadata = GDALGetMetadata( outputDriver, NULL ); 00422 if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) ) 00423 { 00424 return NULL; //driver exist, but it does not support the create operation 00425 } 00426 00427 return outputDriver; 00428 } 00429 00430 GDALDatasetH QgsRelief::openOutputFile( GDALDatasetH inputDataset, GDALDriverH outputDriver ) 00431 { 00432 if ( inputDataset == NULL ) 00433 { 00434 return NULL; 00435 } 00436 00437 int xSize = GDALGetRasterXSize( inputDataset ); 00438 int ySize = GDALGetRasterYSize( inputDataset );; 00439 00440 //open output file 00441 char **papszOptions = NULL; 00442 00443 //use PACKBITS compression for tiffs by default 00444 papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "PACKBITS" ); 00445 00446 //create three band raster (reg, green, blue) 00447 GDALDatasetH outputDataset = GDALCreate( outputDriver, mOutputFile.toLocal8Bit().data(), xSize, ySize, 3, GDT_Byte, papszOptions ); 00448 if ( outputDataset == NULL ) 00449 { 00450 return outputDataset; 00451 } 00452 00453 //get geotransform from inputDataset 00454 double geotransform[6]; 00455 if ( GDALGetGeoTransform( inputDataset, geotransform ) != CE_None ) 00456 { 00457 GDALClose( outputDataset ); 00458 return NULL; 00459 } 00460 GDALSetGeoTransform( outputDataset, geotransform ); 00461 00462 //make sure mCellSizeX and mCellSizeY are always > 0 00463 mCellSizeX = geotransform[1]; 00464 if ( mCellSizeX < 0 ) 00465 { 00466 mCellSizeX = -mCellSizeX; 00467 } 00468 mCellSizeY = geotransform[5]; 00469 if ( mCellSizeY < 0 ) 00470 { 00471 mCellSizeY = -mCellSizeY; 00472 } 00473 00474 const char* projection = GDALGetProjectionRef( inputDataset ); 00475 GDALSetProjection( outputDataset, projection ); 00476 00477 return outputDataset; 00478 } 00479 00480 //this function is mainly there for debugging 00481 bool QgsRelief::exportFrequencyDistributionToCsv( const QString& file ) 00482 { 00483 int nCellsX, nCellsY; 00484 GDALDatasetH inputDataset = openInputFile( nCellsX, nCellsY ); 00485 if ( inputDataset == NULL ) 00486 { 00487 return false; 00488 } 00489 00490 //open first raster band for reading (elevation raster is always single band) 00491 GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset, 1 ); 00492 if ( elevationBand == NULL ) 00493 { 00494 GDALClose( inputDataset ); 00495 return false; 00496 } 00497 00498 //1. get minimum and maximum of elevation raster -> 252 elevation classes 00499 int minOk, maxOk; 00500 double minMax[2]; 00501 minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk ); 00502 minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk ); 00503 00504 if ( !minOk || !maxOk ) 00505 { 00506 GDALComputeRasterMinMax( elevationBand, TRUE, minMax ); 00507 } 00508 00509 //2. go through raster cells and get frequency of classes 00510 00511 //store elevation frequency in 256 elevation classes 00512 double frequency[252]; 00513 double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0; 00514 //initialize to zero 00515 for ( int i = 0; i < 252; ++i ) 00516 { 00517 frequency[i] = 0; 00518 } 00519 00520 float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX ); 00521 int elevationClass = -1; 00522 00523 for ( int i = 0; i < nCellsY; ++i ) 00524 { 00525 GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1, 00526 scanLine, nCellsX, 1, GDT_Float32, 00527 0, 0 ); 00528 for ( int j = 0; j < nCellsX; ++j ) 00529 { 00530 elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange ); 00531 if ( elevationClass >= 0 ) 00532 { 00533 frequency[elevationClass] += 1.0; 00534 } 00535 } 00536 } 00537 00538 CPLFree( scanLine ); 00539 00540 //log10 transformation for all frequency values 00541 for ( int i = 0; i < 252; ++i ) 00542 { 00543 frequency[i] = log10( frequency[i] ); 00544 } 00545 00546 //write out frequency values to csv file for debugging 00547 QFile outFile( file ); 00548 if ( !outFile.open( QIODevice::WriteOnly ) ) 00549 { 00550 return false; 00551 } 00552 00553 QTextStream outstream( &outFile ); 00554 for ( int i = 0; i < 252; ++i ) 00555 { 00556 outstream << QString::number( i ) + "," + QString::number( frequency[i] ) << endl; 00557 } 00558 outFile.close(); 00559 return true; 00560 } 00561 00562 QList< QgsRelief::ReliefColor > QgsRelief::calculateOptimizedReliefClasses() 00563 { 00564 QList< QgsRelief::ReliefColor > resultList; 00565 00566 int nCellsX, nCellsY; 00567 GDALDatasetH inputDataset = openInputFile( nCellsX, nCellsY ); 00568 if ( inputDataset == NULL ) 00569 { 00570 return resultList; 00571 } 00572 00573 //open first raster band for reading (elevation raster is always single band) 00574 GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset, 1 ); 00575 if ( elevationBand == NULL ) 00576 { 00577 GDALClose( inputDataset ); 00578 return resultList; 00579 } 00580 00581 //1. get minimum and maximum of elevation raster -> 252 elevation classes 00582 int minOk, maxOk; 00583 double minMax[2]; 00584 minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk ); 00585 minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk ); 00586 00587 if ( !minOk || !maxOk ) 00588 { 00589 GDALComputeRasterMinMax( elevationBand, TRUE, minMax ); 00590 } 00591 00592 //2. go through raster cells and get frequency of classes 00593 00594 //store elevation frequency in 256 elevation classes 00595 double frequency[252]; 00596 double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0; 00597 //initialize to zero 00598 for ( int i = 0; i < 252; ++i ) 00599 { 00600 frequency[i] = 0; 00601 } 00602 00603 float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX ); 00604 int elevationClass = -1; 00605 00606 for ( int i = 0; i < nCellsY; ++i ) 00607 { 00608 GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1, 00609 scanLine, nCellsX, 1, GDT_Float32, 00610 0, 0 ); 00611 for ( int j = 0; j < nCellsX; ++j ) 00612 { 00613 elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange ); 00614 if ( elevationClass < 0 ) 00615 { 00616 elevationClass = 0; 00617 } 00618 else if ( elevationClass >= 252 ) 00619 { 00620 elevationClass = 251; 00621 } 00622 frequency[elevationClass] += 1.0; 00623 } 00624 } 00625 00626 CPLFree( scanLine ); 00627 00628 //log10 transformation for all frequency values 00629 for ( int i = 0; i < 252; ++i ) 00630 { 00631 frequency[i] = log10( frequency[i] ); 00632 } 00633 00634 //start with 9 uniformly distributed classes 00635 QList<int> classBreaks; 00636 classBreaks.append( 0 ); 00637 classBreaks.append( 28 ); 00638 classBreaks.append( 56 ); 00639 classBreaks.append( 84 ); 00640 classBreaks.append( 112 ); 00641 classBreaks.append( 140 ); 00642 classBreaks.append( 168 ); 00643 classBreaks.append( 196 ); 00644 classBreaks.append( 224 ); 00645 classBreaks.append( 252 ); 00646 00647 for ( int i = 0; i < 10; ++i ) 00648 { 00649 optimiseClassBreaks( classBreaks, frequency ); 00650 } 00651 00652 //debug, print out all the classbreaks 00653 for ( int i = 0; i < classBreaks.size(); ++i ) 00654 { 00655 qWarning( "%d", classBreaks[i] ); 00656 } 00657 00658 //set colors according to optimised class breaks 00659 QList<QColor> colorList; 00660 colorList.push_back( QColor( 7, 165, 144 ) ); 00661 colorList.push_back( QColor( 12, 221, 162 ) ); 00662 colorList.push_back( QColor( 33, 252, 183 ) ); 00663 colorList.push_back( QColor( 247, 252, 152 ) ); 00664 colorList.push_back( QColor( 252, 196, 8 ) ); 00665 colorList.push_back( QColor( 252, 166, 15 ) ); 00666 colorList.push_back( QColor( 175, 101, 15 ) ); 00667 colorList.push_back( QColor( 255, 133, 92 ) ); 00668 colorList.push_back( QColor( 204, 204, 204 ) ); 00669 00670 for ( int i = 1; i < classBreaks.size(); ++i ) 00671 { 00672 double minElevation = minMax[0] + classBreaks[i - 1] * frequencyClassRange; 00673 double maxElevation = minMax[0] + classBreaks[i] * frequencyClassRange; 00674 resultList.push_back( QgsRelief::ReliefColor( colorList.at( i - 1 ), minElevation, maxElevation ) ); 00675 } 00676 00677 return resultList; 00678 } 00679 00680 void QgsRelief::optimiseClassBreaks( QList<int>& breaks, double* frequencies ) 00681 { 00682 int nClasses = breaks.size() - 1; 00683 double* a = new double[nClasses]; //slopes 00684 double* b = new double[nClasses]; //y-offsets 00685 00686 for ( int i = 0; i < nClasses; ++i ) 00687 { 00688 //get all the values between the class breaks into input 00689 QList< QPair < int, double > > regressionInput; 00690 for ( int j = breaks.at( i ); j < breaks.at( i + 1 ); ++j ) 00691 { 00692 regressionInput.push_back( qMakePair( j, frequencies[j] ) ); 00693 } 00694 00695 double aParam, bParam; 00696 if ( regressionInput.size() > 0 && calculateRegression( regressionInput, aParam, bParam ) ) 00697 { 00698 a[i] = aParam; 00699 b[i] = bParam; 00700 } 00701 else 00702 { 00703 a[i] = 0; 00704 b[i] = 0; //better default value 00705 } 00706 } 00707 00708 QList<int> classesToRemove; 00709 00710 //shift class boundaries or eliminate classes which fall together 00711 for ( int i = 1; i < nClasses ; ++i ) 00712 { 00713 if ( breaks[i] == breaks[ i - 1 ] ) 00714 { 00715 continue; 00716 } 00717 00718 if ( doubleNear( a[i - 1 ], a[i] ) ) 00719 { 00720 continue; 00721 } 00722 else 00723 { 00724 int newX = ( b[i - 1] - b[ i ] ) / ( a[ i ] - a[ i - 1 ] ); 00725 00726 if ( newX <= breaks[i - 1] ) 00727 { 00728 newX = breaks[i - 1]; 00729 // classesToRemove.push_back( i );//remove this class later as it falls together with the preceding one 00730 } 00731 else if ( i < nClasses - 1 && newX >= breaks[i + 1] ) 00732 { 00733 newX = breaks[i + 1]; 00734 // classesToRemove.push_back( i );//remove this class later as it falls together with the next one 00735 } 00736 00737 breaks[i] = newX; 00738 } 00739 } 00740 00741 for ( int i = classesToRemove.size() - 1; i >= 0; --i ) 00742 { 00743 breaks.removeAt( classesToRemove.at( i ) ); 00744 } 00745 00746 delete[] a; 00747 delete[] b; 00748 } 00749 00750 int QgsRelief::frequencyClassForElevation( double elevation, double minElevation, double elevationClassRange ) 00751 { 00752 return ( elevation - minElevation ) / elevationClassRange; 00753 } 00754 00755 bool QgsRelief::calculateRegression( const QList< QPair < int, double > >& input, double& a, double& b ) 00756 { 00757 double xMean, yMean; 00758 double xSum = 0; 00759 double ySum = 0; 00760 QList< QPair < int, double > >::const_iterator inputIt = input.constBegin(); 00761 for ( ; inputIt != input.constEnd(); ++inputIt ) 00762 { 00763 xSum += inputIt->first; 00764 ySum += inputIt->second; 00765 } 00766 xMean = xSum / input.size(); 00767 yMean = ySum / input.size(); 00768 00769 double sumCounter = 0; 00770 double sumDenominator = 0; 00771 inputIt = input.constBegin(); 00772 for ( ; inputIt != input.constEnd(); ++inputIt ) 00773 { 00774 sumCounter += (( inputIt->first - xMean ) * ( inputIt->second - yMean ) ); 00775 sumDenominator += (( inputIt->first - xMean ) * ( inputIt->first - xMean ) ); 00776 } 00777 00778 a = sumCounter / sumDenominator; 00779 b = yMean - a * xMean; 00780 00781 return true; 00782 } 00783