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 );
97 GDALDriverH outputDriver = openOutputDriver();
103 GDALDatasetH outputDataset = openOutputFile( inputDataset, outputDriver );
104 if ( !outputDataset )
127 GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, 1 );
130 GDALClose( inputDataset );
131 GDALClose( outputDataset );
134 mInputNodataValue = GDALGetRasterNoDataValue( rasterBand,
nullptr );
141 GDALRasterBandH outputRedBand = GDALGetRasterBand( outputDataset, 1 );
142 GDALRasterBandH outputGreenBand = GDALGetRasterBand( outputDataset, 2 );
143 GDALRasterBandH outputBlueBand = GDALGetRasterBand( outputDataset, 3 );
145 if ( !outputRedBand || !outputGreenBand || !outputBlueBand )
147 GDALClose( inputDataset );
148 GDALClose( outputDataset );
152 GDALSetRasterNoDataValue( outputRedBand, -9999 );
153 GDALSetRasterNoDataValue( outputGreenBand, -9999 );
154 GDALSetRasterNoDataValue( outputBlueBand, -9999 );
155 mOutputNodataValue = GDALGetRasterNoDataValue( outputRedBand,
nullptr );
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;
334 g2 = hillShadeValue315;
335 b2 = hillShadeValue315;
339 r = r * 0.7 + r2 * 0.3;
340 g = g * 0.7 + g2 * 0.3;
341 b = b * 0.7 + b2 * 0.3;
345 float hillShadeValue285 = mHillshadeFilter285->
processNineCellWindow( x1, x2, x3, x4, x5, x6, x7, x8, x9 );
347 if ( hillShadeValue285 != mOutputNodataValue && aspect != mOutputNodataValue )
349 double angle_diff = qAbs( 285 - aspect );
350 if ( angle_diff > 180 )
356 if ( angle_diff < 90 )
358 int aspectVal = ( 1 - cos( angle_diff *
M_PI / 180 ) ) * 255;
359 r3 = 0.5 * 255 + hillShadeValue315 * 0.5;
360 g3 = 0.5 * 255 + hillShadeValue315 * 0.5;
361 b3 = 0.5 * aspectVal + hillShadeValue315 * 0.5;
365 r3 = 0.5 * 255 + hillShadeValue315 * 0.5;
366 g3 = 0.5 * 255 + hillShadeValue315 * 0.5;
367 b3 = 0.5 * 255 + hillShadeValue315 * 0.5;
370 r = r3 * 0.1 + r * 0.9;
371 g = g3 * 0.1 + g * 0.9;
372 b = b3 * 0.1 + b * 0.9;
375 *red = (
unsigned char )r;
376 *green = (
unsigned char )g;
377 *blue = (
unsigned char )b;
381 bool QgsRelief::setElevationColor(
double elevation,
int* red,
int* green,
int* blue )
384 for ( ; reliefColorIt != mReliefColors.constEnd(); ++reliefColorIt )
386 if ( elevation >= reliefColorIt->minElevation && elevation <= reliefColorIt->maxElevation )
388 const QColor& c = reliefColorIt->color;
400 GDALDatasetH QgsRelief::openInputFile(
int& nCellsX,
int& nCellsY )
405 nCellsX = GDALGetRasterXSize( inputDataset );
406 nCellsY = GDALGetRasterYSize( inputDataset );
409 if ( GDALGetRasterCount( inputDataset ) < 1 )
411 GDALClose( inputDataset );
418 GDALDriverH QgsRelief::openOutputDriver()
420 char **driverMetadata;
423 GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.
toLocal8Bit().
data() );
430 driverMetadata = GDALGetMetadata( outputDriver,
nullptr );
431 if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE,
false ) )
446 int xSize = GDALGetRasterXSize( inputDataset );
447 int ySize = GDALGetRasterYSize( inputDataset );
450 char **papszOptions =
nullptr;
453 papszOptions = CSLSetNameValue( papszOptions,
"COMPRESS",
"PACKBITS" );
456 GDALDatasetH outputDataset = GDALCreate( outputDriver,
TO8F( mOutputFile ), xSize, ySize, 3, GDT_Byte, papszOptions );
457 if ( !outputDataset )
459 return outputDataset;
463 double geotransform[6];
464 if ( GDALGetGeoTransform( inputDataset, geotransform ) != CE_None )
466 GDALClose( outputDataset );
469 GDALSetGeoTransform( outputDataset, 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, projection );
486 return outputDataset;
492 int nCellsX, nCellsY;
493 GDALDatasetH inputDataset = openInputFile( nCellsX, nCellsY );
500 GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset, 1 );
501 if ( !elevationBand )
503 GDALClose( inputDataset );
510 minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
511 minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
513 if ( !minOk || !maxOk )
515 GDALComputeRasterMinMax( elevationBand,
true, minMax );
521 double frequency[252];
522 double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
524 for (
int i = 0; i < 252; ++i )
529 float* scanLine = (
float * ) CPLMalloc(
sizeof(
float ) * nCellsX );
530 int elevationClass = -1;
532 for (
int i = 0; i < nCellsY; ++i )
534 GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1,
535 scanLine, nCellsX, 1, GDT_Float32,
537 for (
int j = 0; j < nCellsX; ++j )
539 elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange );
540 if ( elevationClass >= 0 )
542 frequency[elevationClass] += 1.0;
550 for (
int i = 0; i < 252; ++i )
552 frequency[i] = log10( frequency[i] );
556 QFile outFile( file );
557 if ( !outFile.
open( QIODevice::WriteOnly ) )
563 for (
int i = 0; i < 252; ++i )
575 int nCellsX, nCellsY;
576 GDALDatasetH inputDataset = openInputFile( nCellsX, nCellsY );
583 GDALRasterBandH elevationBand = GDALGetRasterBand( inputDataset, 1 );
584 if ( !elevationBand )
586 GDALClose( inputDataset );
593 minMax[0] = GDALGetRasterMinimum( elevationBand, &minOk );
594 minMax[1] = GDALGetRasterMaximum( elevationBand, &maxOk );
596 if ( !minOk || !maxOk )
598 GDALComputeRasterMinMax( elevationBand,
true, minMax );
604 double frequency[252];
605 double frequencyClassRange = ( minMax[1] - minMax[0] ) / 252.0;
607 for (
int i = 0; i < 252; ++i )
612 float* scanLine = (
float * ) CPLMalloc(
sizeof(
float ) * nCellsX );
613 int elevationClass = -1;
615 for (
int i = 0; i < nCellsY; ++i )
617 GDALRasterIO( elevationBand, GF_Read, 0, i, nCellsX, 1,
618 scanLine, nCellsX, 1, GDT_Float32,
620 for (
int j = 0; j < nCellsX; ++j )
622 elevationClass = frequencyClassForElevation( scanLine[j], minMax[0], frequencyClassRange );
623 if ( elevationClass < 0 )
627 else if ( elevationClass >= 252 )
629 elevationClass = 251;
631 frequency[elevationClass] += 1.0;
638 for (
int i = 0; i < 252; ++i )
640 frequency[i] = log10( frequency[i] );
649 classBreaks.
append( 112 );
650 classBreaks.
append( 140 );
651 classBreaks.
append( 168 );
652 classBreaks.
append( 196 );
653 classBreaks.
append( 224 );
654 classBreaks.
append( 252 );
656 for (
int i = 0; i < 10; ++i )
658 optimiseClassBreaks( classBreaks, frequency );
662 for (
int i = 0; i < classBreaks.
size(); ++i )
664 qWarning(
"%d", classBreaks[i] );
680 for (
int i = 1; i < classBreaks.
size(); ++i )
682 double minElevation = minMax[0] + classBreaks[i - 1] * frequencyClassRange;
683 double maxElevation = minMax[0] + classBreaks[i] * frequencyClassRange;
690 void QgsRelief::optimiseClassBreaks(
QList<int>& breaks,
double* frequencies )
692 int nClasses = breaks.
size() - 1;
693 double* a =
new double[nClasses];
694 double* b =
new double[nClasses];
696 for (
int i = 0; i < nClasses; ++i )
700 for (
int j = breaks.
at( i ); j < breaks.
at( i + 1 ); ++j )
702 regressionInput.
push_back( qMakePair( j, frequencies[j] ) );
705 double aParam, bParam;
706 if ( !regressionInput.
isEmpty() && calculateRegression( regressionInput, aParam, bParam ) )
721 for (
int i = 1; i < nClasses ; ++i )
723 if ( breaks[i] == breaks[ i - 1 ] )
734 int newX = ( b[i - 1] - b[ i ] ) / ( a[ i ] - a[ i - 1 ] );
736 if ( newX <= breaks[i - 1] )
738 newX = breaks[i - 1];
741 else if ( i < nClasses - 1 && newX >= breaks[i + 1] )
743 newX = breaks[i + 1];
751 for (
int i = classesToRemove.
size() - 1; i >= 0; --i )
760 int QgsRelief::frequencyClassForElevation(
double elevation,
double minElevation,
double elevationClassRange )
762 return ( elevation - minElevation ) / elevationClassRange;
771 for ( ; inputIt != input.
constEnd(); ++inputIt )
773 xSum += inputIt->
first;
774 ySum += inputIt->second;
776 xMean = xSum / input.
size();
777 yMean = ySum / input.size();
779 double sumCounter = 0;
780 double sumDenominator = 0;
782 for ( ; inputIt != input.
constEnd(); ++inputIt )
784 sumCounter += (( inputIt->
first - xMean ) * ( inputIt->second - yMean ) );
785 sumDenominator += (( inputIt->
first - xMean ) * ( inputIt->
first - xMean ) );
788 a = sumCounter / sumDenominator;
789 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
const T & at(int i) 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)
void push_back(const T &value)
QgsRelief(const QString &inputFile, const QString &outputFile, const QString &outputFormat)
const_iterator constEnd() const
const_iterator constBegin() const