Skip to content

Commit

Permalink
PG Raster: fix nodata values with sparse data
Browse files Browse the repository at this point in the history 
Fix #55784
  • Loading branch information
@elpaso @nyalldawson
elpaso authored and nyalldawson committed Feb 16, 2024
1 parent 32cdc33 commit 79f4364
Show file tree
Hide file tree
Showing 2 changed files with 83 additions and 4 deletions.
52 changes: 48 additions & 4 deletions src/providers/postgres/raster/qgspostgresrasterprovider.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -304,7 +304,7 @@ bool QgsPostgresRasterProvider::readBlock( int bandNo, const QgsRectangle &viewE
}

bool ok;
QString val { result.PQgetvalue( 0, 0 ) };
QString val { result.PQntuples() > 0 ? result.PQgetvalue( 0, 0 ) : QString() };

if ( val.isNull() && mUseSrcNoDataValue[bandNo - 1] )
{
Expand Down Expand Up @@ -402,6 +402,9 @@ bool QgsPostgresRasterProvider::readBlock( int bandNo, const QgsRectangle &viewE
}
else // Fetch block
{
const Qgis::DataType dataType { mDataTypes[ bandNo - 1 ] };
const GDALDataType gdalDataType = QgsGdalUtils::gdalDataTypeFromQgisDataType( dataType );
const double noDataValue { mSrcNoDataValue[ bandNo - 1 ] };

const double xRes = viewExtent.width() / width;
const double yRes = viewExtent.height() / height;
Expand Down Expand Up @@ -456,10 +459,47 @@ bool QgsPostgresRasterProvider::readBlock( int bandNo, const QgsRectangle &viewE
if ( tileResponse.tiles.isEmpty() )
{
// rasters can be sparse by omitting some of the blocks/tiles
// so we should not log an error here but make sure
// the result buffer is filled with nodata
gdal::dataset_unique_ptr dstDS { QgsGdalUtils::createSingleBandMemoryDataset(
gdalDataType, viewExtent, width, height, mCrs ) };
if ( ! dstDS )
{
const QString lastError = QString::fromUtf8( CPLGetLastErrorMsg() ) ;
QgsMessageLog::logMessage( tr( "Unable to create destination raster for tiles from %1: %2" )
.arg( tableToQuery, lastError ), tr( "PostGIS" ), Qgis::MessageLevel::Critical );
return false;
}

GDALSetRasterNoDataValue( GDALGetRasterBand( dstDS.get(), 1 ), noDataValue );
// fill with nodata
GDALFillRaster( GDALGetRasterBand( dstDS.get(), 1 ), noDataValue, 0 );

// copy to the result buffer
CPLErrorReset();
CPLErr err = GDALRasterIO( GDALGetRasterBand( dstDS.get(), 1 ),
GF_Read,
0,
0,
width,
height,
data,
width,
height,
gdalDataType,
0,
0 );
if ( err != CE_None )
{
const QString lastError = QString::fromUtf8( CPLGetLastErrorMsg() ) ;
QgsMessageLog::logMessage( tr( "Unable to write raster to block from %1: %2" )
.arg( mQuery, lastError ), tr( "PostGIS" ), Qgis::MessageLevel::Critical );
return false;
}

return true;
}


// Finally merge the tiles
// We must have at least one tile at this point (we checked for that before)

Expand All @@ -469,8 +509,6 @@ bool QgsPostgresRasterProvider::readBlock( int bandNo, const QgsRectangle &viewE
const int tmpWidth = static_cast<int>( std::round( tilesExtent.width() / tileResponse.tiles.first().scaleX ) );
const int tmpHeight = static_cast<int>( std::round( tilesExtent.height() / std::fabs( tileResponse.tiles.first().scaleY ) ) );

GDALDataType gdalDataType { static_cast<GDALDataType>( sourceDataType( bandNo ) ) };

//qDebug() << "Creating output raster: " << tilesExtent.toString() << tmpWidth << tmpHeight;

gdal::dataset_unique_ptr tmpDS { QgsGdalUtils::createSingleBandMemoryDataset(
Expand All @@ -484,8 +522,11 @@ bool QgsPostgresRasterProvider::readBlock( int bandNo, const QgsRectangle &viewE
}
}

GDALSetRasterNoDataValue( GDALGetRasterBand( tmpDS.get(), 1 ), noDataValue );

// Write tiles to the temporary raster
CPLErrorReset();

for ( auto &tile : std::as_const( tileResponse.tiles ) )
{
// Offset in px from the base raster
Expand Down Expand Up @@ -515,6 +556,7 @@ bool QgsPostgresRasterProvider::readBlock( int bandNo, const QgsRectangle &viewE
}
}


#if 0
// Debug output raster content
double pdfMin;
Expand All @@ -536,6 +578,8 @@ bool QgsPostgresRasterProvider::readBlock( int bandNo, const QgsRectangle &viewE
return false;
}

GDALSetRasterNoDataValue( GDALGetRasterBand( dstDS.get(), 1 ), noDataValue );

// Resample the raster to the final bounds and resolution
if ( ! QgsGdalUtils::resampleSingleBandRaster( tmpDS.get(), dstDS.get(), GDALResampleAlg::GRA_NearestNeighbour, nullptr ) )
{
Expand Down
35 changes: 35 additions & 0 deletions tests/src/python/test_provider_postgresraster.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -630,6 +630,41 @@ def testSparseRaster(self):
critical_postgis_logs = list(filter(lambda log: log[2] == Qgis.Critical and log[1] == "PostGIS", list(log_spy)))
self.assertEqual(len(critical_postgis_logs), 0, list(log_spy))

def testSparseTiles(self):
"""Test issue GH #55784"""

rl = QgsRasterLayer(
self.dbconn + " key=\'rid\' srid=3035 sslmode=disable table={table} schema={schema}".format(
table='raster_sparse_3035', schema='public'), 'pg_layer', 'postgresraster')

self.assertTrue(rl.isValid())

dp = rl.dataProvider()

r = dp.identify(QgsPointXY(4080317.72, 2430635.68), Qgis.RasterIdentifyFormat.Value).results()
self.assertEqual(r[1], -9999.0)

# tile request returned no tiles, check nodata
ext = QgsRectangle.fromCenterAndSize(QgsPointXY(4080317.72, 2430635.68), 1, 1)
b = dp.block(1, ext, 1, 1)
self.assertTrue(b.isValid())
self.assertEqual(b.value(0, 0), -9999.0)
self.assertTrue(b.isNoData(0, 0))

# tile request returned one tile with value and nodata in two adjacent cells
ext = QgsRectangle(4080186, 2430632, 4080216, 2430649)
b = dp.block(1, ext, 2, 1)
self.assertEqual(b.value(0, 1), -9999.0)
self.assertEqual(int(b.value(0, 0)), 223)

# Multiple nodata tiles
ext = QgsRectangle(4080272, 2430686, 4080302, 2430703)
b = dp.block(1, ext, 2, 1)
self.assertEqual(b.value(0, 1), -9999.0)
self.assertEqual(b.value(0, 0), -9999.0)
self.assertTrue(b.isNoData(0, 0))
self.assertTrue(b.isNoData(0, 1))


if __name__ == '__main__':
unittest.main()

0 comments on commit 79f4364

Please sign in to comment.