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