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