26#include "cpl_string.h"
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 ) )
53 mReliefColors.clear();
58 mReliefColors.push_back( color );
61void QgsRelief::setDefaultReliefColors()
83 GDALDriverH outputDriver = openOutputDriver();
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 );
113 GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset.get(), 1 );
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 );
125 GDALRasterBandH outputRedBand = GDALGetRasterBand( outputDataset.get(), 1 );
126 GDALRasterBandH outputGreenBand = GDALGetRasterBand( outputDataset.get(), 2 );
127 GDALRasterBandH outputBlueBand = GDALGetRasterBand( outputDataset.get(), 3 );
129 if ( !outputRedBand || !outputGreenBand || !outputBlueBand )
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 );
150 float *scanLine1 = (
float * ) CPLMalloc(
sizeof(
float ) * xSize );
151 float *scanLine2 = (
float * ) CPLMalloc(
sizeof(
float ) * xSize );
152 float *scanLine3 = (
float * ) CPLMalloc(
sizeof(
float ) * xSize );
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 );
161 for (
int i = 0; i < ySize; ++i )
165 feedback->
setProgress( 100.0 * i /
static_cast< double >( ySize ) );
176 for (
int a = 0; a < xSize; ++a )
178 scanLine1[a] = mInputNodataValue;
180 if ( GDALRasterIO( rasterBand, GF_Read, 0, 0, xSize, 1, scanLine2, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
188 CPLFree( scanLine1 );
189 scanLine1 = scanLine2;
190 scanLine2 = scanLine3;
191 scanLine3 = (
float * ) CPLMalloc(
sizeof(
float ) * xSize );
194 if ( i == ySize - 1 )
196 for (
int a = 0; a < xSize; ++a )
198 scanLine3[a] = mInputNodataValue;
203 if ( GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, scanLine3, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
209 for (
int j = 0; j < xSize; ++j )
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] );
217 else if ( j == xSize - 1 )
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] );
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] );
232 resultRedLine[j] = mOutputNodataValue;
233 resultGreenLine[j] = mOutputNodataValue;
234 resultBlueLine[j] = mOutputNodataValue;
238 if ( GDALRasterIO( outputRedBand, GF_Write, 0, i, xSize, 1, resultRedLine, xSize, 1, GDT_Byte, 0, 0 ) != CE_None )
242 if ( GDALRasterIO( outputGreenBand, GF_Write, 0, i, xSize, 1, resultGreenLine, xSize, 1, GDT_Byte, 0, 0 ) != CE_None )
246 if ( GDALRasterIO( outputBlueBand, GF_Write, 0, i, xSize, 1, resultBlueLine, xSize, 1, GDT_Byte, 0, 0 ) != CE_None )
257 CPLFree( resultRedLine );
258 CPLFree( resultBlueLine );
259 CPLFree( resultGreenLine );
260 CPLFree( scanLine1 );
261 CPLFree( scanLine2 );
262 CPLFree( scanLine3 );
274bool 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 )
282 const float hillShadeValue300 = mHillshadeFilter300->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
283 if ( hillShadeValue300 != mOutputNodataValue )
285 if ( !getElevationColor( *x5, &r, &g, &b ) )
287 r = hillShadeValue300;
288 g = hillShadeValue300;
289 b = hillShadeValue300;
293 r = r / 2.0 + hillShadeValue300 / 2.0;
294 g = g / 2.0 + hillShadeValue300 / 2.0;
295 b = b / 2.0 + hillShadeValue300 / 2.0;
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 )
307 r2 = 0 / 2.0 + hillShadeValue315 / 2.0;
308 g2 = 0 / 2.0 + hillShadeValue315 / 2.0;
309 b2 = 0 / 2.0 + hillShadeValue315 / 2.0;
311 else if ( slope >= 1 )
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;
320 r2 = hillShadeValue315;
321 g2 = hillShadeValue315;
322 b2 = hillShadeValue315;
326 r = r * 0.7 + r2 * 0.3;
327 g = g * 0.7 + g2 * 0.3;
328 b = b * 0.7 + b2 * 0.3;
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 )
336 double angle_diff = std::fabs( 285 - aspect );
337 if ( angle_diff > 180 )
343 if ( angle_diff < 90 )
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;
352 r3 = 0.5 * 255 + hillShadeValue315 * 0.5;
353 g3 = 0.5 * 255 + hillShadeValue315 * 0.5;
354 b3 = 0.5 * 255 + hillShadeValue315 * 0.5;
357 r = r3 * 0.1 + r * 0.9;
358 g = g3 * 0.1 + g * 0.9;
359 b = b3 * 0.1 + b * 0.9;
362 *red = (
unsigned char )r;
363 *green = (
unsigned char )g;
364 *blue = (
unsigned char )b;
368bool QgsRelief::getElevationColor(
double elevation,
int *red,
int *green,
int *blue )
const
370 QList< ReliefColor >::const_iterator reliefColorIt = mReliefColors.constBegin();
371 for ( ; reliefColorIt != mReliefColors.constEnd(); ++reliefColorIt )
373 if ( elevation >= reliefColorIt->minElevation && elevation <= reliefColorIt->maxElevation )
375 const QColor &
c = reliefColorIt->color;
392 nCellsX = GDALGetRasterXSize( inputDataset.get() );
393 nCellsY = GDALGetRasterYSize( inputDataset.get() );
396 if ( GDALGetRasterCount( inputDataset.get() ) < 1 )
404GDALDriverH QgsRelief::openOutputDriver()
407 GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
429 const int xSize = GDALGetRasterXSize( inputDataset );
430 const int ySize = GDALGetRasterYSize( inputDataset );
433 char **papszOptions =
nullptr;
436 papszOptions = CSLSetNameValue( papszOptions,
"COMPRESS",
"PACKBITS" );
439 gdal::dataset_unique_ptr outputDataset( GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), xSize, ySize, 3, GDT_Byte, papszOptions ) );
440 CSLDestroy( papszOptions );
441 papszOptions =
nullptr;
443 if ( !outputDataset )
449 double geotransform[6];
450 if ( GDALGetGeoTransform( inputDataset, geotransform ) != CE_None )
454 GDALSetGeoTransform( outputDataset.get(), geotransform );
457 mCellSizeX = geotransform[1];
458 if ( mCellSizeX < 0 )
460 mCellSizeX = -mCellSizeX;
462 mCellSizeY = geotransform[5];
463 if ( mCellSizeY < 0 )
465 mCellSizeY = -mCellSizeY;
468 const char *projection = GDALGetProjectionRef( inputDataset );
469 GDALSetProjection( outputDataset.get(), projection );
471 return outputDataset;
477 int nCellsX, nCellsY;
485 GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset.get(), 1 );
486 if ( !elevationBand )
494 minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
495 minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
497 if ( !minOk || !maxOk )
499 GDALComputeRasterMinMax( elevationBand,
true, minMax );
505 double frequency[252] = {0};
506 const double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
508 float *scanLine = (
float * ) CPLMalloc(
sizeof(
float ) * nCellsX );
509 int elevationClass = -1;
511 for (
int i = 0; i < nCellsY; ++i )
513 if ( GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1,
514 scanLine, nCellsX, 1, GDT_Float32,
520 for (
int j = 0; j < nCellsX; ++j )
522 elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange );
523 if ( elevationClass >= 0 && elevationClass < 252 )
525 frequency[elevationClass] += 1.0;
533 for (
int i = 0; i < 252; ++i )
535 frequency[i] = std::log10( frequency[i] );
539 QFile outFile( file );
540 if ( !outFile.open( QIODevice::WriteOnly | QIODevice::Truncate ) )
545 QTextStream outstream( &outFile );
546#if QT_VERSION < QT_VERSION_CHECK(6, 0, 0)
547 outstream.setCodec(
"UTF-8" );
549 for (
int i = 0; i < 252; ++i )
551 outstream << QString::number( i ) +
',' + QString::number( frequency[i] ) << Qt::endl;
559 QList< QgsRelief::ReliefColor > resultList;
561 int nCellsX, nCellsY;
569 GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset.get(), 1 );
570 if ( !elevationBand )
578 minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
579 minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
581 if ( !minOk || !maxOk )
583 GDALComputeRasterMinMax( elevationBand,
true, minMax );
589 double frequency[252] = {0};
590 const double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
592 float *scanLine = (
float * ) CPLMalloc(
sizeof(
float ) * nCellsX );
593 int elevationClass = -1;
595 for (
int i = 0; i < nCellsY; ++i )
597 if ( GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1,
598 scanLine, nCellsX, 1, GDT_Float32,
603 for (
int j = 0; j < nCellsX; ++j )
605 elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange );
606 elevationClass = std::max( std::min( elevationClass, 251 ), 0 );
607 frequency[elevationClass] += 1.0;
614 for (
int i = 0; i < 252; ++i )
616 frequency[i] = std::log10( frequency[i] );
620 QList<int> classBreaks;
621 classBreaks.append( 0 );
622 classBreaks.append( 28 );
623 classBreaks.append( 56 );
624 classBreaks.append( 84 );
625 classBreaks.append( 112 );
626 classBreaks.append( 140 );
627 classBreaks.append( 168 );
628 classBreaks.append( 196 );
629 classBreaks.append( 224 );
630 classBreaks.append( 252 );
632 for (
int i = 0; i < 10; ++i )
634 optimiseClassBreaks( classBreaks, frequency );
638 for (
int i = 0; i < classBreaks.size(); ++i )
640 qWarning(
"%d", classBreaks[i] );
644 QVector<QColor> colorList;
645 colorList.reserve( 9 );
646 colorList.push_back( QColor( 7, 165, 144 ) );
647 colorList.push_back( QColor( 12, 221, 162 ) );
648 colorList.push_back( QColor( 33, 252, 183 ) );
649 colorList.push_back( QColor( 247, 252, 152 ) );
650 colorList.push_back( QColor( 252, 196, 8 ) );
651 colorList.push_back( QColor( 252, 166, 15 ) );
652 colorList.push_back( QColor( 175, 101, 15 ) );
653 colorList.push_back( QColor( 255, 133, 92 ) );
654 colorList.push_back( QColor( 204, 204, 204 ) );
656 resultList.reserve( classBreaks.size() );
657 for (
int i = 1; i < classBreaks.size(); ++i )
659 const double minElevation = minMax[0] + classBreaks[i - 1] * frequencyClassRange;
660 const double maxElevation = minMax[0] + classBreaks[i] * frequencyClassRange;
667void QgsRelief::optimiseClassBreaks( QList<int> &breaks,
double *frequencies )
669 const int nClasses = breaks.size() - 1;
670 double *a =
new double[nClasses];
671 double *b =
new double[nClasses];
673 for (
int i = 0; i < nClasses; ++i )
676 QList< QPair < int, double > > regressionInput;
677 regressionInput.reserve( breaks.at( i + 1 ) - breaks.at( i ) );
678 for (
int j = breaks.at( i ); j < breaks.at( i + 1 ); ++j )
680 regressionInput.push_back( qMakePair( j, frequencies[j] ) );
683 double aParam, bParam;
684 if ( !regressionInput.isEmpty() && calculateRegression( regressionInput, aParam, bParam ) )
696 const QList<int> classesToRemove;
699 for (
int i = 1; i < nClasses ; ++i )
701 if ( breaks[i] == breaks[ i - 1 ] )
712 int newX = ( b[i - 1] - b[ i ] ) / ( a[ i ] - a[ i - 1 ] );
714 if ( newX <= breaks[i - 1] )
716 newX = breaks[i - 1];
719 else if ( i < nClasses - 1 && newX >= breaks[i + 1] )
721 newX = breaks[i + 1];
729 for (
int i = classesToRemove.size() - 1; i >= 0; --i )
731 breaks.removeAt( classesToRemove.at( i ) );
738int QgsRelief::frequencyClassForElevation(
double elevation,
double minElevation,
double elevationClassRange )
740 return ( elevation - minElevation ) / elevationClassRange;
743bool QgsRelief::calculateRegression(
const QList< QPair < int, double > > &input,
double &a,
double &b )
748 QList< QPair < int, double > >::const_iterator inputIt = input.constBegin();
749 for ( ; inputIt != input.constEnd(); ++inputIt )
751 xSum += inputIt->first;
752 ySum += inputIt->second;
754 xMean = xSum / input.size();
755 yMean = ySum / input.size();
757 double sumCounter = 0;
758 double sumDenominator = 0;
759 inputIt = input.constBegin();
760 for ( ; inputIt != input.constEnd(); ++inputIt )
762 sumCounter += ( ( inputIt->first - xMean ) * ( inputIt->second - yMean ) );
763 sumDenominator += ( ( inputIt->first - xMean ) * ( inputIt->first - xMean ) );
766 a = sumCounter / sumDenominator;
767 b = yMean - a * xMean;
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.
bool isCanceled() const
Tells whether the operation has been canceled already.
void setProgress(double progress)
Sets the current progress for the feedback object.
static bool supportsRasterCreate(GDALDriverH driver)
Reads whether a driver supports GDALCreate() for raster purposes.
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)
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.
void addReliefColorClass(const QgsRelief::ReliefColor &color)
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)
#define QgsDebugError(str)