QGIS API Documentation 3.41.0-Master (d5b93354e9c)
Loading...
Searching...
No Matches
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
35QgsRelief::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
49QgsRelief::~QgsRelief() = default;
50
52{
53 mReliefColors.clear();
54}
55
57{
58 mReliefColors.push_back( color );
59}
60
61void 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], &scanLine2[j + 1], &mInputNodataValue, &scanLine3[j], &scanLine3[j + 1], &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] );
214 }
215 else if ( j == xSize - 1 )
216 {
217 resultOk = processNineCellWindow( &scanLine1[j - 1], &scanLine1[j], &mInputNodataValue, &scanLine2[j - 1], &scanLine2[j], &mInputNodataValue, &scanLine3[j - 1], &scanLine3[j], &mInputNodataValue, &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] );
218 }
219 else
220 {
221 resultOk = processNineCellWindow( &scanLine1[j - 1], &scanLine1[j], &scanLine1[j + 1], &scanLine2[j - 1], &scanLine2[j], &scanLine2[j + 1], &scanLine3[j - 1], &scanLine3[j], &scanLine3[j + 1], &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] );
222 }
223
224 if ( !resultOk )
225 {
226 resultRedLine[j] = mOutputNodataValue;
227 resultGreenLine[j] = mOutputNodataValue;
228 resultBlueLine[j] = mOutputNodataValue;
229 }
230 }
231
232 if ( GDALRasterIO( outputRedBand, GF_Write, 0, i, xSize, 1, resultRedLine, xSize, 1, GDT_Byte, 0, 0 ) != CE_None )
233 {
234 QgsDebugError( QStringLiteral( "Raster IO Error" ) );
235 }
236 if ( GDALRasterIO( outputGreenBand, GF_Write, 0, i, xSize, 1, resultGreenLine, xSize, 1, GDT_Byte, 0, 0 ) != CE_None )
237 {
238 QgsDebugError( QStringLiteral( "Raster IO Error" ) );
239 }
240 if ( GDALRasterIO( outputBlueBand, GF_Write, 0, i, xSize, 1, resultBlueLine, xSize, 1, GDT_Byte, 0, 0 ) != CE_None )
241 {
242 QgsDebugError( QStringLiteral( "Raster IO Error" ) );
243 }
244 }
245
246 if ( feedback )
247 {
248 feedback->setProgress( 100 );
249 }
250
251 CPLFree( resultRedLine );
252 CPLFree( resultBlueLine );
253 CPLFree( resultGreenLine );
254 CPLFree( scanLine1 );
255 CPLFree( scanLine2 );
256 CPLFree( scanLine3 );
257
258 if ( feedback && feedback->isCanceled() )
259 {
260 //delete the dataset without closing (because it is faster)
261 gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
262 return 7;
263 }
264
265 return 0;
266}
267
268bool QgsRelief::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 )
269{
270 //1. component: color and hillshade from 300 degrees
271 int r = 0;
272 int g = 0;
273 int b = 0;
274
275 const float hillShadeValue300 = mHillshadeFilter300->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
276 if ( hillShadeValue300 != mOutputNodataValue )
277 {
278 if ( !getElevationColor( *x5, &r, &g, &b ) )
279 {
280 r = hillShadeValue300;
281 g = hillShadeValue300;
282 b = hillShadeValue300;
283 }
284 else
285 {
286 r = r / 2.0 + hillShadeValue300 / 2.0;
287 g = g / 2.0 + hillShadeValue300 / 2.0;
288 b = b / 2.0 + hillShadeValue300 / 2.0;
289 }
290 }
291
292 //2. component: hillshade and slope
293 const float hillShadeValue315 = mHillshadeFilter315->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
294 const float slope = mSlopeFilter->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
295 if ( hillShadeValue315 != mOutputNodataValue && slope != mOutputNodataValue )
296 {
297 int r2, g2, b2;
298 if ( slope > 15 )
299 {
300 r2 = 0 / 2.0 + hillShadeValue315 / 2.0;
301 g2 = 0 / 2.0 + hillShadeValue315 / 2.0;
302 b2 = 0 / 2.0 + hillShadeValue315 / 2.0;
303 }
304 else if ( slope >= 1 )
305 {
306 const int slopeValue = 255 - ( slope / 15.0 * 255.0 );
307 r2 = slopeValue / 2.0 + hillShadeValue315 / 2.0;
308 g2 = slopeValue / 2.0 + hillShadeValue315 / 2.0;
309 b2 = slopeValue / 2.0 + hillShadeValue315 / 2.0;
310 }
311 else
312 {
313 r2 = hillShadeValue315;
314 g2 = hillShadeValue315;
315 b2 = hillShadeValue315;
316 }
317
318 //combine with r,g,b with 70 percentage coverage
319 r = r * 0.7 + r2 * 0.3;
320 g = g * 0.7 + g2 * 0.3;
321 b = b * 0.7 + b2 * 0.3;
322 }
323
324 //3. combine yellow aspect with 10% transparency, illumination from 285 degrees
325 const float hillShadeValue285 = mHillshadeFilter285->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
326 const float aspect = mAspectFilter->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
327 if ( hillShadeValue285 != mOutputNodataValue && aspect != mOutputNodataValue )
328 {
329 double angle_diff = std::fabs( 285 - aspect );
330 if ( angle_diff > 180 )
331 {
332 angle_diff -= 180;
333 }
334
335 int r3, g3, b3;
336 if ( angle_diff < 90 )
337 {
338 const int aspectVal = ( 1 - std::cos( angle_diff * M_PI / 180 ) ) * 255;
339 r3 = 0.5 * 255 + hillShadeValue315 * 0.5;
340 g3 = 0.5 * 255 + hillShadeValue315 * 0.5;
341 b3 = 0.5 * aspectVal + hillShadeValue315 * 0.5;
342 }
343 else //white
344 {
345 r3 = 0.5 * 255 + hillShadeValue315 * 0.5;
346 g3 = 0.5 * 255 + hillShadeValue315 * 0.5;
347 b3 = 0.5 * 255 + hillShadeValue315 * 0.5;
348 }
349
350 r = r3 * 0.1 + r * 0.9;
351 g = g3 * 0.1 + g * 0.9;
352 b = b3 * 0.1 + b * 0.9;
353 }
354
355 *red = ( unsigned char ) r;
356 *green = ( unsigned char ) g;
357 *blue = ( unsigned char ) b;
358 return true;
359}
360
361bool QgsRelief::getElevationColor( double elevation, int *red, int *green, int *blue ) const
362{
363 QList<ReliefColor>::const_iterator reliefColorIt = mReliefColors.constBegin();
364 for ( ; reliefColorIt != mReliefColors.constEnd(); ++reliefColorIt )
365 {
366 if ( elevation >= reliefColorIt->minElevation && elevation <= reliefColorIt->maxElevation )
367 {
368 const QColor &c = reliefColorIt->color;
369 *red = c.red();
370 *green = c.green();
371 *blue = c.blue();
372
373 return true;
374 }
375 }
376 return false;
377}
378
379//duplicated from QgsNineCellFilter. Todo: make common base class
380gdal::dataset_unique_ptr QgsRelief::openInputFile( int &nCellsX, int &nCellsY )
381{
382 gdal::dataset_unique_ptr inputDataset( GDALOpen( mInputFile.toUtf8().constData(), GA_ReadOnly ) );
383 if ( inputDataset )
384 {
385 nCellsX = GDALGetRasterXSize( inputDataset.get() );
386 nCellsY = GDALGetRasterYSize( inputDataset.get() );
387
388 //we need at least one band
389 if ( GDALGetRasterCount( inputDataset.get() ) < 1 )
390 {
391 return nullptr;
392 }
393 }
394 return inputDataset;
395}
396
397GDALDriverH QgsRelief::openOutputDriver()
398{
399 //open driver
400 GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
401
402 if ( !outputDriver )
403 {
404 return outputDriver; //return nullptr, driver does not exist
405 }
406
407 if ( !QgsGdalUtils::supportsRasterCreate( outputDriver ) )
408 {
409 return nullptr; //driver exist, but it does not support the create operation
410 }
411
412 return outputDriver;
413}
414
415gdal::dataset_unique_ptr QgsRelief::openOutputFile( GDALDatasetH inputDataset, GDALDriverH outputDriver )
416{
417 if ( !inputDataset )
418 {
419 return nullptr;
420 }
421
422 const int xSize = GDALGetRasterXSize( inputDataset );
423 const int ySize = GDALGetRasterYSize( inputDataset );
424
425 //open output file
426 char **papszOptions = nullptr;
427
428 //use PACKBITS compression for tiffs by default
429 papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "PACKBITS" );
430
431 //create three band raster (red, green, blue)
432 gdal::dataset_unique_ptr outputDataset( GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), xSize, ySize, 3, GDT_Byte, papszOptions ) );
433 CSLDestroy( papszOptions );
434 papszOptions = nullptr;
435
436 if ( !outputDataset )
437 {
438 return nullptr;
439 }
440
441 //get geotransform from inputDataset
442 double geotransform[6];
443 if ( GDALGetGeoTransform( inputDataset, geotransform ) != CE_None )
444 {
445 return nullptr;
446 }
447 GDALSetGeoTransform( outputDataset.get(), geotransform );
448
449 //make sure mCellSizeX and mCellSizeY are always > 0
450 mCellSizeX = geotransform[1];
451 if ( mCellSizeX < 0 )
452 {
453 mCellSizeX = -mCellSizeX;
454 }
455 mCellSizeY = geotransform[5];
456 if ( mCellSizeY < 0 )
457 {
458 mCellSizeY = -mCellSizeY;
459 }
460
461 const char *projection = GDALGetProjectionRef( inputDataset );
462 GDALSetProjection( outputDataset.get(), projection );
463
464 return outputDataset;
465}
466
467//this function is mainly there for debugging
469{
470 int nCellsX, nCellsY;
471 const gdal::dataset_unique_ptr inputDataset = openInputFile( nCellsX, nCellsY );
472 if ( !inputDataset )
473 {
474 return false;
475 }
476
477 //open first raster band for reading (elevation raster is always single band)
478 GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset.get(), 1 );
479 if ( !elevationBand )
480 {
481 return false;
482 }
483
484 //1. get minimum and maximum of elevation raster -> 252 elevation classes
485 int minOk, maxOk;
486 double minMax[2];
487 minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
488 minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
489
490 if ( !minOk || !maxOk )
491 {
492 GDALComputeRasterMinMax( elevationBand, true, minMax );
493 }
494
495 //2. go through raster cells and get frequency of classes
496
497 //store elevation frequency in 256 elevation classes
498 double frequency[252] = { 0 };
499 const double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
500
501 float *scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
502 int elevationClass = -1;
503
504 for ( int i = 0; i < nCellsY; ++i )
505 {
506 if ( GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 ) != CE_None )
507 {
508 QgsDebugError( QStringLiteral( "Raster IO Error" ) );
509 }
510
511 for ( int j = 0; j < nCellsX; ++j )
512 {
513 elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange );
514 if ( elevationClass >= 0 && elevationClass < 252 )
515 {
516 frequency[elevationClass] += 1.0;
517 }
518 }
519 }
520
521 CPLFree( scanLine );
522
523 //log10 transformation for all frequency values
524 for ( int i = 0; i < 252; ++i )
525 {
526 frequency[i] = std::log10( frequency[i] );
527 }
528
529 //write out frequency values to csv file for debugging
530 QFile outFile( file );
531 if ( !outFile.open( QIODevice::WriteOnly | QIODevice::Truncate ) )
532 {
533 return false;
534 }
535
536 QTextStream outstream( &outFile );
537#if QT_VERSION < QT_VERSION_CHECK( 6, 0, 0 )
538 outstream.setCodec( "UTF-8" );
539#endif
540 for ( int i = 0; i < 252; ++i )
541 {
542 outstream << QString::number( i ) + ',' + QString::number( frequency[i] ) << Qt::endl;
543 }
544 outFile.close();
545 return true;
546}
547
548QList<QgsRelief::ReliefColor> QgsRelief::calculateOptimizedReliefClasses()
549{
550 QList<QgsRelief::ReliefColor> resultList;
551
552 int nCellsX, nCellsY;
553 const gdal::dataset_unique_ptr inputDataset = openInputFile( nCellsX, nCellsY );
554 if ( !inputDataset )
555 {
556 return resultList;
557 }
558
559 //open first raster band for reading (elevation raster is always single band)
560 GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset.get(), 1 );
561 if ( !elevationBand )
562 {
563 return resultList;
564 }
565
566 //1. get minimum and maximum of elevation raster -> 252 elevation classes
567 int minOk, maxOk;
568 double minMax[2];
569 minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
570 minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
571
572 if ( !minOk || !maxOk )
573 {
574 GDALComputeRasterMinMax( elevationBand, true, minMax );
575 }
576
577 //2. go through raster cells and get frequency of classes
578
579 //store elevation frequency in 256 elevation classes
580 double frequency[252] = { 0 };
581 const double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
582
583 float *scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
584 int elevationClass = -1;
585
586 for ( int i = 0; i < nCellsY; ++i )
587 {
588 if ( GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 ) != CE_None )
589 {
590 QgsDebugError( QStringLiteral( "Raster IO Error" ) );
591 }
592 for ( int j = 0; j < nCellsX; ++j )
593 {
594 elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange );
595 elevationClass = std::max( std::min( elevationClass, 251 ), 0 );
596 frequency[elevationClass] += 1.0;
597 }
598 }
599
600 CPLFree( scanLine );
601
602 //log10 transformation for all frequency values
603 for ( int i = 0; i < 252; ++i )
604 {
605 frequency[i] = std::log10( frequency[i] );
606 }
607
608 //start with 9 uniformly distributed classes
609 QList<int> classBreaks;
610 classBreaks.append( 0 );
611 classBreaks.append( 28 );
612 classBreaks.append( 56 );
613 classBreaks.append( 84 );
614 classBreaks.append( 112 );
615 classBreaks.append( 140 );
616 classBreaks.append( 168 );
617 classBreaks.append( 196 );
618 classBreaks.append( 224 );
619 classBreaks.append( 252 );
620
621 for ( int i = 0; i < 10; ++i )
622 {
623 optimiseClassBreaks( classBreaks, frequency );
624 }
625
626 //debug, print out all the classbreaks
627 for ( int i = 0; i < classBreaks.size(); ++i )
628 {
629 qWarning( "%d", classBreaks[i] );
630 }
631
632 //set colors according to optimised class breaks
633 QVector<QColor> colorList;
634 colorList.reserve( 9 );
635 colorList.push_back( QColor( 7, 165, 144 ) );
636 colorList.push_back( QColor( 12, 221, 162 ) );
637 colorList.push_back( QColor( 33, 252, 183 ) );
638 colorList.push_back( QColor( 247, 252, 152 ) );
639 colorList.push_back( QColor( 252, 196, 8 ) );
640 colorList.push_back( QColor( 252, 166, 15 ) );
641 colorList.push_back( QColor( 175, 101, 15 ) );
642 colorList.push_back( QColor( 255, 133, 92 ) );
643 colorList.push_back( QColor( 204, 204, 204 ) );
644
645 resultList.reserve( classBreaks.size() );
646 for ( int i = 1; i < classBreaks.size(); ++i )
647 {
648 const double minElevation = minMax[0] + classBreaks[i - 1] * frequencyClassRange;
649 const double maxElevation = minMax[0] + classBreaks[i] * frequencyClassRange;
650 resultList.push_back( QgsRelief::ReliefColor( colorList.at( i - 1 ), minElevation, maxElevation ) );
651 }
652
653 return resultList;
654}
655
656void QgsRelief::optimiseClassBreaks( QList<int> &breaks, double *frequencies )
657{
658 const int nClasses = breaks.size() - 1;
659 double *a = new double[nClasses]; //slopes
660 double *b = new double[nClasses]; //y-offsets
661
662 for ( int i = 0; i < nClasses; ++i )
663 {
664 //get all the values between the class breaks into input
665 QList<QPair<int, double>> regressionInput;
666 regressionInput.reserve( breaks.at( i + 1 ) - breaks.at( i ) );
667 for ( int j = breaks.at( i ); j < breaks.at( i + 1 ); ++j )
668 {
669 regressionInput.push_back( qMakePair( j, frequencies[j] ) );
670 }
671
672 double aParam, bParam;
673 if ( !regressionInput.isEmpty() && calculateRegression( regressionInput, aParam, bParam ) )
674 {
675 a[i] = aParam;
676 b[i] = bParam;
677 }
678 else
679 {
680 a[i] = 0;
681 b[i] = 0; //better default value
682 }
683 }
684
685 const QList<int> classesToRemove;
686
687 //shift class boundaries or eliminate classes which fall together
688 for ( int i = 1; i < nClasses; ++i )
689 {
690 if ( breaks[i] == breaks[i - 1] )
691 {
692 continue;
693 }
694
695 if ( qgsDoubleNear( a[i - 1], a[i] ) )
696 {
697 continue;
698 }
699 else
700 {
701 int newX = ( b[i - 1] - b[i] ) / ( a[i] - a[i - 1] );
702
703 if ( newX <= breaks[i - 1] )
704 {
705 newX = breaks[i - 1];
706 // classesToRemove.push_back( i );//remove this class later as it falls together with the preceding one
707 }
708 else if ( i < nClasses - 1 && newX >= breaks[i + 1] )
709 {
710 newX = breaks[i + 1];
711 // classesToRemove.push_back( i );//remove this class later as it falls together with the next one
712 }
713
714 breaks[i] = newX;
715 }
716 }
717
718 for ( int i = classesToRemove.size() - 1; i >= 0; --i )
719 {
720 breaks.removeAt( classesToRemove.at( i ) ); // cppcheck-suppress containerOutOfBounds
721 }
722
723 delete[] a;
724 delete[] b;
725}
726
727int QgsRelief::frequencyClassForElevation( double elevation, double minElevation, double elevationClassRange )
728{
729 return ( elevation - minElevation ) / elevationClassRange;
730}
731
732bool QgsRelief::calculateRegression( const QList<QPair<int, double>> &input, double &a, double &b )
733{
734 double xMean, yMean;
735 double xSum = 0;
736 double ySum = 0;
737 QList<QPair<int, double>>::const_iterator inputIt = input.constBegin();
738 for ( ; inputIt != input.constEnd(); ++inputIt )
739 {
740 xSum += inputIt->first;
741 ySum += inputIt->second;
742 }
743 xMean = xSum / input.size();
744 yMean = ySum / input.size();
745
746 double sumCounter = 0;
747 double sumDenominator = 0;
748 inputIt = input.constBegin();
749 for ( ; inputIt != input.constEnd(); ++inputIt )
750 {
751 sumCounter += ( ( inputIt->first - xMean ) * ( inputIt->second - yMean ) );
752 sumDenominator += ( ( inputIt->first - xMean ) * ( inputIt->first - xMean ) );
753 }
754
755 a = sumCounter / sumDenominator;
756 b = yMean - a * xMean;
757
758 return true;
759}
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.
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 ...
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.
std::unique_ptr< std::remove_pointer< GDALDatasetH >::type, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
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:6024
void * GDALDatasetH
#define QgsDebugError(str)
Definition qgslogger.h:38