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()
37 : mInputFile( inputFile )
38 , mOutputFile( outputFile )
39 , mOutputFormat( outputFormat )
42 , mInputNodataValue( -1 )
43 , mOutputNodataValue( -1 )
46 mSlopeFilter =
new QgsSlopeFilter( inputFile, outputFile, outputFormat );
47 mAspectFilter =
new QgsAspectFilter( inputFile, outputFile, outputFormat );
48 mHillshadeFilter285 =
new QgsHillshadeFilter( inputFile, outputFile, outputFormat, 285, 30 );
49 mHillshadeFilter300 =
new QgsHillshadeFilter( inputFile, outputFile, outputFormat, 300, 30 );
50 mHillshadeFilter315 =
new QgsHillshadeFilter( inputFile, outputFile, outputFormat, 315, 30 );
60 delete mHillshadeFilter285;
61 delete mHillshadeFilter300;
62 delete mHillshadeFilter315;
67 mReliefColors.clear();
72 mReliefColors.push_back( color );
75 void QgsRelief::setDefaultReliefColors()
90 GDALDatasetH inputDataset = openInputFile( xSize, ySize );
91 if ( inputDataset == NULL )
97 GDALDriverH outputDriver = openOutputDriver();
98 if ( outputDriver == 0 )
103 GDALDatasetH outputDataset = openOutputFile( inputDataset, outputDriver );
104 if ( outputDataset == NULL )
127 GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, 1 );
128 if ( rasterBand == NULL )
130 GDALClose( inputDataset );
131 GDALClose( outputDataset );
134 mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, NULL );
141 GDALRasterBandH outputRedBand = GDALGetRasterBand( outputDataset, 1 );
142 GDALRasterBandH outputGreenBand = GDALGetRasterBand( outputDataset, 2 );
143 GDALRasterBandH outputBlueBand = GDALGetRasterBand( outputDataset, 3 );
145 if ( outputRedBand == NULL || outputGreenBand == NULL || outputBlueBand == NULL )
147 GDALClose( inputDataset );
148 GDALClose( outputDataset );
152 GDALSetRasterNoDataValue( outputRedBand, -9999 );
153 GDALSetRasterNoDataValue( outputGreenBand, -9999 );
154 GDALSetRasterNoDataValue( outputBlueBand, -9999 );
155 mOutputNodataValue = GDALGetRasterNoDataValue( outputRedBand, NULL );
164 GDALClose( inputDataset );
165 GDALClose( outputDataset );
170 float* scanLine1 = (
float * ) CPLMalloc(
sizeof(
float ) * xSize );
171 float* scanLine2 = (
float * ) CPLMalloc(
sizeof(
float ) * xSize );
172 float* scanLine3 = (
float * ) CPLMalloc(
sizeof(
float ) * xSize );
174 unsigned char* resultRedLine = (
unsigned char * ) CPLMalloc(
sizeof(
unsigned char ) * xSize );
175 unsigned char* resultGreenLine = (
unsigned char * ) CPLMalloc(
sizeof(
unsigned char ) * xSize );
176 unsigned char* resultBlueLine = (
unsigned char * ) CPLMalloc(
sizeof(
unsigned char ) * xSize );
186 for (
int i = 0; i < ySize; ++i )
201 for (
int a = 0; a < xSize; ++a )
203 scanLine1[a] = mInputNodataValue;
205 GDALRasterIO( rasterBand, GF_Read, 0, 0, xSize, 1, scanLine2, xSize, 1, GDT_Float32, 0, 0 );
210 CPLFree( scanLine1 );
211 scanLine1 = scanLine2;
212 scanLine2 = scanLine3;
213 scanLine3 = (
float * ) CPLMalloc(
sizeof(
float ) * xSize );
216 if ( i == ySize - 1 )
218 for (
int a = 0; a < xSize; ++a )
220 scanLine3[a] = mInputNodataValue;
225 GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, scanLine3, xSize, 1, GDT_Float32, 0, 0 );
228 for (
int j = 0; j < xSize; ++j )
232 resultOk = processNineCellWindow( &mInputNodataValue, &scanLine1[j], &scanLine1[j+1], &mInputNodataValue, &scanLine2[j], \
233 &scanLine2[j+1], &mInputNodataValue, &scanLine3[j], &scanLine3[j+1], \
234 &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] );
236 else if ( j == xSize - 1 )
238 resultOk = processNineCellWindow( &scanLine1[j-1], &scanLine1[j], &mInputNodataValue, &scanLine2[j-1], &scanLine2[j], \
239 &mInputNodataValue, &scanLine3[j-1], &scanLine3[j], &mInputNodataValue, \
240 &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] );
244 resultOk = processNineCellWindow( &scanLine1[j-1], &scanLine1[j], &scanLine1[j+1], &scanLine2[j-1], &scanLine2[j], \
245 &scanLine2[j+1], &scanLine3[j-1], &scanLine3[j], &scanLine3[j+1], \
246 &resultRedLine[j], &resultGreenLine[j], &resultBlueLine[j] );
251 resultRedLine[j] = mOutputNodataValue;
252 resultGreenLine[j] = mOutputNodataValue;
253 resultBlueLine[j] = mOutputNodataValue;
257 GDALRasterIO( outputRedBand, GF_Write, 0, i, xSize, 1, resultRedLine, xSize, 1, GDT_Byte, 0, 0 );
258 GDALRasterIO( outputGreenBand, GF_Write, 0, i, xSize, 1, resultGreenLine, xSize, 1, GDT_Byte, 0, 0 );
259 GDALRasterIO( outputBlueBand, GF_Write, 0, i, xSize, 1, resultBlueLine, xSize, 1, GDT_Byte, 0, 0 );
267 CPLFree( resultRedLine );
268 CPLFree( resultBlueLine );
269 CPLFree( resultGreenLine );
270 CPLFree( scanLine1 );
271 CPLFree( scanLine2 );
272 CPLFree( scanLine3 );
274 GDALClose( inputDataset );
279 GDALDeleteDataset( outputDriver,
TO8F( mOutputFile ) );
282 GDALClose( outputDataset );
287 bool QgsRelief::processNineCellWindow(
float* x1,
float* x2,
float* x3,
float* x4,
float* x5,
float* x6,
float* x7,
float* x8,
float* x9,
288 unsigned char* red,
unsigned char* green,
unsigned char* blue )
295 float hillShadeValue300 = mHillshadeFilter300->
processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
296 if ( hillShadeValue300 != mOutputNodataValue )
298 if ( !setElevationColor( *x5, &r, &g, &b ) )
300 r = hillShadeValue300;
301 g = hillShadeValue300;
302 b = hillShadeValue300;
306 r = r / 2.0 + hillShadeValue300 / 2.0;
307 g = g / 2.0 + hillShadeValue300 / 2.0;
308 b = b / 2.0 + hillShadeValue300 / 2.0;
313 float hillShadeValue315 = mHillshadeFilter315->
processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
315 if ( hillShadeValue315 != mOutputNodataValue && slope != mOutputNodataValue )
320 r2 = 0 / 2.0 + hillShadeValue315 / 2.0;
321 g2 = 0 / 2.0 + hillShadeValue315 / 2.0;
322 b2 = 0 / 2.0 + hillShadeValue315 / 2.0;
324 else if ( slope >= 1 )
326 int slopeValue = 255 - ( slope / 15.0 * 255.0 );
327 r2 = slopeValue / 2.0 + hillShadeValue315 / 2.0;
328 g2 = slopeValue / 2.0 + hillShadeValue315 / 2.0;
329 b2 = slopeValue / 2.0 + hillShadeValue315 / 2.0;
333 r2 = hillShadeValue315; g2 = hillShadeValue315; b2 = hillShadeValue315;
337 r = r * 0.7 + r2 * 0.3;
338 g = g * 0.7 + g2 * 0.3;
339 b = b * 0.7 + b2 * 0.3;
343 float hillShadeValue285 = mHillshadeFilter285->
processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
345 if ( hillShadeValue285 != mOutputNodataValue && aspect != mOutputNodataValue )
347 double angle_diff = qAbs( 285 - aspect );
348 if ( angle_diff > 180 )
354 if ( angle_diff < 90 )
356 int aspectVal = ( 1 - cos( angle_diff *
M_PI / 180 ) ) * 255;
357 r3 = 0.5 * 255 + hillShadeValue315 * 0.5;
358 g3 = 0.5 * 255 + hillShadeValue315 * 0.5;
359 b3 = 0.5 * aspectVal + hillShadeValue315 * 0.5;
363 r3 = 0.5 * 255 + hillShadeValue315 * 0.5;
364 g3 = 0.5 * 255 + hillShadeValue315 * 0.5;
365 b3 = 0.5 * 255 + hillShadeValue315 * 0.5;
368 r = r3 * 0.1 + r * 0.9;
369 g = g3 * 0.1 + g * 0.9;
370 b = b3 * 0.1 + b * 0.9;
373 *red = (
unsigned char )r;
374 *green = (
unsigned char )g;
375 *blue = (
unsigned char )b;
379 bool QgsRelief::setElevationColor(
double elevation,
int* red,
int* green,
int* blue )
382 for ( ; reliefColorIt != mReliefColors.constEnd(); ++reliefColorIt )
384 if ( elevation >= reliefColorIt->minElevation && elevation <= reliefColorIt->maxElevation )
386 const QColor& c = reliefColorIt->color;
398 GDALDatasetH QgsRelief::openInputFile(
int& nCellsX,
int& nCellsY )
401 if ( inputDataset != NULL )
403 nCellsX = GDALGetRasterXSize( inputDataset );
404 nCellsY = GDALGetRasterYSize( inputDataset );
407 if ( GDALGetRasterCount( inputDataset ) < 1 )
409 GDALClose( inputDataset );
416 GDALDriverH QgsRelief::openOutputDriver()
418 char **driverMetadata;
421 GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.
toLocal8Bit().
data() );
423 if ( outputDriver == NULL )
428 driverMetadata = GDALGetMetadata( outputDriver, NULL );
429 if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE,
false ) )
439 if ( inputDataset == NULL )
444 int xSize = GDALGetRasterXSize( inputDataset );
445 int ySize = GDALGetRasterYSize( inputDataset );
448 char **papszOptions = NULL;
451 papszOptions = CSLSetNameValue( papszOptions,
"COMPRESS",
"PACKBITS" );
454 GDALDatasetH outputDataset = GDALCreate( outputDriver,
TO8F( mOutputFile ), xSize, ySize, 3, GDT_Byte, papszOptions );
455 if ( outputDataset == NULL )
457 return outputDataset;
461 double geotransform[6];
462 if ( GDALGetGeoTransform( inputDataset, geotransform ) != CE_None )
464 GDALClose( outputDataset );
467 GDALSetGeoTransform( outputDataset, geotransform );
470 mCellSizeX = geotransform[1];
471 if ( mCellSizeX < 0 )
473 mCellSizeX = -mCellSizeX;
475 mCellSizeY = geotransform[5];
476 if ( mCellSizeY < 0 )
478 mCellSizeY = -mCellSizeY;
481 const char* projection = GDALGetProjectionRef( inputDataset );
482 GDALSetProjection( outputDataset, projection );
484 return outputDataset;
490 int nCellsX, nCellsY;
491 GDALDatasetH inputDataset = openInputFile( nCellsX, nCellsY );
492 if ( inputDataset == NULL )
498 GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset, 1 );
499 if ( elevationBand == NULL )
501 GDALClose( inputDataset );
508 minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
509 minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
511 if ( !minOk || !maxOk )
513 GDALComputeRasterMinMax( elevationBand, TRUE, minMax );
519 double frequency[252];
520 double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
522 for (
int i = 0; i < 252; ++i )
527 float* scanLine = (
float * ) CPLMalloc(
sizeof(
float ) * nCellsX );
528 int elevationClass = -1;
530 for (
int i = 0; i < nCellsY; ++i )
532 GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1,
533 scanLine, nCellsX, 1, GDT_Float32,
535 for (
int j = 0; j < nCellsX; ++j )
537 elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange );
538 if ( elevationClass >= 0 )
540 frequency[elevationClass] += 1.0;
548 for (
int i = 0; i < 252; ++i )
550 frequency[i] = log10( frequency[i] );
554 QFile outFile( file );
555 if ( !outFile.
open( QIODevice::WriteOnly ) )
561 for (
int i = 0; i < 252; ++i )
573 int nCellsX, nCellsY;
574 GDALDatasetH inputDataset = openInputFile( nCellsX, nCellsY );
575 if ( inputDataset == NULL )
581 GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset, 1 );
582 if ( elevationBand == NULL )
584 GDALClose( inputDataset );
591 minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
592 minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
594 if ( !minOk || !maxOk )
596 GDALComputeRasterMinMax( elevationBand, TRUE, minMax );
602 double frequency[252];
603 double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
605 for (
int i = 0; i < 252; ++i )
610 float* scanLine = (
float * ) CPLMalloc(
sizeof(
float ) * nCellsX );
611 int elevationClass = -1;
613 for (
int i = 0; i < nCellsY; ++i )
615 GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1,
616 scanLine, nCellsX, 1, GDT_Float32,
618 for (
int j = 0; j < nCellsX; ++j )
620 elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange );
621 if ( elevationClass < 0 )
625 else if ( elevationClass >= 252 )
627 elevationClass = 251;
629 frequency[elevationClass] += 1.0;
636 for (
int i = 0; i < 252; ++i )
638 frequency[i] = log10( frequency[i] );
647 classBreaks.
append( 112 );
648 classBreaks.
append( 140 );
649 classBreaks.
append( 168 );
650 classBreaks.
append( 196 );
651 classBreaks.
append( 224 );
652 classBreaks.
append( 252 );
654 for (
int i = 0; i < 10; ++i )
656 optimiseClassBreaks( classBreaks, frequency );
660 for (
int i = 0; i < classBreaks.
size(); ++i )
662 qWarning(
"%d", classBreaks[i] );
678 for (
int i = 1; i < classBreaks.
size(); ++i )
680 double minElevation = minMax[0] + classBreaks[i - 1] * frequencyClassRange;
681 double maxElevation = minMax[0] + classBreaks[i] * frequencyClassRange;
688 void QgsRelief::optimiseClassBreaks(
QList<int>& breaks,
double* frequencies )
690 int nClasses = breaks.
size() - 1;
691 double* a =
new double[nClasses];
692 double* b =
new double[nClasses];
694 for (
int i = 0; i < nClasses; ++i )
698 for (
int j = breaks.
at( i ); j < breaks.
at( i + 1 ); ++j )
700 regressionInput.
push_back( qMakePair( j, frequencies[j] ) );
703 double aParam, bParam;
704 if ( regressionInput.
size() > 0 && calculateRegression( regressionInput, aParam, bParam ) )
719 for (
int i = 1; i < nClasses ; ++i )
721 if ( breaks[i] == breaks[ i - 1 ] )
732 int newX = ( b[i - 1] - b[ i ] ) / ( a[ i ] - a[ i - 1 ] );
734 if ( newX <= breaks[i - 1] )
736 newX = breaks[i - 1];
739 else if ( i < nClasses - 1 && newX >= breaks[i + 1] )
741 newX = breaks[i + 1];
749 for (
int i = classesToRemove.
size() - 1; i >= 0; --i )
758 int QgsRelief::frequencyClassForElevation(
double elevation,
double minElevation,
double elevationClassRange )
760 return ( elevation - minElevation ) / elevationClassRange;
769 for ( ; inputIt != input.
constEnd(); ++inputIt )
771 xSum += inputIt->
first;
772 ySum += inputIt->second;
774 xMean = xSum / input.
size();
775 yMean = ySum / input.size();
777 double sumCounter = 0;
778 double sumDenominator = 0;
780 for ( ; inputIt != input.
constEnd(); ++inputIt )
782 sumCounter += (( inputIt->
first - xMean ) * ( inputIt->second - yMean ) );
783 sumDenominator += (( inputIt->
first - xMean ) * ( inputIt->
first - xMean ) );
786 a = sumCounter / sumDenominator;
787 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)
void addReliefColorClass(const ReliefColor &color)
int processRaster(QProgressDialog *p)
Starts the calculation, reads from mInputFile and stores the result in mOutputFile.
void setMaximum(int maximum)
void push_back(const T &value)
QList< ReliefColor > calculateOptimizedReliefClasses()
Calculates class breaks according with the method of Buenzli (2011) using an iterative algorithm for ...
void setCellSizeY(double size)
const T & at(int i) const
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
void setValue(int progress)
bool exportFrequencyDistributionToCsv(const QString &file)
Write frequency of elevation values to file for manual inspection.
QString number(int n, int base)
void append(const T &value)
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)
virtual bool open(QFlags< QIODevice::OpenModeFlag > mode)
QByteArray toLocal8Bit() const
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)
const_iterator constEnd() const
const_iterator constBegin() const