QGIS API Documentation  2.17.0-Master (973e4b0)
qgsrelief.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsrelief.cpp - description
3  ---------------------------
4  begin : November 2011
5  copyright : (C) 2011 by Marco Hugentobler
6  email : marco dot hugentobler at sourcepole dot ch
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #include "qgslogger.h"
19 #include "qgsrelief.h"
20 #include "qgsaspectfilter.h"
21 #include "qgshillshadefilter.h"
22 #include "qgsslopefilter.h"
23 #include "qgis.h"
24 #include "cpl_string.h"
25 #include <QProgressDialog>
26 #include <cfloat>
27 
28 #include <QFile>
29 #include <QTextStream>
30 
31 #if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800
32 #define TO8F(x) (x).toUtf8().constData()
33 #else
34 #define TO8F(x) QFile::encodeName( x ).constData()
35 #endif
36 
37 QgsRelief::QgsRelief( const QString& inputFile, const QString& outputFile, const QString& outputFormat )
38  : mInputFile( inputFile )
39  , mOutputFile( outputFile )
40  , mOutputFormat( outputFormat )
41  , mCellSizeX( 0.0 )
42  , mCellSizeY( 0.0 )
43  , mInputNodataValue( -1 )
44  , mOutputNodataValue( -1 )
45  , mZFactor( 1.0 )
46 {
47  mSlopeFilter = new QgsSlopeFilter( inputFile, outputFile, outputFormat );
48  mAspectFilter = new QgsAspectFilter( inputFile, outputFile, outputFormat );
49  mHillshadeFilter285 = new QgsHillshadeFilter( inputFile, outputFile, outputFormat, 285, 30 );
50  mHillshadeFilter300 = new QgsHillshadeFilter( inputFile, outputFile, outputFormat, 300, 30 );
51  mHillshadeFilter315 = new QgsHillshadeFilter( inputFile, outputFile, outputFormat, 315, 30 );
52 
53  /*mReliefColors = calculateOptimizedReliefClasses();
54  setDefaultReliefColors();*/
55 }
56 
58 {
59  delete mSlopeFilter;
60  delete mAspectFilter;
61  delete mHillshadeFilter285;
62  delete mHillshadeFilter300;
63  delete mHillshadeFilter315;
64 }
65 
67 {
68  mReliefColors.clear();
69 }
70 
72 {
73  mReliefColors.push_back( color );
74 }
75 
76 void QgsRelief::setDefaultReliefColors()
77 {
79  addReliefColorClass( ReliefColor( QColor( 9, 176, 76 ), 0, 200 ) );
80  addReliefColorClass( ReliefColor( QColor( 20, 228, 128 ), 200, 500 ) );
81  addReliefColorClass( ReliefColor( QColor( 167, 239, 153 ), 500, 1000 ) );
82  addReliefColorClass( ReliefColor( QColor( 218, 188, 143 ), 1000, 2000 ) );
83  addReliefColorClass( ReliefColor( QColor( 233, 158, 91 ), 2000, 4000 ) );
84  addReliefColorClass( ReliefColor( QColor( 255, 255, 255 ), 4000, 9000 ) );
85 }
86 
88 {
89  //open input file
90  int xSize, ySize;
91  GDALDatasetH inputDataset = openInputFile( xSize, ySize );
92  if ( !inputDataset )
93  {
94  return 1; //opening of input file failed
95  }
96 
97  //output driver
98  GDALDriverH outputDriver = openOutputDriver();
99  if ( !outputDriver )
100  {
101  return 2;
102  }
103 
104  GDALDatasetH outputDataset = openOutputFile( inputDataset, outputDriver );
105  if ( !outputDataset )
106  {
107  return 3; //create operation on output file failed
108  }
109 
110  //initialize dependency filters with cell sizes
111  mHillshadeFilter285->setCellSizeX( mCellSizeX );
112  mHillshadeFilter285->setCellSizeY( mCellSizeY );
113  mHillshadeFilter285->setZFactor( mZFactor );
114  mHillshadeFilter300->setCellSizeX( mCellSizeX );
115  mHillshadeFilter300->setCellSizeY( mCellSizeY );
116  mHillshadeFilter300->setZFactor( mZFactor );
117  mHillshadeFilter315->setCellSizeX( mCellSizeX );
118  mHillshadeFilter315->setCellSizeY( mCellSizeY );
119  mHillshadeFilter315->setZFactor( mZFactor );
120  mSlopeFilter->setCellSizeX( mCellSizeX );
121  mSlopeFilter->setCellSizeY( mCellSizeY );
122  mSlopeFilter->setZFactor( mZFactor );
123  mAspectFilter->setCellSizeX( mCellSizeX );
124  mAspectFilter->setCellSizeY( mCellSizeY );
125  mAspectFilter->setZFactor( mZFactor );
126 
127  //open first raster band for reading (operation is only for single band raster)
128  GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, 1 );
129  if ( !rasterBand )
130  {
131  GDALClose( inputDataset );
132  GDALClose( outputDataset );
133  return 4;
134  }
135  mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, nullptr );
136  mSlopeFilter->setInputNodataValue( mInputNodataValue );
137  mAspectFilter->setInputNodataValue( mInputNodataValue );
138  mHillshadeFilter285->setInputNodataValue( mInputNodataValue );
139  mHillshadeFilter300->setInputNodataValue( mInputNodataValue );
140  mHillshadeFilter315->setInputNodataValue( mInputNodataValue );
141 
142  GDALRasterBandH outputRedBand = GDALGetRasterBand( outputDataset, 1 );
143  GDALRasterBandH outputGreenBand = GDALGetRasterBand( outputDataset, 2 );
144  GDALRasterBandH outputBlueBand = GDALGetRasterBand( outputDataset, 3 );
145 
146  if ( !outputRedBand || !outputGreenBand || !outputBlueBand )
147  {
148  GDALClose( inputDataset );
149  GDALClose( outputDataset );
150  return 5;
151  }
152  //try to set -9999 as nodata value
153  GDALSetRasterNoDataValue( outputRedBand, -9999 );
154  GDALSetRasterNoDataValue( outputGreenBand, -9999 );
155  GDALSetRasterNoDataValue( outputBlueBand, -9999 );
156  mOutputNodataValue = GDALGetRasterNoDataValue( outputRedBand, nullptr );
157  mSlopeFilter->setOutputNodataValue( mOutputNodataValue );
158  mAspectFilter->setOutputNodataValue( mOutputNodataValue );
159  mHillshadeFilter285->setOutputNodataValue( mOutputNodataValue );
160  mHillshadeFilter300->setOutputNodataValue( mOutputNodataValue );
161  mHillshadeFilter315->setOutputNodataValue( mOutputNodataValue );
162 
163  if ( ySize < 3 ) //we require at least three rows (should be true for most datasets)
164  {
165  GDALClose( inputDataset );
166  GDALClose( outputDataset );
167  return 6;
168  }
169 
170  //keep only three scanlines in memory at a time
171  float* scanLine1 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
172  float* scanLine2 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
173  float* scanLine3 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
174 
175  unsigned char* resultRedLine = ( unsigned char * ) CPLMalloc( sizeof( unsigned char ) * xSize );
176  unsigned char* resultGreenLine = ( unsigned char * ) CPLMalloc( sizeof( unsigned char ) * xSize );
177  unsigned char* resultBlueLine = ( unsigned char * ) CPLMalloc( sizeof( unsigned char ) * xSize );
178 
179  if ( p )
180  {
181  p->setMaximum( ySize );
182  }
183 
184  bool resultOk;
185 
186  //values outside the layer extent (if the 3x3 window is on the border) are sent to the processing method as (input) nodata values
187  for ( int i = 0; i < ySize; ++i )
188  {
189  if ( p )
190  {
191  p->setValue( i );
192  }
193 
194  if ( p && p->wasCanceled() )
195  {
196  break;
197  }
198 
199  if ( i == 0 )
200  {
201  //fill scanline 1 with (input) nodata for the values above the first row and feed scanline2 with the first row
202  for ( int a = 0; a < xSize; ++a )
203  {
204  scanLine1[a] = mInputNodataValue;
205  }
206  if ( GDALRasterIO( rasterBand, GF_Read, 0, 0, xSize, 1, scanLine2, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
207  {
208  QgsDebugMsg( "Raster IO Error" );
209  }
210  }
211  else
212  {
213  //normally fetch only scanLine3 and release scanline 1 if we move forward one row
214  CPLFree( scanLine1 );
215  scanLine1 = scanLine2;
216  scanLine2 = scanLine3;
217  scanLine3 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
218  }
219 
220  if ( i == ySize - 1 ) //fill the row below the bottom with nodata values
221  {
222  for ( int a = 0; a < xSize; ++a )
223  {
224  scanLine3[a] = mInputNodataValue;
225  }
226  }
227  else
228  {
229  if ( GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, scanLine3, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
230  {
231  QgsDebugMsg( "Raster IO Error" );
232  }
233  }
234 
235  for ( int j = 0; j < xSize; ++j )
236  {
237  if ( j == 0 )
238  {
239  resultOk = processNineCellWindow( &mInputNodataValue, &scanLine1[j], &scanLine1[j+1], &mInputNodataValue, &scanLine2[j], \
240  &scanLine2[j+1], &mInputNodataValue, &scanLine3[j], &scanLine3[j+1], \
241  &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] );
242  }
243  else if ( j == xSize - 1 )
244  {
245  resultOk = processNineCellWindow( &scanLine1[j-1], &scanLine1[j], &mInputNodataValue, &scanLine2[j-1], &scanLine2[j], \
246  &mInputNodataValue, &scanLine3[j-1], &scanLine3[j], &mInputNodataValue, \
247  &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] );
248  }
249  else
250  {
251  resultOk = processNineCellWindow( &scanLine1[j-1], &scanLine1[j], &scanLine1[j+1], &scanLine2[j-1], &scanLine2[j], \
252  &scanLine2[j+1], &scanLine3[j-1], &scanLine3[j], &scanLine3[j+1], \
253  &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] );
254  }
255 
256  if ( !resultOk )
257  {
258  resultRedLine[j] = mOutputNodataValue;
259  resultGreenLine[j] = mOutputNodataValue;
260  resultBlueLine[j] = mOutputNodataValue;
261  }
262  }
263 
264  if ( GDALRasterIO( outputRedBand, GF_Write, 0, i, xSize, 1, resultRedLine, xSize, 1, GDT_Byte, 0, 0 ) != CE_None )
265  {
266  QgsDebugMsg( "Raster IO Error" );
267  }
268  if ( GDALRasterIO( outputGreenBand, GF_Write, 0, i, xSize, 1, resultGreenLine, xSize, 1, GDT_Byte, 0, 0 ) != CE_None )
269  {
270  QgsDebugMsg( "Raster IO Error" );
271  }
272  if ( GDALRasterIO( outputBlueBand, GF_Write, 0, i, xSize, 1, resultBlueLine, xSize, 1, GDT_Byte, 0, 0 ) != CE_None )
273  {
274  QgsDebugMsg( "Raster IO Error" );
275  }
276  }
277 
278  if ( p )
279  {
280  p->setValue( ySize );
281  }
282 
283  CPLFree( resultRedLine );
284  CPLFree( resultBlueLine );
285  CPLFree( resultGreenLine );
286  CPLFree( scanLine1 );
287  CPLFree( scanLine2 );
288  CPLFree( scanLine3 );
289 
290  GDALClose( inputDataset );
291 
292  if ( p && p->wasCanceled() )
293  {
294  //delete the dataset without closing (because it is faster)
295  GDALDeleteDataset( outputDriver, TO8F( mOutputFile ) );
296  return 7;
297  }
298  GDALClose( outputDataset );
299 
300  return 0;
301 }
302 
303 bool QgsRelief::processNineCellWindow( float* x1, float* x2, float* x3, float* x4, float* x5, float* x6, float* x7, float* x8, float* x9,
304  unsigned char* red, unsigned char* green, unsigned char* blue )
305 {
306  //1. component: color and hillshade from 300 degrees
307  int r = 0;
308  int g = 0;
309  int b = 0;
310 
311  float hillShadeValue300 = mHillshadeFilter300->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
312  if ( hillShadeValue300 != mOutputNodataValue )
313  {
314  if ( !setElevationColor( *x5, &r, &g, &b ) )
315  {
316  r = hillShadeValue300;
317  g = hillShadeValue300;
318  b = hillShadeValue300;
319  }
320  else
321  {
322  r = r / 2.0 + hillShadeValue300 / 2.0;
323  g = g / 2.0 + hillShadeValue300 / 2.0;
324  b = b / 2.0 + hillShadeValue300 / 2.0;
325  }
326  }
327 
328  //2. component: hillshade and slope
329  float hillShadeValue315 = mHillshadeFilter315->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
330  float slope = mSlopeFilter->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
331  if ( hillShadeValue315 != mOutputNodataValue && slope != mOutputNodataValue )
332  {
333  int r2, g2, b2;
334  if ( slope > 15 )
335  {
336  r2 = 0 / 2.0 + hillShadeValue315 / 2.0;
337  g2 = 0 / 2.0 + hillShadeValue315 / 2.0;
338  b2 = 0 / 2.0 + hillShadeValue315 / 2.0;
339  }
340  else if ( slope >= 1 )
341  {
342  int slopeValue = 255 - ( slope / 15.0 * 255.0 );
343  r2 = slopeValue / 2.0 + hillShadeValue315 / 2.0;
344  g2 = slopeValue / 2.0 + hillShadeValue315 / 2.0;
345  b2 = slopeValue / 2.0 + hillShadeValue315 / 2.0;
346  }
347  else
348  {
349  r2 = hillShadeValue315;
350  g2 = hillShadeValue315;
351  b2 = hillShadeValue315;
352  }
353 
354  //combine with r,g,b with 70 percentage coverage
355  r = r * 0.7 + r2 * 0.3;
356  g = g * 0.7 + g2 * 0.3;
357  b = b * 0.7 + b2 * 0.3;
358  }
359 
360  //3. combine yellow aspect with 10% transparency, illumination from 285 degrees
361  float hillShadeValue285 = mHillshadeFilter285->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
362  float aspect = mAspectFilter->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
363  if ( hillShadeValue285 != mOutputNodataValue && aspect != mOutputNodataValue )
364  {
365  double angle_diff = qAbs( 285 - aspect );
366  if ( angle_diff > 180 )
367  {
368  angle_diff -= 180;
369  }
370 
371  int r3, g3, b3;
372  if ( angle_diff < 90 )
373  {
374  int aspectVal = ( 1 - cos( angle_diff * M_PI / 180 ) ) * 255;
375  r3 = 0.5 * 255 + hillShadeValue315 * 0.5;
376  g3 = 0.5 * 255 + hillShadeValue315 * 0.5;
377  b3 = 0.5 * aspectVal + hillShadeValue315 * 0.5;
378  }
379  else //white
380  {
381  r3 = 0.5 * 255 + hillShadeValue315 * 0.5;
382  g3 = 0.5 * 255 + hillShadeValue315 * 0.5;
383  b3 = 0.5 * 255 + hillShadeValue315 * 0.5;
384  }
385 
386  r = r3 * 0.1 + r * 0.9;
387  g = g3 * 0.1 + g * 0.9;
388  b = b3 * 0.1 + b * 0.9;
389  }
390 
391  *red = ( unsigned char )r;
392  *green = ( unsigned char )g;
393  *blue = ( unsigned char )b;
394  return true;
395 }
396 
397 bool QgsRelief::setElevationColor( double elevation, int* red, int* green, int* blue )
398 {
399  QList< ReliefColor >::const_iterator reliefColorIt = mReliefColors.constBegin();
400  for ( ; reliefColorIt != mReliefColors.constEnd(); ++reliefColorIt )
401  {
402  if ( elevation >= reliefColorIt->minElevation && elevation <= reliefColorIt->maxElevation )
403  {
404  const QColor& c = reliefColorIt->color;
405  *red = c.red();
406  *green = c.green();
407  *blue = c.blue();
408 
409  return true;
410  }
411  }
412  return false;
413 }
414 
415 //duplicated from QgsNineCellFilter. Todo: make common base class
416 GDALDatasetH QgsRelief::openInputFile( int& nCellsX, int& nCellsY )
417 {
418  GDALDatasetH inputDataset = GDALOpen( TO8F( mInputFile ), GA_ReadOnly );
419  if ( inputDataset )
420  {
421  nCellsX = GDALGetRasterXSize( inputDataset );
422  nCellsY = GDALGetRasterYSize( inputDataset );
423 
424  //we need at least one band
425  if ( GDALGetRasterCount( inputDataset ) < 1 )
426  {
427  GDALClose( inputDataset );
428  return nullptr;
429  }
430  }
431  return inputDataset;
432 }
433 
434 GDALDriverH QgsRelief::openOutputDriver()
435 {
436  char **driverMetadata;
437 
438  //open driver
439  GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
440 
441  if ( !outputDriver )
442  {
443  return outputDriver; //return nullptr, driver does not exist
444  }
445 
446  driverMetadata = GDALGetMetadata( outputDriver, nullptr );
447  if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
448  {
449  return nullptr; //driver exist, but it does not support the create operation
450  }
451 
452  return outputDriver;
453 }
454 
455 GDALDatasetH QgsRelief::openOutputFile( GDALDatasetH inputDataset, GDALDriverH outputDriver )
456 {
457  if ( !inputDataset )
458  {
459  return nullptr;
460  }
461 
462  int xSize = GDALGetRasterXSize( inputDataset );
463  int ySize = GDALGetRasterYSize( inputDataset );
464 
465  //open output file
466  char **papszOptions = nullptr;
467 
468  //use PACKBITS compression for tiffs by default
469  papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "PACKBITS" );
470 
471  //create three band raster (reg, green, blue)
472  GDALDatasetH outputDataset = GDALCreate( outputDriver, TO8F( mOutputFile ), xSize, ySize, 3, GDT_Byte, papszOptions );
473  if ( !outputDataset )
474  {
475  return outputDataset;
476  }
477 
478  //get geotransform from inputDataset
479  double geotransform[6];
480  if ( GDALGetGeoTransform( inputDataset, geotransform ) != CE_None )
481  {
482  GDALClose( outputDataset );
483  return nullptr;
484  }
485  GDALSetGeoTransform( outputDataset, geotransform );
486 
487  //make sure mCellSizeX and mCellSizeY are always > 0
488  mCellSizeX = geotransform[1];
489  if ( mCellSizeX < 0 )
490  {
491  mCellSizeX = -mCellSizeX;
492  }
493  mCellSizeY = geotransform[5];
494  if ( mCellSizeY < 0 )
495  {
496  mCellSizeY = -mCellSizeY;
497  }
498 
499  const char* projection = GDALGetProjectionRef( inputDataset );
500  GDALSetProjection( outputDataset, projection );
501 
502  return outputDataset;
503 }
504 
505 //this function is mainly there for debugging
507 {
508  int nCellsX, nCellsY;
509  GDALDatasetH inputDataset = openInputFile( nCellsX, nCellsY );
510  if ( !inputDataset )
511  {
512  return false;
513  }
514 
515  //open first raster band for reading (elevation raster is always single band)
516  GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset, 1 );
517  if ( !elevationBand )
518  {
519  GDALClose( inputDataset );
520  return false;
521  }
522 
523  //1. get minimum and maximum of elevation raster -> 252 elevation classes
524  int minOk, maxOk;
525  double minMax[2];
526  minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
527  minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
528 
529  if ( !minOk || !maxOk )
530  {
531  GDALComputeRasterMinMax( elevationBand, true, minMax );
532  }
533 
534  //2. go through raster cells and get frequency of classes
535 
536  //store elevation frequency in 256 elevation classes
537  double frequency[252];
538  double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
539  //initialize to zero
540  for ( int i = 0; i < 252; ++i )
541  {
542  frequency[i] = 0;
543  }
544 
545  float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
546  int elevationClass = -1;
547 
548  for ( int i = 0; i < nCellsY; ++i )
549  {
550  if ( GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1,
551  scanLine, nCellsX, 1, GDT_Float32,
552  0, 0 ) != CE_None )
553  {
554  QgsDebugMsg( "Raster IO Error" );
555  }
556 
557  for ( int j = 0; j < nCellsX; ++j )
558  {
559  elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange );
560  if ( elevationClass >= 0 )
561  {
562  frequency[elevationClass] += 1.0;
563  }
564  }
565  }
566 
567  CPLFree( scanLine );
568 
569  //log10 transformation for all frequency values
570  for ( int i = 0; i < 252; ++i )
571  {
572  frequency[i] = log10( frequency[i] );
573  }
574 
575  //write out frequency values to csv file for debugging
576  QFile outFile( file );
577  if ( !outFile.open( QIODevice::WriteOnly ) )
578  {
579  return false;
580  }
581 
582  QTextStream outstream( &outFile );
583  for ( int i = 0; i < 252; ++i )
584  {
585  outstream << QString::number( i ) + ',' + QString::number( frequency[i] ) << endl;
586  }
587  outFile.close();
588  return true;
589 }
590 
592 {
594 
595  int nCellsX, nCellsY;
596  GDALDatasetH inputDataset = openInputFile( nCellsX, nCellsY );
597  if ( !inputDataset )
598  {
599  return resultList;
600  }
601 
602  //open first raster band for reading (elevation raster is always single band)
603  GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset, 1 );
604  if ( !elevationBand )
605  {
606  GDALClose( inputDataset );
607  return resultList;
608  }
609 
610  //1. get minimum and maximum of elevation raster -> 252 elevation classes
611  int minOk, maxOk;
612  double minMax[2];
613  minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
614  minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
615 
616  if ( !minOk || !maxOk )
617  {
618  GDALComputeRasterMinMax( elevationBand, true, minMax );
619  }
620 
621  //2. go through raster cells and get frequency of classes
622 
623  //store elevation frequency in 256 elevation classes
624  double frequency[252];
625  double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
626  //initialize to zero
627  for ( int i = 0; i < 252; ++i )
628  {
629  frequency[i] = 0;
630  }
631 
632  float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
633  int elevationClass = -1;
634 
635  for ( int i = 0; i < nCellsY; ++i )
636  {
637  if ( GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1,
638  scanLine, nCellsX, 1, GDT_Float32,
639  0, 0 ) != CE_None )
640  {
641  QgsDebugMsg( "Raster IO Error" );
642  }
643  for ( int j = 0; j < nCellsX; ++j )
644  {
645  elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange );
646  if ( elevationClass < 0 )
647  {
648  elevationClass = 0;
649  }
650  else if ( elevationClass >= 252 )
651  {
652  elevationClass = 251;
653  }
654  frequency[elevationClass] += 1.0;
655  }
656  }
657 
658  CPLFree( scanLine );
659 
660  //log10 transformation for all frequency values
661  for ( int i = 0; i < 252; ++i )
662  {
663  frequency[i] = log10( frequency[i] );
664  }
665 
666  //start with 9 uniformly distributed classes
667  QList<int> classBreaks;
668  classBreaks.append( 0 );
669  classBreaks.append( 28 );
670  classBreaks.append( 56 );
671  classBreaks.append( 84 );
672  classBreaks.append( 112 );
673  classBreaks.append( 140 );
674  classBreaks.append( 168 );
675  classBreaks.append( 196 );
676  classBreaks.append( 224 );
677  classBreaks.append( 252 );
678 
679  for ( int i = 0; i < 10; ++i )
680  {
681  optimiseClassBreaks( classBreaks, frequency );
682  }
683 
684  //debug, print out all the classbreaks
685  for ( int i = 0; i < classBreaks.size(); ++i )
686  {
687  qWarning( "%d", classBreaks[i] );
688  }
689 
690  //set colors according to optimised class breaks
691  QVector<QColor> colorList;
692  colorList.push_back( QColor( 7, 165, 144 ) );
693  colorList.push_back( QColor( 12, 221, 162 ) );
694  colorList.push_back( QColor( 33, 252, 183 ) );
695  colorList.push_back( QColor( 247, 252, 152 ) );
696  colorList.push_back( QColor( 252, 196, 8 ) );
697  colorList.push_back( QColor( 252, 166, 15 ) );
698  colorList.push_back( QColor( 175, 101, 15 ) );
699  colorList.push_back( QColor( 255, 133, 92 ) );
700  colorList.push_back( QColor( 204, 204, 204 ) );
701 
702  resultList.reserve( classBreaks.size() );
703  for ( int i = 1; i < classBreaks.size(); ++i )
704  {
705  double minElevation = minMax[0] + classBreaks[i - 1] * frequencyClassRange;
706  double maxElevation = minMax[0] + classBreaks[i] * frequencyClassRange;
707  resultList.push_back( QgsRelief::ReliefColor( colorList.at( i - 1 ), minElevation, maxElevation ) );
708  }
709 
710  return resultList;
711 }
712 
713 void QgsRelief::optimiseClassBreaks( QList<int>& breaks, double* frequencies )
714 {
715  int nClasses = breaks.size() - 1;
716  double* a = new double[nClasses]; //slopes
717  double* b = new double[nClasses]; //y-offsets
718 
719  for ( int i = 0; i < nClasses; ++i )
720  {
721  //get all the values between the class breaks into input
722  QList< QPair < int, double > > regressionInput;
723  for ( int j = breaks.at( i ); j < breaks.at( i + 1 ); ++j )
724  {
725  regressionInput.push_back( qMakePair( j, frequencies[j] ) );
726  }
727 
728  double aParam, bParam;
729  if ( !regressionInput.isEmpty() && calculateRegression( regressionInput, aParam, bParam ) )
730  {
731  a[i] = aParam;
732  b[i] = bParam;
733  }
734  else
735  {
736  a[i] = 0;
737  b[i] = 0; //better default value
738  }
739  }
740 
741  QList<int> classesToRemove;
742 
743  //shift class boundaries or eliminate classes which fall together
744  for ( int i = 1; i < nClasses ; ++i )
745  {
746  if ( breaks[i] == breaks[ i - 1 ] )
747  {
748  continue;
749  }
750 
751  if ( qgsDoubleNear( a[i - 1 ], a[i] ) )
752  {
753  continue;
754  }
755  else
756  {
757  int newX = ( b[i - 1] - b[ i ] ) / ( a[ i ] - a[ i - 1 ] );
758 
759  if ( newX <= breaks[i - 1] )
760  {
761  newX = breaks[i - 1];
762  // classesToRemove.push_back( i );//remove this class later as it falls together with the preceding one
763  }
764  else if ( i < nClasses - 1 && newX >= breaks[i + 1] )
765  {
766  newX = breaks[i + 1];
767  // classesToRemove.push_back( i );//remove this class later as it falls together with the next one
768  }
769 
770  breaks[i] = newX;
771  }
772  }
773 
774  for ( int i = classesToRemove.size() - 1; i >= 0; --i )
775  {
776  breaks.removeAt( classesToRemove.at( i ) );
777  }
778 
779  delete[] a;
780  delete[] b;
781 }
782 
783 int QgsRelief::frequencyClassForElevation( double elevation, double minElevation, double elevationClassRange )
784 {
785  return ( elevation - minElevation ) / elevationClassRange;
786 }
787 
788 bool QgsRelief::calculateRegression( const QList< QPair < int, double > >& input, double& a, double& b )
789 {
790  double xMean, yMean;
791  double xSum = 0;
792  double ySum = 0;
793  QList< QPair < int, double > >::const_iterator inputIt = input.constBegin();
794  for ( ; inputIt != input.constEnd(); ++inputIt )
795  {
796  xSum += inputIt->first;
797  ySum += inputIt->second;
798  }
799  xMean = xSum / input.size();
800  yMean = ySum / input.size();
801 
802  double sumCounter = 0;
803  double sumDenominator = 0;
804  inputIt = input.constBegin();
805  for ( ; inputIt != input.constEnd(); ++inputIt )
806  {
807  sumCounter += (( inputIt->first - xMean ) * ( inputIt->second - yMean ) );
808  sumDenominator += (( inputIt->first - xMean ) * ( inputIt->first - xMean ) );
809  }
810 
811  a = sumCounter / sumDenominator;
812  b = yMean - a * xMean;
813 
814  return true;
815 }
816 
float processNineCellWindow(float *x11, float *x21, float *x31, float *x12, float *x22, float *x32, float *x13, float *x23, float *x33) override
Calculates output value from nine input values.
void setZFactor(double factor)
void addReliefColorClass(const ReliefColor &color)
Definition: qgsrelief.cpp:71
int processRaster(QProgressDialog *p)
Starts the calculation, reads from mInputFile and stores the result in mOutputFile.
Definition: qgsrelief.cpp:87
void setMaximum(int maximum)
void push_back(const T &value)
QList< ReliefColor > calculateOptimizedReliefClasses()
Calculates class breaks according with the method of Buenzli (2011) using an iterative algorithm for ...
Definition: qgsrelief.cpp:591
void setCellSizeY(double size)
#define QgsDebugMsg(str)
Definition: qgslogger.h:33
void reserve(int alloc)
const T & at(int i) const
void removeAt(int i)
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Compare two doubles (but allow some difference)
Definition: qgis.h:353
int size() const
#define TO8F(x)
Definition: qgsrelief.cpp:34
void setValue(int progress)
bool exportFrequencyDistributionToCsv(const QString &file)
Write frequency of elevation values to file for manual inspection.
Definition: qgsrelief.cpp:506
QString number(int n, int base)
void append(const T &value)
float processNineCellWindow(float *x11, float *x21, float *x31, float *x12, float *x22, float *x32, float *x13, float *x23, float *x33) override
Calculates output value from nine input values.
int red() const
void setInputNodataValue(double value)
Calculates aspect values in a window of 3x3 cells based on first order derivatives in x- and y- direc...
bool isEmpty() const
#define M_PI
void setCellSizeX(double size)
T & first()
void * GDALDatasetH
virtual bool open(QFlags< QIODevice::OpenModeFlag > mode)
int green() const
QByteArray toLocal8Bit() const
void clearReliefColors()
Definition: qgsrelief.cpp:66
virtual void close()
int blue() const
const T & at(int i) const
float processNineCellWindow(float *x11, float *x21, float *x31, float *x12, float *x22, float *x32, float *x13, float *x23, float *x33) override
Calculates output value from nine input values.
Calculates slope values in a window of 3x3 cells based on first order derivatives in x- and y- direct...
void setOutputNodataValue(double value)
char * data()
void push_back(const T &value)
QgsRelief(const QString &inputFile, const QString &outputFile, const QString &outputFormat)
Definition: qgsrelief.cpp:37
const_iterator constEnd() const
const_iterator constBegin() const