25 #include "cpl_string.h" 32 #include <QTextStream> 35 : mInputFile( inputFile )
36 , mOutputFile( outputFile )
37 , mOutputFormat( outputFormat )
40 , mInputNodataValue( -1 )
41 , mOutputNodataValue( -1 )
44 mSlopeFilter =
new QgsSlopeFilter( inputFile, outputFile, outputFormat );
45 mAspectFilter =
new QgsAspectFilter( inputFile, outputFile, outputFormat );
46 mHillshadeFilter285 =
new QgsHillshadeFilter( inputFile, outputFile, outputFormat, 285, 30 );
47 mHillshadeFilter300 =
new QgsHillshadeFilter( inputFile, outputFile, outputFormat, 300, 30 );
48 mHillshadeFilter315 =
new QgsHillshadeFilter( inputFile, outputFile, outputFormat, 315, 30 );
58 delete mHillshadeFilter285;
59 delete mHillshadeFilter300;
60 delete mHillshadeFilter315;
65 mReliefColors.clear();
70 mReliefColors.push_back( color );
73 void QgsRelief::setDefaultReliefColors()
95 GDALDriverH outputDriver = openOutputDriver();
102 if ( !outputDataset )
125 GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset.get(), 1 );
130 mInputNodataValue = GDALGetRasterNoDataValue( rasterBand,
nullptr );
137 GDALRasterBandH outputRedBand = GDALGetRasterBand( outputDataset.get(), 1 );
138 GDALRasterBandH outputGreenBand = GDALGetRasterBand( outputDataset.get(), 2 );
139 GDALRasterBandH outputBlueBand = GDALGetRasterBand( outputDataset.get(), 3 );
141 if ( !outputRedBand || !outputGreenBand || !outputBlueBand )
146 GDALSetRasterNoDataValue( outputRedBand, -9999 );
147 GDALSetRasterNoDataValue( outputGreenBand, -9999 );
148 GDALSetRasterNoDataValue( outputBlueBand, -9999 );
149 mOutputNodataValue = GDALGetRasterNoDataValue( outputRedBand,
nullptr );
162 float *scanLine1 = (
float * ) CPLMalloc(
sizeof(
float ) * xSize );
163 float *scanLine2 = (
float * ) CPLMalloc(
sizeof(
float ) * xSize );
164 float *scanLine3 = (
float * ) CPLMalloc(
sizeof(
float ) * xSize );
166 unsigned char *resultRedLine = (
unsigned char * ) CPLMalloc(
sizeof(
unsigned char ) * xSize );
167 unsigned char *resultGreenLine = (
unsigned char * ) CPLMalloc(
sizeof(
unsigned char ) * xSize );
168 unsigned char *resultBlueLine = (
unsigned char * ) CPLMalloc(
sizeof(
unsigned char ) * xSize );
173 for (
int i = 0; i < ySize; ++i )
177 feedback->
setProgress( 100.0 * i / static_cast< double >( ySize ) );
188 for (
int a = 0; a < xSize; ++a )
190 scanLine1[a] = mInputNodataValue;
192 if ( GDALRasterIO( rasterBand, GF_Read, 0, 0, xSize, 1, scanLine2, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
200 CPLFree( scanLine1 );
201 scanLine1 = scanLine2;
202 scanLine2 = scanLine3;
203 scanLine3 = (
float * ) CPLMalloc(
sizeof(
float ) * xSize );
206 if ( i == ySize - 1 )
208 for (
int a = 0; a < xSize; ++a )
210 scanLine3[a] = mInputNodataValue;
215 if ( GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, scanLine3, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
221 for (
int j = 0; j < xSize; ++j )
225 resultOk = processNineCellWindow( &mInputNodataValue, &scanLine1[j], &scanLine1[j + 1], &mInputNodataValue, &scanLine2[j], \
226 &scanLine2[j + 1], &mInputNodataValue, &scanLine3[j], &scanLine3[j + 1], \
227 &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] );
229 else if ( j == xSize - 1 )
231 resultOk = processNineCellWindow( &scanLine1[j - 1], &scanLine1[j], &mInputNodataValue, &scanLine2[j - 1], &scanLine2[j], \
232 &mInputNodataValue, &scanLine3[j - 1], &scanLine3[j], &mInputNodataValue, \
233 &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] );
237 resultOk = processNineCellWindow( &scanLine1[j - 1], &scanLine1[j], &scanLine1[j + 1], &scanLine2[j - 1], &scanLine2[j], \
238 &scanLine2[j + 1], &scanLine3[j - 1], &scanLine3[j], &scanLine3[j + 1], \
239 &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] );
244 resultRedLine[j] = mOutputNodataValue;
245 resultGreenLine[j] = mOutputNodataValue;
246 resultBlueLine[j] = mOutputNodataValue;
250 if ( GDALRasterIO( outputRedBand, GF_Write, 0, i, xSize, 1, resultRedLine, xSize, 1, GDT_Byte, 0, 0 ) != CE_None )
254 if ( GDALRasterIO( outputGreenBand, GF_Write, 0, i, xSize, 1, resultGreenLine, xSize, 1, GDT_Byte, 0, 0 ) != CE_None )
258 if ( GDALRasterIO( outputBlueBand, GF_Write, 0, i, xSize, 1, resultBlueLine, xSize, 1, GDT_Byte, 0, 0 ) != CE_None )
269 CPLFree( resultRedLine );
270 CPLFree( resultBlueLine );
271 CPLFree( resultGreenLine );
272 CPLFree( scanLine1 );
273 CPLFree( scanLine2 );
274 CPLFree( scanLine3 );
286 bool QgsRelief::processNineCellWindow(
float *x1,
float *x2,
float *x3,
float *x4,
float *x5,
float *x6,
float *x7,
float *x8,
float *x9,
287 unsigned char *red,
unsigned char *green,
unsigned char *blue )
294 float hillShadeValue300 = mHillshadeFilter300->
processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
295 if ( hillShadeValue300 != mOutputNodataValue )
297 if ( !setElevationColor( *x5, &r, &g, &b ) )
299 r = hillShadeValue300;
300 g = hillShadeValue300;
301 b = hillShadeValue300;
305 r = r / 2.0 + hillShadeValue300 / 2.0;
306 g = g / 2.0 + hillShadeValue300 / 2.0;
307 b = b / 2.0 + hillShadeValue300 / 2.0;
312 float hillShadeValue315 = mHillshadeFilter315->
processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
314 if ( hillShadeValue315 != mOutputNodataValue && slope != mOutputNodataValue )
319 r2 = 0 / 2.0 + hillShadeValue315 / 2.0;
320 g2 = 0 / 2.0 + hillShadeValue315 / 2.0;
321 b2 = 0 / 2.0 + hillShadeValue315 / 2.0;
323 else if ( slope >= 1 )
325 int slopeValue = 255 - ( slope / 15.0 * 255.0 );
326 r2 = slopeValue / 2.0 + hillShadeValue315 / 2.0;
327 g2 = slopeValue / 2.0 + hillShadeValue315 / 2.0;
328 b2 = slopeValue / 2.0 + hillShadeValue315 / 2.0;
332 r2 = hillShadeValue315;
333 g2 = hillShadeValue315;
334 b2 = hillShadeValue315;
338 r = r * 0.7 + r2 * 0.3;
339 g = g * 0.7 + g2 * 0.3;
340 b = b * 0.7 + b2 * 0.3;
344 float hillShadeValue285 = mHillshadeFilter285->
processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
346 if ( hillShadeValue285 != mOutputNodataValue && aspect != mOutputNodataValue )
348 double angle_diff = std::fabs( 285 - aspect );
349 if ( angle_diff > 180 )
355 if ( angle_diff < 90 )
357 int aspectVal = ( 1 - std::cos( angle_diff * M_PI / 180 ) ) * 255;
358 r3 = 0.5 * 255 + hillShadeValue315 * 0.5;
359 g3 = 0.5 * 255 + hillShadeValue315 * 0.5;
360 b3 = 0.5 * aspectVal + hillShadeValue315 * 0.5;
364 r3 = 0.5 * 255 + hillShadeValue315 * 0.5;
365 g3 = 0.5 * 255 + hillShadeValue315 * 0.5;
366 b3 = 0.5 * 255 + hillShadeValue315 * 0.5;
369 r = r3 * 0.1 + r * 0.9;
370 g = g3 * 0.1 + g * 0.9;
371 b = b3 * 0.1 + b * 0.9;
374 *red = (
unsigned char )r;
375 *green = (
unsigned char )g;
376 *blue = (
unsigned char )b;
380 bool QgsRelief::setElevationColor(
double elevation,
int *red,
int *green,
int *blue )
382 QList< ReliefColor >::const_iterator reliefColorIt = mReliefColors.constBegin();
383 for ( ; reliefColorIt != mReliefColors.constEnd(); ++reliefColorIt )
385 if ( elevation >= reliefColorIt->minElevation && elevation <= reliefColorIt->maxElevation )
387 const QColor &c = reliefColorIt->color;
404 nCellsX = GDALGetRasterXSize( inputDataset.get() );
405 nCellsY = GDALGetRasterYSize( inputDataset.get() );
408 if ( GDALGetRasterCount( inputDataset.get() ) < 1 )
416 GDALDriverH QgsRelief::openOutputDriver()
418 char **driverMetadata =
nullptr;
421 GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
428 driverMetadata = GDALGetMetadata( outputDriver,
nullptr );
429 if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE,
false ) )
444 int xSize = GDALGetRasterXSize( inputDataset );
445 int ySize = GDALGetRasterYSize( inputDataset );
448 char **papszOptions =
nullptr;
451 papszOptions = CSLSetNameValue( papszOptions,
"COMPRESS",
"PACKBITS" );
454 gdal::dataset_unique_ptr outputDataset( GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), xSize, ySize, 3, GDT_Byte, papszOptions ) );
455 CSLDestroy( papszOptions );
456 papszOptions =
nullptr;
458 if ( !outputDataset )
464 double geotransform[6];
465 if ( GDALGetGeoTransform( inputDataset, geotransform ) != CE_None )
469 GDALSetGeoTransform( outputDataset.get(), geotransform );
472 mCellSizeX = geotransform[1];
473 if ( mCellSizeX < 0 )
475 mCellSizeX = -mCellSizeX;
477 mCellSizeY = geotransform[5];
478 if ( mCellSizeY < 0 )
480 mCellSizeY = -mCellSizeY;
483 const char *projection = GDALGetProjectionRef( inputDataset );
484 GDALSetProjection( outputDataset.get(), projection );
486 return outputDataset;
492 int nCellsX, nCellsY;
500 GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset.get(), 1 );
501 if ( !elevationBand )
509 minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
510 minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
512 if ( !minOk || !maxOk )
514 GDALComputeRasterMinMax( elevationBand,
true, minMax );
520 double frequency[252];
521 double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
523 for (
int i = 0; i < 252; ++i )
528 float *scanLine = (
float * ) CPLMalloc(
sizeof(
float ) * nCellsX );
529 int elevationClass = -1;
531 for (
int i = 0; i < nCellsY; ++i )
533 if ( GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1,
534 scanLine, nCellsX, 1, GDT_Float32,
540 for (
int j = 0; j < nCellsX; ++j )
542 elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange );
543 if ( elevationClass >= 0 )
545 frequency[elevationClass] += 1.0;
553 for (
int i = 0; i < 252; ++i )
555 frequency[i] = std::log10( frequency[i] );
559 QFile outFile( file );
560 if ( !outFile.open( QIODevice::WriteOnly | QIODevice::Truncate ) )
565 QTextStream outstream( &outFile );
566 for (
int i = 0; i < 252; ++i )
568 outstream << QString::number( i ) +
',' + QString::number( frequency[i] ) << endl;
576 QList< QgsRelief::ReliefColor > resultList;
578 int nCellsX, nCellsY;
586 GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset.get(), 1 );
587 if ( !elevationBand )
595 minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
596 minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
598 if ( !minOk || !maxOk )
600 GDALComputeRasterMinMax( elevationBand,
true, minMax );
606 double frequency[252];
607 double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
609 for (
int i = 0; i < 252; ++i )
614 float *scanLine = (
float * ) CPLMalloc(
sizeof(
float ) * nCellsX );
615 int elevationClass = -1;
617 for (
int i = 0; i < nCellsY; ++i )
619 if ( GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1,
620 scanLine, nCellsX, 1, GDT_Float32,
625 for (
int j = 0; j < nCellsX; ++j )
627 elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange );
628 if ( elevationClass < 0 )
632 else if ( elevationClass >= 252 )
634 elevationClass = 251;
636 frequency[elevationClass] += 1.0;
643 for (
int i = 0; i < 252; ++i )
645 frequency[i] = std::log10( frequency[i] );
649 QList<int> classBreaks;
650 classBreaks.append( 0 );
651 classBreaks.append( 28 );
652 classBreaks.append( 56 );
653 classBreaks.append( 84 );
654 classBreaks.append( 112 );
655 classBreaks.append( 140 );
656 classBreaks.append( 168 );
657 classBreaks.append( 196 );
658 classBreaks.append( 224 );
659 classBreaks.append( 252 );
661 for (
int i = 0; i < 10; ++i )
663 optimiseClassBreaks( classBreaks, frequency );
667 for (
int i = 0; i < classBreaks.size(); ++i )
669 qWarning(
"%d", classBreaks[i] );
673 QVector<QColor> colorList;
674 colorList.push_back( QColor( 7, 165, 144 ) );
675 colorList.push_back( QColor( 12, 221, 162 ) );
676 colorList.push_back( QColor( 33, 252, 183 ) );
677 colorList.push_back( QColor( 247, 252, 152 ) );
678 colorList.push_back( QColor( 252, 196, 8 ) );
679 colorList.push_back( QColor( 252, 166, 15 ) );
680 colorList.push_back( QColor( 175, 101, 15 ) );
681 colorList.push_back( QColor( 255, 133, 92 ) );
682 colorList.push_back( QColor( 204, 204, 204 ) );
684 resultList.reserve( classBreaks.size() );
685 for (
int i = 1; i < classBreaks.size(); ++i )
687 double minElevation = minMax[0] + classBreaks[i - 1] * frequencyClassRange;
688 double maxElevation = minMax[0] + classBreaks[i] * frequencyClassRange;
695 void QgsRelief::optimiseClassBreaks( QList<int> &breaks,
double *frequencies )
697 int nClasses = breaks.size() - 1;
698 double *a =
new double[nClasses];
699 double *b =
new double[nClasses];
701 for (
int i = 0; i < nClasses; ++i )
704 QList< QPair < int, double > > regressionInput;
705 for (
int j = breaks.at( i ); j < breaks.at( i + 1 ); ++j )
707 regressionInput.push_back( qMakePair( j, frequencies[j] ) );
710 double aParam, bParam;
711 if ( !regressionInput.isEmpty() && calculateRegression( regressionInput, aParam, bParam ) )
723 QList<int> classesToRemove;
726 for (
int i = 1; i < nClasses ; ++i )
728 if ( breaks[i] == breaks[ i - 1 ] )
739 int newX = ( b[i - 1] - b[ i ] ) / ( a[ i ] - a[ i - 1 ] );
741 if ( newX <= breaks[i - 1] )
743 newX = breaks[i - 1];
746 else if ( i < nClasses - 1 && newX >= breaks[i + 1] )
748 newX = breaks[i + 1];
756 for (
int i = classesToRemove.size() - 1; i >= 0; --i )
758 breaks.removeAt( classesToRemove.at( i ) );
765 int QgsRelief::frequencyClassForElevation(
double elevation,
double minElevation,
double elevationClassRange )
767 return ( elevation - minElevation ) / elevationClassRange;
770 bool QgsRelief::calculateRegression(
const QList< QPair < int, double > > &input,
double &a,
double &b )
775 QList< QPair < int, double > >::const_iterator inputIt = input.constBegin();
776 for ( ; inputIt != input.constEnd(); ++inputIt )
778 xSum += inputIt->first;
779 ySum += inputIt->second;
781 xMean = xSum / input.size();
782 yMean = ySum / input.size();
784 double sumCounter = 0;
785 double sumDenominator = 0;
786 inputIt = input.constBegin();
787 for ( ; inputIt != input.constEnd(); ++inputIt )
789 sumCounter += ( ( inputIt->first - xMean ) * ( inputIt->second - yMean ) );
790 sumDenominator += ( ( inputIt->first - xMean ) * ( inputIt->first - xMean ) );
793 a = sumCounter / sumDenominator;
794 b = yMean - a * xMean;
float processNineCellWindow(float *x11, float *x21, float *x31, float *x12, float *x22, float *x32, float *x13, float *x23, float *x33) override
Calculates output value from nine input values.
void setZFactor(double factor)
QList< QgsRelief::ReliefColor > calculateOptimizedReliefClasses()
Calculates class breaks according with the method of Buenzli (2011) using an iterative algorithm for ...
void setCellSizeY(double size)
void setProgress(double progress)
Sets the current progress for the feedback object.
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Compare two doubles (but allow some difference)
Base class for feedback objects to be used for cancelation of something running in a worker thread...
bool exportFrequencyDistributionToCsv(const QString &file)
Write frequency of elevation values to file for manual inspection.
int processRaster(QgsFeedback *feedback=nullptr)
Starts the calculation, reads from mInputFile and stores the result in mOutputFile.
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...
float processNineCellWindow(float *x11, float *x21, float *x31, float *x12, float *x22, float *x32, float *x13, float *x23, float *x33) override
Calculates output value from nine input values.
void setInputNodataValue(double value)
Calculates aspect values in a window of 3x3 cells based on first order derivatives in x- and y- direc...
void setCellSizeX(double size)
std::unique_ptr< void, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
bool isCanceled() const
Tells whether the operation has been canceled already.
float processNineCellWindow(float *x11, float *x21, float *x31, float *x12, float *x22, float *x32, float *x13, float *x23, float *x33) override
Calculates output value from nine input values.
Calculates slope values in a window of 3x3 cells based on first order derivatives in x- and y- direct...
void setOutputNodataValue(double value)
QgsRelief(const QString &inputFile, const QString &outputFile, const QString &outputFormat)
void addReliefColorClass(const QgsRelief::ReliefColor &color)