From b381abb79a1095dfca8a398197f6457e031e9a60 Mon Sep 17 00:00:00 2001 From: Alessandro Pasotti Date: Thu, 16 Sep 2021 12:30:49 +0200 Subject: [PATCH 1/2] GDAL Use proper type when sampling raster Fixes #44902 --- src/core/providers/gdal/qgsgdalprovider.cpp | 77 ++++++++++- tests/src/python/test_qgsrasterlayer.py | 137 ++++++++++++++++++++ 2 files changed, 211 insertions(+), 3 deletions(-) diff --git a/src/core/providers/gdal/qgsgdalprovider.cpp b/src/core/providers/gdal/qgsgdalprovider.cpp index 679af95d3779..840298d66a7f 100644 --- a/src/core/providers/gdal/qgsgdalprovider.cpp +++ b/src/core/providers/gdal/qgsgdalprovider.cpp @@ -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::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( 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( 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( 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( 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( 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( 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::quiet_NaN(); diff --git a/tests/src/python/test_qgsrasterlayer.py b/tests/src/python/test_qgsrasterlayer.py index 30fcbe4b9f57..703626e9ba93 100644 --- a/tests/src/python/test_qgsrasterlayer.py +++ b/tests/src/python/test_qgsrasterlayer.py @@ -114,6 +114,143 @@ def testIdentify(self): myMessage = 'Expected: %s\nGot: %s' % (myValues, myExpectedValues) self.assertEqual(myValues, myExpectedValues, myMessage) + def testSampleIdentify(self): + """Test that sample() and identify() return the same values, GH #44902""" + + tempdir = QTemporaryDir() + temppath = os.path.join(tempdir.path(), 'test_sample.tif') + + def _test(qgis_data_type): + rlayer = QgsRasterLayer(temppath, 'test_sample') + self.assertTrue(rlayer.isValid()) + self.assertEqual(rlayer.dataProvider().dataType(1), qgis_data_type) + + x = 252290.5 + for y in [5022000.5, 5022001.5]: + pos = QgsPointXY(x, y) + value_sample = rlayer.dataProvider().sample(pos, 1)[0] + value_identify = rlayer.dataProvider().identify(pos, QgsRaster.IdentifyFormatValue).results()[1] + # Check values for UInt32 + if qgis_data_type == Qgis.UInt32: + if y == 5022000.5: + self.assertEqual(value_sample, 4294967000.0) + else: + self.assertEqual(value_sample, 4294967293.0) + self.assertEqual(value_sample, value_identify) + # print(value_sample, value_identify) + + # Test GDT_UInt32 + driver = gdal.GetDriverByName('GTiff') + outRaster = driver.Create(temppath, 3, 3, 1, gdal.GDT_UInt32) + outRaster.SetGeoTransform((252290.0, 1.0, 0.0, 5022002.0, 0.0, -1.0)) + outband = outRaster.GetRasterBand(1) + npdata = np.array([[4294967293, 1000, 5], + [4294967000, 50, 4]], dtype=np.uint32) + outband.WriteArray(npdata) + outband.FlushCache() + outRaster.FlushCache() + del outRaster + + _test(Qgis.UInt32) + + # Test GDT_Int32 + driver = gdal.GetDriverByName('GTiff') + outRaster = driver.Create(temppath, 3, 3, 1, gdal.GDT_Int32) + outRaster.SetGeoTransform((252290.0, 1.0, 0.0, 5022002.0, 0.0, -1.0)) + outband = outRaster.GetRasterBand(1) + npdata = np.array([[1294967293, 1000, 5], + [1294967000, 50, 4]], dtype=np.int32) + outband.WriteArray(npdata) + outband.FlushCache() + outRaster.FlushCache() + del outRaster + + _test(Qgis.Int32) + + # Test GDT_Float32 + driver = gdal.GetDriverByName('GTiff') + outRaster = driver.Create(temppath, 3, 3, 1, gdal.GDT_Float32) + outRaster.SetGeoTransform((252290.0, 1.0, 0.0, 5022002.0, 0.0, -1.0)) + outband = outRaster.GetRasterBand(1) + npdata = np.array([[1294967293, 1000, 5], + [1294967000, 50, 4]], dtype=np.float32) + outband.WriteArray(npdata) + outband.FlushCache() + outRaster.FlushCache() + del outRaster + + _test(Qgis.Float32) + + # Test GDT_Float64 + driver = gdal.GetDriverByName('GTiff') + outRaster = driver.Create(temppath, 3, 3, 1, gdal.GDT_Float64) + outRaster.SetGeoTransform((252290.0, 1.0, 0.0, 5022002.0, 0.0, -1.0)) + outband = outRaster.GetRasterBand(1) + npdata = np.array([[1294967293, 1000, 5], + [1294967000, 50, 4]], dtype=np.float64) + outband.WriteArray(npdata) + outband.FlushCache() + outRaster.FlushCache() + del outRaster + + _test(Qgis.Float64) + + # Test GDT_Uint16 + driver = gdal.GetDriverByName('GTiff') + outRaster = driver.Create(temppath, 3, 3, 1, gdal.GDT_UInt16) + outRaster.SetGeoTransform((252290.0, 1.0, 0.0, 5022002.0, 0.0, -1.0)) + outband = outRaster.GetRasterBand(1) + npdata = np.array([[1294967293, 1000, 5], + [1294967000, 50, 4]], dtype=np.uint16) + outband.WriteArray(npdata) + outband.FlushCache() + outRaster.FlushCache() + del outRaster + + _test(Qgis.UInt16) + + # Test GDT_Int16 + driver = gdal.GetDriverByName('GTiff') + outRaster = driver.Create(temppath, 3, 3, 1, gdal.GDT_Int16) + outRaster.SetGeoTransform((252290.0, 1.0, 0.0, 5022002.0, 0.0, -1.0)) + outband = outRaster.GetRasterBand(1) + npdata = np.array([[31768, 1000, 5], + [12345, 50, 4]], dtype=np.int16) + outband.WriteArray(npdata) + outband.FlushCache() + outRaster.FlushCache() + del outRaster + + _test(Qgis.Int16) + + # Test GDT_Int32 + driver = gdal.GetDriverByName('GTiff') + outRaster = driver.Create(temppath, 3, 3, 1, gdal.GDT_Int32) + outRaster.SetGeoTransform((252290.0, 1.0, 0.0, 5022002.0, 0.0, -1.0)) + outband = outRaster.GetRasterBand(1) + npdata = np.array([[31768, 1000, 5], + [12345, 50, 4]], dtype=np.int32) + outband.WriteArray(npdata) + outband.FlushCache() + outRaster.FlushCache() + del outRaster + + _test(Qgis.Int32) + + # Test GDT_Byte + driver = gdal.GetDriverByName('GTiff') + outRaster = driver.Create(temppath, 3, 3, 1, gdal.GDT_Byte) + outRaster.SetGeoTransform((252290.0, 1.0, 0.0, 5022002.0, 0.0, -1.0)) + outband = outRaster.GetRasterBand(1) + npdata = np.array([[123, 255, 5], + [1, 50, 4]], dtype=np.byte) + outband.WriteArray(npdata) + outband.FlushCache() + outRaster.FlushCache() + del outRaster + + _test(Qgis.Byte) + def testTransparency(self): myPath = os.path.join(unitTestDataPath('raster'), 'band1_float32_noct_epsg4326.tif') From c59cffbfd160a6e2fe8cd03c40d193d2b11bf906 Mon Sep 17 00:00:00 2001 From: Alessandro Pasotti Date: Thu, 16 Sep 2021 12:31:36 +0200 Subject: [PATCH 2/2] Don't crash collecting stats when band is not valid unreported --- src/core/raster/qgspalettedrasterrenderer.cpp | 196 +++++++++--------- 1 file changed, 100 insertions(+), 96 deletions(-) diff --git a/src/core/raster/qgspalettedrasterrenderer.cpp b/src/core/raster/qgspalettedrasterrenderer.cpp index df47b8674422..b202eed70e24 100644 --- a/src/core/raster/qgspalettedrasterrenderer.cpp +++ b/src/core/raster/qgspalettedrasterrenderer.cpp @@ -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 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 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( 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( 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( 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( 1, 100 * ( i + 1 ) / numClasses ) ); + } + cIt->color = ramp->color( i / static_cast( numClasses ) ); + i++; } - cIt->color = ramp->color( i / static_cast( numClasses ) ); - i++; } } return data;