Skip to content
Permalink
Browse files
GDAL Use proper type when sampling raster
Fixes #44902
  • Loading branch information
elpaso authored and nyalldawson committed Nov 1, 2021
1 parent 4caf8bb commit 25aaacf9be88a475444284c09ae32ea8d4f4fc3f
Showing with 211 additions and 3 deletions.
  1. +74 −3 src/core/providers/gdal/qgsgdalprovider.cpp
  2. +137 −0 tests/src/python/test_qgsrasterlayer.py
@@ -1457,9 +1457,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();

@@ -110,6 +110,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')

0 comments on commit 25aaacf

Please sign in to comment.