2

I am importing a TIFF file into my programm and I am expecting a maximum value for the pixel values of 5.686229 because that's what I have when I read the file with QGIS. But then I get a maximum of 6.171190 when I import it with GDAL (but only with the GDALComputeRasterMinMax function). The values are indeed stretching to 6.171190 when I check manually, and the TIFF I create from this one (without changing the values) have a maximum of 6.171190.

Do you know how the change of values/maximum can be explained?

Here is my code:

GDALAllRegister();
GDALDataset* dataset = ( GDALDataset* ) GDALOpen( dataset_filename.c_str(), GA_ReadOnly);
if( dataset == NULL ){
    printf( "Error: The %s dataset has been incorrectly opened. "
            "Path might be incorrect, or file corrupt.\n",
            dataset_filename.c_str() );
}

printf( "Data name: %s;\n", dataset_filename.c_str() );
printf( "Driver: %s/%s;\n", dataset->GetDriver()->GetDescription(),
        dataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) );

n = dataset->GetRasterXSize();
m = dataset->GetRasterYSize();
printf( "Raster size: %i x %i x %i;\n", n, m, dataset->GetRasterCount() );
int cell_nb = n * m;
printf( "Total number of cells: %i;\n", cell_nb );

const char* projection_ref; // Projection reference of the DEM input.
projection_ref = dataset->GetProjectionRef();
printf( "Projection: `%s´;\n", projection_ref );

importUsingGdalGeostransform( dataset, dataset_filename );

GDALRasterBand* poBand = dataset->GetRasterBand( 1 );
int nBlockXSize, nBlockYSize;
poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
printf( "Block size: %i x %i", nBlockXSize, nBlockYSize );
printf( ", Data type: %s", GDALGetDataTypeName( poBand->GetRasterDataType() ) );
printf( ", Color interpretation: %s;\n", GDALGetColorInterpretationName( poBand->GetColorInterpretation() ) );

With:

void Grid::importUsingGdalGeostransform( GDALDataset* dataset, const string& dataset_filename ){

    double* geo_transform = new double[6];
    if( dataset->GetGeoTransform( geo_transform ) != CE_None ){
         printf( "Error: The %s dataset has no GetGeoTransform in it.\n", dataset_filename.c_str() );
    }

    printf( "Origin: (%f,%f);\n", geo_transform[0], geo_transform[3] );

    float pixel_x_size = fabs( geo_transform[1] );
    float pixel_y_size = fabs( geo_transform[5] );
    printf( "Pixel Size: (%f,%f);\n", pixel_x_size, pixel_y_size );

    if( pixel_x_size * pixel_y_size < 23 ){
        printf( "Warning: The pixel area is of %fm. Make sure the input is correct.\n",
            pixel_x_size * pixel_y_size );
    }

    delete [] geo_transform;
}

And the function triggering different answers:

void Grid::importUsingGdalMinMaxStats( GDALRasterBand* poBand ){
    double adfMinMax[2];
    int bGotMin, bGotMax;
    adfMinMax[0] = poBand->GetMinimum( &bGotMin );
    adfMinMax[1] = poBand->GetMaximum( &bGotMax );
    printf( "Before computation, Min = %f, Max = %f.\n", adfMinMax[0], adfMinMax[1] );
    GDALComputeRasterMinMax( ( GDALRasterBandH )poBand, FALSE, adfMinMax );
    printf( "After computation, Min = %f, Max = %f.\n", adfMinMax[0], adfMinMax[1] );
    adfMinMax[0] = poBand->GetMinimum( &bGotMin );
    adfMinMax[1] = poBand->GetMaximum( &bGotMax );
    printf( "After computation, Min = %f, Max = %f.\n", adfMinMax[0], adfMinMax[1] );
    adfMinMax[1] = GDALGetRasterMaximum( poBand, &bGotMax );
    printf( "After computation, Min = %f, Max = %f, bGotMax %i.\n", adfMinMax[0], adfMinMax[1], bGotMax );
}

And the ouput being:

Data name: /home/name/eclipse-workspace/Library/Artificial/DEM1Step8.tif;
Driver: GTiff/GeoTIFF;
Raster size: 1024 x 512 x 1;
Total number of cells: 524288;
Projection: `PROJCS["WGS 84 / UTM zone 32N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32632"]]´;
Origin: (0.000000,512.000000);
Pixel Size: (1.000000,1.000000);
Warning: The pixel area is of 1.000000m. Make sure the input is correct.
Block size: 1024 x 2, Data type: Float32, Color interpretation: Gray;
Before computation, Min = 0.000000, Max = 5.686229.
After computation, Min = 0.000000, Max = 6.171190.
After computation, Min = 0.000000, Max = 5.686229.
After computation, Min = 0.000000, Max = 5.686229, bGotMax 1.
  • 1
    I've seen this before in Esri where the statistics are calculated with a skip factor, for speed often near enough is good enough so only a subset of pixels are actually calculated, specified by skip factor, for example skip of 3 only uses every 3rd row and column but runs in less than half the time. This may be the case here, or it could be calculating from the pyramids but I can't say for sure. Why not implement your own by reading the raster band and use SetStatistics(double dfMin, double dfMax, double dfMean, double dfStdDev) to enforce the values. – Michael Stimson Nov 28 '18 at 23:21
  • If I read right you did not set the computed min/max back to image as adviced in the documentation of the same kind of function ComputeStatistics. Once computed, the statistics will generally be "set" back on the raster band using SetStatistics().. – user30184 Nov 30 '18 at 7:07
  • There does not seem to be anything like SetRasterMinMax so perhaps you should try with ComputeRasterStatistics and SetRasterStatistics for overwriting the min/max values which are stored into image metadata. – user30184 Nov 30 '18 at 7:18
0

@Michael Stimson, you were right. I am not a very proficient user of QGIS yet, and although I saw the difference between "Cumulative count cut" and "Min / Max", I missed the option "Accuracy" that you are talking about when investigating the values of the image in input.

enter image description here

Strangely, even if I did not set the Accuracy to "Actual" for the image I had in output, it still gave me the real range of values.

I am now sure I have the same value range in input and output, thank you all.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.