23 #include "cpl_string.h"
24 #include <QProgressDialog>
28 #include <QTextStream>
30 #if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800
31 #define TO8F(x) (x).toUtf8().constData()
33 #define TO8F(x) QFile::encodeName( x ).constData()
36 QgsRelief::QgsRelief(
const QString& inputFile,
const QString& outputFile,
const QString& outputFormat ): \
37 mInputFile( inputFile ), mOutputFile( outputFile ), mOutputFormat( outputFormat ), mZFactor( 1.0 )
84 if ( inputDataset == NULL )
91 if ( outputDriver == 0 )
96 GDALDatasetH outputDataset =
openOutputFile( inputDataset, outputDriver );
97 if ( outputDataset == NULL )
120 GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, 1 );
121 if ( rasterBand == NULL )
123 GDALClose( inputDataset );
124 GDALClose( outputDataset );
134 GDALRasterBandH outputRedBand = GDALGetRasterBand( outputDataset, 1 );
135 GDALRasterBandH outputGreenBand = GDALGetRasterBand( outputDataset, 2 );
136 GDALRasterBandH outputBlueBand = GDALGetRasterBand( outputDataset, 3 );
138 if ( outputRedBand == NULL || outputGreenBand == NULL || outputBlueBand == NULL )
140 GDALClose( inputDataset );
141 GDALClose( outputDataset );
145 GDALSetRasterNoDataValue( outputRedBand, -9999 );
146 GDALSetRasterNoDataValue( outputGreenBand, -9999 );
147 GDALSetRasterNoDataValue( outputBlueBand, -9999 );
157 GDALClose( inputDataset );
158 GDALClose( outputDataset );
163 float* scanLine1 = (
float * ) CPLMalloc(
sizeof(
float ) * xSize );
164 float* scanLine2 = (
float * ) CPLMalloc(
sizeof(
float ) * xSize );
165 float* scanLine3 = (
float * ) CPLMalloc(
sizeof(
float ) * xSize );
167 unsigned char* resultRedLine = (
unsigned char * ) CPLMalloc(
sizeof(
unsigned char ) * xSize );
168 unsigned char* resultGreenLine = (
unsigned char * ) CPLMalloc(
sizeof(
unsigned char ) * xSize );
169 unsigned char* resultBlueLine = (
unsigned char * ) CPLMalloc(
sizeof(
unsigned char ) * xSize );
173 p->setMaximum( ySize );
179 for (
int i = 0; i < ySize; ++i )
186 if ( p && p->wasCanceled() )
194 for (
int a = 0; a < xSize; ++a )
198 GDALRasterIO( rasterBand, GF_Read, 0, 0, xSize, 1, scanLine2, xSize, 1, GDT_Float32, 0, 0 );
203 CPLFree( scanLine1 );
204 scanLine1 = scanLine2;
205 scanLine2 = scanLine3;
206 scanLine3 = (
float * ) CPLMalloc(
sizeof(
float ) * xSize );
209 if ( i == ySize - 1 )
211 for (
int a = 0; a < xSize; ++a )
218 GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, scanLine3, xSize, 1, GDT_Float32, 0, 0 );
221 for (
int j = 0; j < xSize; ++j )
227 &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] );
229 else if ( j == xSize - 1 )
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] );
250 GDALRasterIO( outputRedBand, GF_Write, 0, i, xSize, 1, resultRedLine, xSize, 1, GDT_Byte, 0, 0 );
251 GDALRasterIO( outputGreenBand, GF_Write, 0, i, xSize, 1, resultGreenLine, xSize, 1, GDT_Byte, 0, 0 );
252 GDALRasterIO( outputBlueBand, GF_Write, 0, i, xSize, 1, resultBlueLine, xSize, 1, GDT_Byte, 0, 0 );
257 p->setValue( ySize );
260 CPLFree( resultRedLine );
261 CPLFree( resultBlueLine );
262 CPLFree( resultGreenLine );
263 CPLFree( scanLine1 );
264 CPLFree( scanLine2 );
265 CPLFree( scanLine3 );
267 GDALClose( inputDataset );
269 if ( p && p->wasCanceled() )
275 GDALClose( outputDataset );
281 unsigned char* red,
unsigned char* green,
unsigned char* blue )
293 r = hillShadeValue300;
294 g = hillShadeValue300;
295 b = hillShadeValue300;
299 r = r / 2.0 + hillShadeValue300 / 2.0;
300 g = g / 2.0 + hillShadeValue300 / 2.0;
301 b = b / 2.0 + hillShadeValue300 / 2.0;
313 r2 = 0 / 2.0 + hillShadeValue315 / 2.0;
314 g2 = 0 / 2.0 + hillShadeValue315 / 2.0;
315 b2 = 0 / 2.0 + hillShadeValue315 / 2.0;
317 else if ( slope >= 1 )
319 int slopeValue = 255 - ( slope / 15.0 * 255.0 );
320 r2 = slopeValue / 2.0 + hillShadeValue315 / 2.0;
321 g2 = slopeValue / 2.0 + hillShadeValue315 / 2.0;
322 b2 = slopeValue / 2.0 + hillShadeValue315 / 2.0;
326 r2 = hillShadeValue315; g2 = hillShadeValue315; b2 = hillShadeValue315;
330 r = r * 0.7 + r2 * 0.3;
331 g = g * 0.7 + g2 * 0.3;
332 b = b * 0.7 + b2 * 0.3;
340 double angle_diff = qAbs( 285 - aspect );
341 if ( angle_diff > 180 )
347 if ( angle_diff < 90 )
349 int aspectVal = ( 1 - cos( angle_diff *
M_PI / 180 ) ) * 255;
350 r3 = 0.5 * 255 + hillShadeValue315 * 0.5;
351 g3 = 0.5 * 255 + hillShadeValue315 * 0.5;
352 b3 = 0.5 * aspectVal + hillShadeValue315 * 0.5;
356 r3 = 0.5 * 255 + hillShadeValue315 * 0.5;
357 g3 = 0.5 * 255 + hillShadeValue315 * 0.5;
358 b3 = 0.5 * 255 + hillShadeValue315 * 0.5;
361 r = r3 * 0.1 + r * 0.9;
362 g = g3 * 0.1 + g * 0.9;
363 b = b3 * 0.1 + b * 0.9;
366 *red = (
unsigned char )r;
367 *green = (
unsigned char )g;
368 *blue = (
unsigned char )b;
374 QList< ReliefColor >::const_iterator reliefColorIt =
mReliefColors.constBegin();
375 for ( ; reliefColorIt !=
mReliefColors.constEnd(); ++reliefColorIt )
377 if ( elevation >= reliefColorIt->minElevation && elevation <= reliefColorIt->maxElevation )
379 const QColor& c = reliefColorIt->color;
393 GDALDatasetH inputDataset = GDALOpen(
TO8F(
mInputFile ), GA_ReadOnly );
394 if ( inputDataset != NULL )
396 nCellsX = GDALGetRasterXSize( inputDataset );
397 nCellsY = GDALGetRasterYSize( inputDataset );
400 if ( GDALGetRasterCount( inputDataset ) < 1 )
402 GDALClose( inputDataset );
411 char **driverMetadata;
414 GDALDriverH outputDriver = GDALGetDriverByName(
mOutputFormat.toLocal8Bit().data() );
416 if ( outputDriver == NULL )
421 driverMetadata = GDALGetMetadata( outputDriver, NULL );
422 if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE,
false ) )
432 if ( inputDataset == NULL )
437 int xSize = GDALGetRasterXSize( inputDataset );
438 int ySize = GDALGetRasterYSize( inputDataset );;
441 char **papszOptions = NULL;
444 papszOptions = CSLSetNameValue( papszOptions,
"COMPRESS",
"PACKBITS" );
447 GDALDatasetH outputDataset = GDALCreate( outputDriver,
TO8F(
mOutputFile ), xSize, ySize, 3, GDT_Byte, papszOptions );
448 if ( outputDataset == NULL )
450 return outputDataset;
454 double geotransform[6];
455 if ( GDALGetGeoTransform( inputDataset, geotransform ) != CE_None )
457 GDALClose( outputDataset );
460 GDALSetGeoTransform( outputDataset, geotransform );
474 const char* projection = GDALGetProjectionRef( inputDataset );
475 GDALSetProjection( outputDataset, projection );
477 return outputDataset;
483 int nCellsX, nCellsY;
484 GDALDatasetH inputDataset =
openInputFile( nCellsX, nCellsY );
485 if ( inputDataset == NULL )
491 GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset, 1 );
492 if ( elevationBand == NULL )
494 GDALClose( inputDataset );
501 minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
502 minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
504 if ( !minOk || !maxOk )
506 GDALComputeRasterMinMax( elevationBand, TRUE, minMax );
512 double frequency[252];
513 double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
515 for (
int i = 0; i < 252; ++i )
520 float* scanLine = (
float * ) CPLMalloc(
sizeof(
float ) * nCellsX );
521 int elevationClass = -1;
523 for (
int i = 0; i < nCellsY; ++i )
525 GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1,
526 scanLine, nCellsX, 1, GDT_Float32,
528 for (
int j = 0; j < nCellsX; ++j )
531 if ( elevationClass >= 0 )
533 frequency[elevationClass] += 1.0;
541 for (
int i = 0; i < 252; ++i )
543 frequency[i] = log10( frequency[i] );
547 QFile outFile( file );
548 if ( !outFile.open( QIODevice::WriteOnly ) )
553 QTextStream outstream( &outFile );
554 for (
int i = 0; i < 252; ++i )
556 outstream << QString::number( i ) +
"," + QString::number( frequency[i] ) << endl;
564 QList< QgsRelief::ReliefColor > resultList;
566 int nCellsX, nCellsY;
567 GDALDatasetH inputDataset =
openInputFile( nCellsX, nCellsY );
568 if ( inputDataset == NULL )
574 GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset, 1 );
575 if ( elevationBand == NULL )
577 GDALClose( inputDataset );
584 minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
585 minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
587 if ( !minOk || !maxOk )
589 GDALComputeRasterMinMax( elevationBand, TRUE, minMax );
595 double frequency[252];
596 double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
598 for (
int i = 0; i < 252; ++i )
603 float* scanLine = (
float * ) CPLMalloc(
sizeof(
float ) * nCellsX );
604 int elevationClass = -1;
606 for (
int i = 0; i < nCellsY; ++i )
608 GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1,
609 scanLine, nCellsX, 1, GDT_Float32,
611 for (
int j = 0; j < nCellsX; ++j )
614 if ( elevationClass < 0 )
618 else if ( elevationClass >= 252 )
620 elevationClass = 251;
622 frequency[elevationClass] += 1.0;
629 for (
int i = 0; i < 252; ++i )
631 frequency[i] = log10( frequency[i] );
635 QList<int> classBreaks;
636 classBreaks.append( 0 );
637 classBreaks.append( 28 );
638 classBreaks.append( 56 );
639 classBreaks.append( 84 );
640 classBreaks.append( 112 );
641 classBreaks.append( 140 );
642 classBreaks.append( 168 );
643 classBreaks.append( 196 );
644 classBreaks.append( 224 );
645 classBreaks.append( 252 );
647 for (
int i = 0; i < 10; ++i )
653 for (
int i = 0; i < classBreaks.size(); ++i )
655 qWarning(
"%d", classBreaks[i] );
659 QList<QColor> colorList;
660 colorList.push_back( QColor( 7, 165, 144 ) );
661 colorList.push_back( QColor( 12, 221, 162 ) );
662 colorList.push_back( QColor( 33, 252, 183 ) );
663 colorList.push_back( QColor( 247, 252, 152 ) );
664 colorList.push_back( QColor( 252, 196, 8 ) );
665 colorList.push_back( QColor( 252, 166, 15 ) );
666 colorList.push_back( QColor( 175, 101, 15 ) );
667 colorList.push_back( QColor( 255, 133, 92 ) );
668 colorList.push_back( QColor( 204, 204, 204 ) );
670 for (
int i = 1; i < classBreaks.size(); ++i )
672 double minElevation = minMax[0] + classBreaks[i - 1] * frequencyClassRange;
673 double maxElevation = minMax[0] + classBreaks[i] * frequencyClassRange;
682 int nClasses = breaks.size() - 1;
683 double* a =
new double[nClasses];
684 double* b =
new double[nClasses];
686 for (
int i = 0; i < nClasses; ++i )
689 QList< QPair < int, double > > regressionInput;
690 for (
int j = breaks.at( i ); j < breaks.at( i + 1 ); ++j )
692 regressionInput.push_back( qMakePair( j, frequencies[j] ) );
695 double aParam, bParam;
696 if ( regressionInput.size() > 0 &&
calculateRegression( regressionInput, aParam, bParam ) )
708 QList<int> classesToRemove;
711 for (
int i = 1; i < nClasses ; ++i )
713 if ( breaks[i] == breaks[ i - 1 ] )
724 int newX = ( b[i - 1] - b[ i ] ) / ( a[ i ] - a[ i - 1 ] );
726 if ( newX <= breaks[i - 1] )
728 newX = breaks[i - 1];
731 else if ( i < nClasses - 1 && newX >= breaks[i + 1] )
733 newX = breaks[i + 1];
741 for (
int i = classesToRemove.size() - 1; i >= 0; --i )
743 breaks.removeAt( classesToRemove.at( i ) );
752 return ( elevation - minElevation ) / elevationClassRange;
760 QList< QPair < int, double > >::const_iterator inputIt = input.constBegin();
761 for ( ; inputIt != input.constEnd(); ++inputIt )
763 xSum += inputIt->first;
764 ySum += inputIt->second;
766 xMean = xSum / input.size();
767 yMean = ySum / input.size();
769 double sumCounter = 0;
770 double sumDenominator = 0;
771 inputIt = input.constBegin();
772 for ( ; inputIt != input.constEnd(); ++inputIt )
774 sumCounter += (( inputIt->first - xMean ) * ( inputIt->second - yMean ) );
775 sumDenominator += (( inputIt->first - xMean ) * ( inputIt->first - xMean ) );
778 a = sumCounter / sumDenominator;
779 b = yMean - a * xMean;
GDALDriverH openOutputDriver()
Opens the output driver and tests if it supports the creation of a new dataset.
QgsHillshadeFilter * mHillshadeFilter300
void setZFactor(double factor)
void addReliefColorClass(const ReliefColor &color)
int processRaster(QProgressDialog *p)
Starts the calculation, reads from mInputFile and stores the result in mOutputFile.
GDALDatasetH openOutputFile(GDALDatasetH inputDataset, GDALDriverH outputDriver)
Opens the output file and sets the same geotransform and CRS as the input data.
float processNineCellWindow(float *x11, float *x21, float *x31, float *x12, float *x22, float *x32, float *x13, float *x23, float *x33)
Calculates output value from nine input values.
QList< ReliefColor > calculateOptimizedReliefClasses()
Calculates class breaks according with the method of Buenzli (2011) using an iterative algorithm for ...
void setCellSizeY(double size)
QgsSlopeFilter * mSlopeFilter
void setDefaultReliefColors()
Sets relief colors.
float mOutputNodataValue
The nodata value of the output layer.
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
bool exportFrequencyDistributionToCsv(const QString &file)
Write frequency of elevation values to file for manual inspection.
bool calculateRegression(const QList< QPair< int, double > > &input, double &a, double &b)
Calculates coefficients a and b.
float processNineCellWindow(float *x11, float *x21, float *x31, float *x12, float *x22, float *x32, float *x13, float *x23, float *x33)
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...
float mInputNodataValue
The nodata value of the input layer.
void setCellSizeX(double size)
QgsAspectFilter * mAspectFilter
QList< ReliefColor > mReliefColors
bool setElevationColor(double elevation, int *red, int *green, int *blue)
Set elevation color.
QgsHillshadeFilter * mHillshadeFilter285
GDALDatasetH openInputFile(int &nCellsX, int &nCellsY)
Opens the input file and returns the dataset handle and the number of pixels in x-/y- direction...
bool 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)
Calculates slope values in a window of 3x3 cells based on first order derivatives in x- and y- direct...
QgsHillshadeFilter * mHillshadeFilter315
void optimiseClassBreaks(QList< int > &breaks, double *frequencies)
Do one iteration of class break optimisation (algorithm from Garcia and Rodriguez) ...
void setOutputNodataValue(double value)
QgsRelief(const QString &inputFile, const QString &outputFile, const QString &outputFormat)
float processNineCellWindow(float *x11, float *x21, float *x31, float *x12, float *x22, float *x32, float *x13, float *x23, float *x33)
Calculates output value from nine input values.
int frequencyClassForElevation(double elevation, double minElevation, double elevationClassRange)
Returns class (0-255) for an elevation value.