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