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
QgsRelief(const QString & inputFile,const QString & outputFile,const QString & outputFormat)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
clearReliefColors()51 void QgsRelief::clearReliefColors()
52 {
53 mReliefColors.clear();
54 }
55
addReliefColorClass(const ReliefColor & color)56 void QgsRelief::addReliefColorClass( const ReliefColor &color )
57 {
58 mReliefColors.push_back( color );
59 }
60
setDefaultReliefColors()61 void QgsRelief::setDefaultReliefColors()
62 {
63 clearReliefColors();
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
processRaster(QgsFeedback * feedback)72 int QgsRelief::processRaster( QgsFeedback *feedback )
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 QgsDebugMsg( 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 QgsDebugMsg( 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 QgsDebugMsg( 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 QgsDebugMsg( 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 QgsDebugMsg( 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
processNineCellWindow(float * x1,float * x2,float * x3,float * x4,float * x5,float * x6,float * x7,float * x8,float * x9,unsigned char * red,unsigned char * green,unsigned char * blue)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
setElevationColor(double elevation,int * red,int * green,int * blue)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
openInputFile(int & nCellsX,int & nCellsY)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
openOutputDriver()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
openOutputFile(GDALDatasetH inputDataset,GDALDriverH outputDriver)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
exportFrequencyDistributionToCsv(const QString & file)475 bool QgsRelief::exportFrequencyDistributionToCsv( const QString &file )
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 QgsDebugMsg( 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 #if QT_VERSION < QT_VERSION_CHECK(5, 14, 0)
552 outstream << QString::number( i ) + ',' + QString::number( frequency[i] ) << endl;
553 #else
554 outstream << QString::number( i ) + ',' + QString::number( frequency[i] ) << Qt::endl;
555 #endif
556 }
557 outFile.close();
558 return true;
559 }
560
calculateOptimizedReliefClasses()561 QList< QgsRelief::ReliefColor > QgsRelief::calculateOptimizedReliefClasses()
562 {
563 QList< QgsRelief::ReliefColor > resultList;
564
565 int nCellsX, nCellsY;
566 const gdal::dataset_unique_ptr inputDataset = openInputFile( nCellsX, nCellsY );
567 if ( !inputDataset )
568 {
569 return resultList;
570 }
571
572 //open first raster band for reading (elevation raster is always single band)
573 GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset.get(), 1 );
574 if ( !elevationBand )
575 {
576 return resultList;
577 }
578
579 //1. get minimum and maximum of elevation raster -> 252 elevation classes
580 int minOk, maxOk;
581 double minMax[2];
582 minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
583 minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
584
585 if ( !minOk || !maxOk )
586 {
587 GDALComputeRasterMinMax( elevationBand, true, minMax );
588 }
589
590 //2. go through raster cells and get frequency of classes
591
592 //store elevation frequency in 256 elevation classes
593 double frequency[252] = {0};
594 const double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
595
596 float *scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
597 int elevationClass = -1;
598
599 for ( int i = 0; i < nCellsY; ++i )
600 {
601 if ( GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1,
602 scanLine, nCellsX, 1, GDT_Float32,
603 0, 0 ) != CE_None )
604 {
605 QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
606 }
607 for ( int j = 0; j < nCellsX; ++j )
608 {
609 elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange );
610 elevationClass = std::max( std::min( elevationClass, 251 ), 0 );
611 frequency[elevationClass] += 1.0;
612 }
613 }
614
615 CPLFree( scanLine );
616
617 //log10 transformation for all frequency values
618 for ( int i = 0; i < 252; ++i )
619 {
620 frequency[i] = std::log10( frequency[i] );
621 }
622
623 //start with 9 uniformly distributed classes
624 QList<int> classBreaks;
625 classBreaks.append( 0 );
626 classBreaks.append( 28 );
627 classBreaks.append( 56 );
628 classBreaks.append( 84 );
629 classBreaks.append( 112 );
630 classBreaks.append( 140 );
631 classBreaks.append( 168 );
632 classBreaks.append( 196 );
633 classBreaks.append( 224 );
634 classBreaks.append( 252 );
635
636 for ( int i = 0; i < 10; ++i )
637 {
638 optimiseClassBreaks( classBreaks, frequency );
639 }
640
641 //debug, print out all the classbreaks
642 for ( int i = 0; i < classBreaks.size(); ++i )
643 {
644 qWarning( "%d", classBreaks[i] );
645 }
646
647 //set colors according to optimised class breaks
648 QVector<QColor> colorList;
649 colorList.reserve( 9 );
650 colorList.push_back( QColor( 7, 165, 144 ) );
651 colorList.push_back( QColor( 12, 221, 162 ) );
652 colorList.push_back( QColor( 33, 252, 183 ) );
653 colorList.push_back( QColor( 247, 252, 152 ) );
654 colorList.push_back( QColor( 252, 196, 8 ) );
655 colorList.push_back( QColor( 252, 166, 15 ) );
656 colorList.push_back( QColor( 175, 101, 15 ) );
657 colorList.push_back( QColor( 255, 133, 92 ) );
658 colorList.push_back( QColor( 204, 204, 204 ) );
659
660 resultList.reserve( classBreaks.size() );
661 for ( int i = 1; i < classBreaks.size(); ++i )
662 {
663 const double minElevation = minMax[0] + classBreaks[i - 1] * frequencyClassRange;
664 const double maxElevation = minMax[0] + classBreaks[i] * frequencyClassRange;
665 resultList.push_back( QgsRelief::ReliefColor( colorList.at( i - 1 ), minElevation, maxElevation ) );
666 }
667
668 return resultList;
669 }
670
optimiseClassBreaks(QList<int> & breaks,double * frequencies)671 void QgsRelief::optimiseClassBreaks( QList<int> &breaks, double *frequencies )
672 {
673 const int nClasses = breaks.size() - 1;
674 double *a = new double[nClasses]; //slopes
675 double *b = new double[nClasses]; //y-offsets
676
677 for ( int i = 0; i < nClasses; ++i )
678 {
679 //get all the values between the class breaks into input
680 QList< QPair < int, double > > regressionInput;
681 regressionInput.reserve( breaks.at( i + 1 ) - breaks.at( i ) );
682 for ( int j = breaks.at( i ); j < breaks.at( i + 1 ); ++j )
683 {
684 regressionInput.push_back( qMakePair( j, frequencies[j] ) );
685 }
686
687 double aParam, bParam;
688 if ( !regressionInput.isEmpty() && calculateRegression( regressionInput, aParam, bParam ) )
689 {
690 a[i] = aParam;
691 b[i] = bParam;
692 }
693 else
694 {
695 a[i] = 0;
696 b[i] = 0; //better default value
697 }
698 }
699
700 const QList<int> classesToRemove;
701
702 //shift class boundaries or eliminate classes which fall together
703 for ( int i = 1; i < nClasses ; ++i )
704 {
705 if ( breaks[i] == breaks[ i - 1 ] )
706 {
707 continue;
708 }
709
710 if ( qgsDoubleNear( a[i - 1 ], a[i] ) )
711 {
712 continue;
713 }
714 else
715 {
716 int newX = ( b[i - 1] - b[ i ] ) / ( a[ i ] - a[ i - 1 ] );
717
718 if ( newX <= breaks[i - 1] )
719 {
720 newX = breaks[i - 1];
721 // classesToRemove.push_back( i );//remove this class later as it falls together with the preceding one
722 }
723 else if ( i < nClasses - 1 && newX >= breaks[i + 1] )
724 {
725 newX = breaks[i + 1];
726 // classesToRemove.push_back( i );//remove this class later as it falls together with the next one
727 }
728
729 breaks[i] = newX;
730 }
731 }
732
733 for ( int i = classesToRemove.size() - 1; i >= 0; --i )
734 {
735 breaks.removeAt( classesToRemove.at( i ) );
736 }
737
738 delete[] a;
739 delete[] b;
740 }
741
frequencyClassForElevation(double elevation,double minElevation,double elevationClassRange)742 int QgsRelief::frequencyClassForElevation( double elevation, double minElevation, double elevationClassRange )
743 {
744 return ( elevation - minElevation ) / elevationClassRange;
745 }
746
calculateRegression(const QList<QPair<int,double>> & input,double & a,double & b)747 bool QgsRelief::calculateRegression( const QList< QPair < int, double > > &input, double &a, double &b )
748 {
749 double xMean, yMean;
750 double xSum = 0;
751 double ySum = 0;
752 QList< QPair < int, double > >::const_iterator inputIt = input.constBegin();
753 for ( ; inputIt != input.constEnd(); ++inputIt )
754 {
755 xSum += inputIt->first;
756 ySum += inputIt->second;
757 }
758 xMean = xSum / input.size();
759 yMean = ySum / input.size();
760
761 double sumCounter = 0;
762 double sumDenominator = 0;
763 inputIt = input.constBegin();
764 for ( ; inputIt != input.constEnd(); ++inputIt )
765 {
766 sumCounter += ( ( inputIt->first - xMean ) * ( inputIt->second - yMean ) );
767 sumDenominator += ( ( inputIt->first - xMean ) * ( inputIt->first - xMean ) );
768 }
769
770 a = sumCounter / sumDenominator;
771 b = yMean - a * xMean;
772
773 return true;
774 }
775
776