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

PG Raster: fix nodata values with sparse data #56363

Merged
merged 1 commit into from
Feb 15, 2024
Merged

PG Raster: fix nodata values with sparse data #56363

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
52 changes: 48 additions & 4 deletions src/providers/postgres/raster/qgspostgresrasterprovider.cpp
View file
Open in desktop
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
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -629,6 +629,41 @@ def testSparseRaster(self):
critical_postgis_logs = list(filter(lambda log: log[2] == Qgis.MessageLevel.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()
Loading