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