Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix gdal sample identify/sample precision issue #45086

Merged
merged 2 commits into from
Sep 17, 2021
Merged

Fix gdal sample identify/sample precision issue #45086

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b381abb
GDAL Use proper type when sampling raster
elpaso Sep 16, 2021
c59cffb
Don't crash collecting stats when band is not valid
elpaso Sep 16, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
77 changes: 74 additions & 3 deletions src/core/providers/gdal/qgsgdalprovider.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1440,9 +1440,80 @@ double QgsGdalProvider::sample( const QgsPointXY &point, int band, bool *ok, con
if ( !worldToPixel( point.x(), point.y(), col, row ) )
return std::numeric_limits<double>::quiet_NaN();

float value = 0;
CPLErr err = GDALRasterIO( hBand, GF_Read, col, row, 1, 1,
&value, 1, 1, GDT_Float32, 0, 0 );
double value{0};
CPLErr err {CE_Failure};

const GDALDataType dataType {GDALGetRasterDataType( hBand )};
switch ( dataType )
{
case GDT_Byte:
{
unsigned char tempVal{0};
err = GDALRasterIO( hBand, GF_Read, col, row, 1, 1,
&tempVal, 1, 1, dataType, 0, 0 );
value = static_cast<double>( tempVal );
break;
}
case GDT_UInt16:
{
uint16_t tempVal{0};
err = GDALRasterIO( hBand, GF_Read, col, row, 1, 1,
&tempVal, 1, 1, dataType, 0, 0 );
value = static_cast<double>( tempVal );
break;
}
case GDT_Int16:
{
int16_t tempVal{0};
err = GDALRasterIO( hBand, GF_Read, col, row, 1, 1,
&tempVal, 1, 1, dataType, 0, 0 );
value = static_cast<double>( tempVal );
break;
}
case GDT_UInt32:
{
uint32_t tempVal{0};
err = GDALRasterIO( hBand, GF_Read, col, row, 1, 1,
&tempVal, 1, 1, dataType, 0, 0 );
value = static_cast<double>( tempVal );
break;
}
case GDT_Int32:
{
int32_t tempVal{0};
err = GDALRasterIO( hBand, GF_Read, col, row, 1, 1,
&tempVal, 1, 1, dataType, 0, 0 );
value = static_cast<double>( tempVal );
break;
}
case GDT_Float32:
{
float tempVal{0};
err = GDALRasterIO( hBand, GF_Read, col, row, 1, 1,
&tempVal, 1, 1, dataType, 0, 0 );
value = static_cast<double>( tempVal );
break;
}
case GDT_Float64:
{
// No need to cast for double
err = GDALRasterIO( hBand, GF_Read, col, row, 1, 1,
&value, 1, 1, dataType, 0, 0 );
break;
}
case GDT_CInt16:
case GDT_CInt32:
case GDT_CFloat32:
case GDT_CFloat64:
QgsDebugMsg( QStringLiteral( "Raster complex data type is not supported!" ) );
break;
case GDT_TypeCount:
QgsDebugMsg( QStringLiteral( "Raster data type GDT_TypeCount is not supported!" ) );
break;
case GDT_Unknown:
QgsDebugMsg( QStringLiteral( "Raster data type is unknown!" ) );
break;
}
if ( err != CE_None )
return std::numeric_limits<double>::quiet_NaN();

Expand Down
196 changes: 100 additions & 96 deletions src/core/raster/qgspalettedrasterrenderer.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -509,135 +509,139 @@ QgsPalettedRasterRenderer::ClassData QgsPalettedRasterRenderer::classDataFromRas
return ClassData();

ClassData data;
qlonglong numClasses = 0;

if ( feedback )
feedback->setProgress( 0 );

// Collect unique values for float rasters
if ( raster->dataType( bandNumber ) == Qgis::DataType::Float32 || raster->dataType( bandNumber ) == Qgis::DataType::Float64 )
if ( bandNumber > 0 && bandNumber <= raster->bandCount() )
{
qlonglong numClasses = 0;

if ( feedback && feedback->isCanceled() )
{
return data;
}
if ( feedback )
feedback->setProgress( 0 );

std::set<double> values;
// Collect unique values for float rasters
if ( raster->dataType( bandNumber ) == Qgis::DataType::Float32 || raster->dataType( bandNumber ) == Qgis::DataType::Float64 )
{

const int maxWidth = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_WIDTH;
const int maxHeight = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_HEIGHT;
if ( feedback && feedback->isCanceled() )
{
return data;
}

QgsRasterIterator iter( raster );
iter.startRasterRead( bandNumber, raster->xSize(), raster->ySize(), raster->extent(), feedback );
std::set<double> values;

const int nbBlocksWidth = static_cast< int >( std::ceil( 1.0 * raster->xSize() / maxWidth ) );
const int nbBlocksHeight = static_cast< int >( std::ceil( 1.0 * raster->ySize() / maxHeight ) );
const int nbBlocks = nbBlocksWidth * nbBlocksHeight;
const int maxWidth = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_WIDTH;
const int maxHeight = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_HEIGHT;

int iterLeft = 0;
int iterTop = 0;
int iterCols = 0;
int iterRows = 0;
std::unique_ptr< QgsRasterBlock > rasterBlock;
QgsRectangle blockExtent;
bool isNoData = false;
while ( iter.readNextRasterPart( bandNumber, iterCols, iterRows, rasterBlock, iterLeft, iterTop, &blockExtent ) )
{
if ( feedback )
feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
QgsRasterIterator iter( raster );
iter.startRasterRead( bandNumber, raster->xSize(), raster->ySize(), raster->extent(), feedback );

if ( feedback && feedback->isCanceled() )
break;
const int nbBlocksWidth = static_cast< int >( std::ceil( 1.0 * raster->xSize() / maxWidth ) );
const int nbBlocksHeight = static_cast< int >( std::ceil( 1.0 * raster->ySize() / maxHeight ) );
const int nbBlocks = nbBlocksWidth * nbBlocksHeight;

for ( int row = 0; row < iterRows; row++ )
int iterLeft = 0;
int iterTop = 0;
int iterCols = 0;
int iterRows = 0;
std::unique_ptr< QgsRasterBlock > rasterBlock;
QgsRectangle blockExtent;
bool isNoData = false;
while ( iter.readNextRasterPart( bandNumber, iterCols, iterRows, rasterBlock, iterLeft, iterTop, &blockExtent ) )
{
if ( feedback )
feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );

if ( feedback && feedback->isCanceled() )
break;

for ( int column = 0; column < iterCols; column++ )
for ( int row = 0; row < iterRows; row++ )
{
if ( feedback && feedback->isCanceled() )
break;

const double currentValue = rasterBlock->valueAndNoData( row, column, isNoData );
if ( numClasses >= MAX_FLOAT_CLASSES )
{
QgsMessageLog::logMessage( QStringLiteral( "Number of classes exceeded maximum (%1)." ).arg( MAX_FLOAT_CLASSES ), QStringLiteral( "Raster" ) );
break;
}
if ( !isNoData && values.find( currentValue ) == values.end() )
for ( int column = 0; column < iterCols; column++ )
{
values.insert( currentValue );
data.push_back( Class( currentValue, QColor(), QLocale().toString( currentValue ) ) );
numClasses++;
if ( feedback && feedback->isCanceled() )
break;

const double currentValue = rasterBlock->valueAndNoData( row, column, isNoData );
if ( numClasses >= MAX_FLOAT_CLASSES )
{
QgsMessageLog::logMessage( QStringLiteral( "Number of classes exceeded maximum (%1)." ).arg( MAX_FLOAT_CLASSES ), QStringLiteral( "Raster" ) );
break;
}
if ( !isNoData && values.find( currentValue ) == values.end() )
{
values.insert( currentValue );
data.push_back( Class( currentValue, QColor(), QLocale().toString( currentValue ) ) );
numClasses++;
}
}
}
}
// must be sorted
std::sort( data.begin(), data.end(), []( const Class & a, const Class & b ) -> bool
{
return a.value < b.value;
} );
}
// must be sorted
std::sort( data.begin(), data.end(), []( const Class & a, const Class & b ) -> bool
{
return a.value < b.value;
} );
}
else
{
// get min and max value from raster
const QgsRasterBandStats stats = raster->bandStatistics( bandNumber, QgsRasterBandStats::Min | QgsRasterBandStats::Max, QgsRectangle(), 0, feedback );
if ( feedback && feedback->isCanceled() )
return ClassData();

const double min = stats.minimumValue;
const double max = stats.maximumValue;
// need count of every individual value
const int bins = std::ceil( max - min ) + 1;
if ( bins <= 0 )
return ClassData();

const QgsRasterHistogram histogram = raster->histogram( bandNumber, bins, min, max, QgsRectangle(), 0, false, feedback );
if ( feedback && feedback->isCanceled() )
return ClassData();

const double interval = ( histogram.maximum - histogram.minimum + 1 ) / histogram.binCount;
double currentValue = histogram.minimum;
for ( int idx = 0; idx < histogram.binCount; ++idx )
else
{
const int count = histogram.histogramVector.at( idx );
if ( count > 0 )
// get min and max value from raster
const QgsRasterBandStats stats = raster->bandStatistics( bandNumber, QgsRasterBandStats::Min | QgsRasterBandStats::Max, QgsRectangle(), 0, feedback );
if ( feedback && feedback->isCanceled() )
return ClassData();

const double min = stats.minimumValue;
const double max = stats.maximumValue;
// need count of every individual value
const int bins = std::ceil( max - min ) + 1;
if ( bins <= 0 )
return ClassData();

const QgsRasterHistogram histogram = raster->histogram( bandNumber, bins, min, max, QgsRectangle(), 0, false, feedback );
if ( feedback && feedback->isCanceled() )
return ClassData();

const double interval = ( histogram.maximum - histogram.minimum + 1 ) / histogram.binCount;
double currentValue = histogram.minimum;
for ( int idx = 0; idx < histogram.binCount; ++idx )
{
data << Class( currentValue, QColor(), QLocale().toString( currentValue ) );
numClasses++;
const int count = histogram.histogramVector.at( idx );
if ( count > 0 )
{
data << Class( currentValue, QColor(), QLocale().toString( currentValue ) );
numClasses++;
}
currentValue += interval;
}
currentValue += interval;
}
}

// assign colors from ramp
if ( ramp && numClasses > 0 )
{
int i = 0;

if ( QgsRandomColorRamp *randomRamp = dynamic_cast<QgsRandomColorRamp *>( ramp ) )
// assign colors from ramp
if ( ramp && numClasses > 0 )
{
//ramp is a random colors ramp, so inform it of the total number of required colors
//this allows the ramp to pregenerate a set of visually distinctive colors
randomRamp->setTotalColorCount( data.count() );
}
int i = 0;

if ( numClasses > 1 )
numClasses -= 1; //avoid duplicate first color
if ( QgsRandomColorRamp *randomRamp = dynamic_cast<QgsRandomColorRamp *>( ramp ) )
{
//ramp is a random colors ramp, so inform it of the total number of required colors
//this allows the ramp to pregenerate a set of visually distinctive colors
randomRamp->setTotalColorCount( data.count() );
}

QgsPalettedRasterRenderer::ClassData::iterator cIt = data.begin();
for ( ; cIt != data.end(); ++cIt )
{
if ( feedback )
if ( numClasses > 1 )
numClasses -= 1; //avoid duplicate first color

QgsPalettedRasterRenderer::ClassData::iterator cIt = data.begin();
for ( ; cIt != data.end(); ++cIt )
{
// Show no less than 1%, then the max between class fill and real progress
feedback->setProgress( std::max<int>( 1, 100 * ( i + 1 ) / numClasses ) );
if ( feedback )
{
// Show no less than 1%, then the max between class fill and real progress
feedback->setProgress( std::max<int>( 1, 100 * ( i + 1 ) / numClasses ) );
}
cIt->color = ramp->color( i / static_cast<double>( numClasses ) );
i++;
}
cIt->color = ramp->color( i / static_cast<double>( numClasses ) );
i++;
}
}
return data;
Expand Down
Loading