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