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], &scanLine2[j + 1], &mInputNodataValue, &scanLine3[j], &scanLine3[j + 1], &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] );
215 else if ( j == xSize - 1 )
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] );
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] );
226 resultRedLine[j] = mOutputNodataValue;
227 resultGreenLine[j] = mOutputNodataValue;
228 resultBlueLine[j] = mOutputNodataValue;
232 if ( GDALRasterIO( outputRedBand, GF_Write, 0, i, xSize, 1, resultRedLine, xSize, 1, GDT_Byte, 0, 0 ) != CE_None )
236 if ( GDALRasterIO( outputGreenBand, GF_Write, 0, i, xSize, 1, resultGreenLine, xSize, 1, GDT_Byte, 0, 0 ) != CE_None )
240 if ( GDALRasterIO( outputBlueBand, GF_Write, 0, i, xSize, 1, resultBlueLine, xSize, 1, GDT_Byte, 0, 0 ) != CE_None )
251 CPLFree( resultRedLine );
252 CPLFree( resultBlueLine );
253 CPLFree( resultGreenLine );
254 CPLFree( scanLine1 );
255 CPLFree( scanLine2 );
256 CPLFree( scanLine3 );
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 )
275 const float hillShadeValue300 = mHillshadeFilter300->processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
276 if ( hillShadeValue300 != mOutputNodataValue )
278 if ( !getElevationColor( *x5, &r, &g, &b ) )
280 r = hillShadeValue300;
281 g = hillShadeValue300;
282 b = hillShadeValue300;
286 r = r / 2.0 + hillShadeValue300 / 2.0;
287 g = g / 2.0 + hillShadeValue300 / 2.0;
288 b = b / 2.0 + hillShadeValue300 / 2.0;
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 )
300 r2 = 0 / 2.0 + hillShadeValue315 / 2.0;
301 g2 = 0 / 2.0 + hillShadeValue315 / 2.0;
302 b2 = 0 / 2.0 + hillShadeValue315 / 2.0;
304 else if ( slope >= 1 )
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;
313 r2 = hillShadeValue315;
314 g2 = hillShadeValue315;
315 b2 = hillShadeValue315;
319 r = r * 0.7 + r2 * 0.3;
320 g = g * 0.7 + g2 * 0.3;
321 b = b * 0.7 + b2 * 0.3;
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 )
329 double angle_diff = std::fabs( 285 - aspect );
330 if ( angle_diff > 180 )
336 if ( angle_diff < 90 )
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;
345 r3 = 0.5 * 255 + hillShadeValue315 * 0.5;
346 g3 = 0.5 * 255 + hillShadeValue315 * 0.5;
347 b3 = 0.5 * 255 + hillShadeValue315 * 0.5;
350 r = r3 * 0.1 + r * 0.9;
351 g = g3 * 0.1 + g * 0.9;
352 b = b3 * 0.1 + b * 0.9;
355 *red = (
unsigned char ) r;
356 *green = (
unsigned char ) g;
357 *blue = (
unsigned char ) b;
361bool QgsRelief::getElevationColor(
double elevation,
int *red,
int *green,
int *blue )
const
363 QList<ReliefColor>::const_iterator reliefColorIt = mReliefColors.constBegin();
364 for ( ; reliefColorIt != mReliefColors.constEnd(); ++reliefColorIt )
366 if ( elevation >= reliefColorIt->minElevation && elevation <= reliefColorIt->maxElevation )
368 const QColor &
c = reliefColorIt->color;
385 nCellsX = GDALGetRasterXSize( inputDataset.get() );
386 nCellsY = GDALGetRasterYSize( inputDataset.get() );
389 if ( GDALGetRasterCount( inputDataset.get() ) < 1 )
397GDALDriverH QgsRelief::openOutputDriver()
400 GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
422 const int xSize = GDALGetRasterXSize( inputDataset );
423 const int ySize = GDALGetRasterYSize( inputDataset );
426 char **papszOptions =
nullptr;
429 papszOptions = CSLSetNameValue( papszOptions,
"COMPRESS",
"PACKBITS" );
432 gdal::dataset_unique_ptr outputDataset( GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), xSize, ySize, 3, GDT_Byte, papszOptions ) );
433 CSLDestroy( papszOptions );
434 papszOptions =
nullptr;
436 if ( !outputDataset )
442 double geotransform[6];
443 if ( GDALGetGeoTransform( inputDataset, geotransform ) != CE_None )
447 GDALSetGeoTransform( outputDataset.get(), geotransform );
450 mCellSizeX = geotransform[1];
451 if ( mCellSizeX < 0 )
453 mCellSizeX = -mCellSizeX;
455 mCellSizeY = geotransform[5];
456 if ( mCellSizeY < 0 )
458 mCellSizeY = -mCellSizeY;
461 const char *projection = GDALGetProjectionRef( inputDataset );
462 GDALSetProjection( outputDataset.get(), projection );
464 return outputDataset;
470 int nCellsX, nCellsY;
478 GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset.get(), 1 );
479 if ( !elevationBand )
487 minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
488 minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
490 if ( !minOk || !maxOk )
492 GDALComputeRasterMinMax( elevationBand,
true, minMax );
498 double frequency[252] = { 0 };
499 const double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
501 float *scanLine = (
float * ) CPLMalloc(
sizeof(
float ) * nCellsX );
502 int elevationClass = -1;
504 for (
int i = 0; i < nCellsY; ++i )
506 if ( GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 ) != CE_None )
511 for (
int j = 0; j < nCellsX; ++j )
513 elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange );
514 if ( elevationClass >= 0 && elevationClass < 252 )
516 frequency[elevationClass] += 1.0;
524 for (
int i = 0; i < 252; ++i )
526 frequency[i] = std::log10( frequency[i] );
530 QFile outFile( file );
531 if ( !outFile.open( QIODevice::WriteOnly | QIODevice::Truncate ) )
536 QTextStream outstream( &outFile );
537#if QT_VERSION < QT_VERSION_CHECK( 6, 0, 0 )
538 outstream.setCodec(
"UTF-8" );
540 for (
int i = 0; i < 252; ++i )
542 outstream << QString::number( i ) +
',' + QString::number( frequency[i] ) << Qt::endl;
550 QList<QgsRelief::ReliefColor> resultList;
552 int nCellsX, nCellsY;
560 GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset.get(), 1 );
561 if ( !elevationBand )
569 minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
570 minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
572 if ( !minOk || !maxOk )
574 GDALComputeRasterMinMax( elevationBand,
true, minMax );
580 double frequency[252] = { 0 };
581 const double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
583 float *scanLine = (
float * ) CPLMalloc(
sizeof(
float ) * nCellsX );
584 int elevationClass = -1;
586 for (
int i = 0; i < nCellsY; ++i )
588 if ( GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 ) != CE_None )
592 for (
int j = 0; j < nCellsX; ++j )
594 elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange );
595 elevationClass = std::max( std::min( elevationClass, 251 ), 0 );
596 frequency[elevationClass] += 1.0;
603 for (
int i = 0; i < 252; ++i )
605 frequency[i] = std::log10( frequency[i] );
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 );
621 for (
int i = 0; i < 10; ++i )
623 optimiseClassBreaks( classBreaks, frequency );
627 for (
int i = 0; i < classBreaks.size(); ++i )
629 qWarning(
"%d", classBreaks[i] );
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 ) );
645 resultList.reserve( classBreaks.size() );
646 for (
int i = 1; i < classBreaks.size(); ++i )
648 const double minElevation = minMax[0] + classBreaks[i - 1] * frequencyClassRange;
649 const double maxElevation = minMax[0] + classBreaks[i] * frequencyClassRange;
656void QgsRelief::optimiseClassBreaks( QList<int> &breaks,
double *frequencies )
658 const int nClasses = breaks.size() - 1;
659 double *a =
new double[nClasses];
660 double *b =
new double[nClasses];
662 for (
int i = 0; i < nClasses; ++i )
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 )
669 regressionInput.push_back( qMakePair( j, frequencies[j] ) );
672 double aParam, bParam;
673 if ( !regressionInput.isEmpty() && calculateRegression( regressionInput, aParam, bParam ) )
685 const QList<int> classesToRemove;
688 for (
int i = 1; i < nClasses; ++i )
690 if ( breaks[i] == breaks[i - 1] )
701 int newX = ( b[i - 1] - b[i] ) / ( a[i] - a[i - 1] );
703 if ( newX <= breaks[i - 1] )
705 newX = breaks[i - 1];
708 else if ( i < nClasses - 1 && newX >= breaks[i + 1] )
710 newX = breaks[i + 1];
718 for (
int i = classesToRemove.size() - 1; i >= 0; --i )
720 breaks.removeAt( classesToRemove.at( i ) );
727int QgsRelief::frequencyClassForElevation(
double elevation,
double minElevation,
double elevationClassRange )
729 return ( elevation - minElevation ) / elevationClassRange;
732bool QgsRelief::calculateRegression(
const QList<QPair<int, double>> &input,
double &a,
double &b )
737 QList<QPair<int, double>>::const_iterator inputIt = input.constBegin();
738 for ( ; inputIt != input.constEnd(); ++inputIt )
740 xSum += inputIt->first;
741 ySum += inputIt->second;
743 xMean = xSum / input.size();
744 yMean = ySum / input.size();
746 double sumCounter = 0;
747 double sumDenominator = 0;
748 inputIt = input.constBegin();
749 for ( ; inputIt != input.constEnd(); ++inputIt )
751 sumCounter += ( ( inputIt->first - xMean ) * ( inputIt->second - yMean ) );
752 sumDenominator += ( ( inputIt->first - xMean ) * ( inputIt->first - xMean ) );
755 a = sumCounter / sumDenominator;
756 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)