Quantum GIS API Documentation  1.8
src/analysis/raster/qgsrelief.cpp
Go to the documentation of this file.
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines