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