gpkg_2d_gridded_coverage_ancillary table uncompletely filled when encoding elevation values in Int16 with a scale factor

[gcebelieu]

What is the bug?

I'm leading some experimentations with the use integer encoding for DEM in GeoPackage + coverage gridded data.

For this I'm using gdal_translate to convert a GeoTIFF file with Float32 encoded elevation values to a GeoPackage coverage gridded data layer encoded with Integer16 values within PNG tiles with a scale factor of 0.01 in order to keep a centimetric accuracy.

I'm using this command line:

gdal_translate arctic_60_41_2_2_2m_v4-1_dem-5041-3855-lzw3.tif arctic_60_41_2_2_2m_v4-1_dem-5041-3855.gpkg -co RASTER_TABLE=Artic_DEM -co VERSION=1.3 -co GRID_CELL_ENCODING=grid-value-is-center -co CRS_WKT_EXTENSION=YES -ot Int16 -co TILE_FORMAT=PNG -scale 0 1 0 100 -ovr NONE -co BLOCKSIZE=251 -co UOM=m

As a result, the values are correctly scaled and encoded in the PNG Tiles but the gpkg_2d_gridded_coverage_ancillary table doesn't reflect the scale factor that was applied.

As a consequence when visualizing data into QGiS, the values of elevations of the GeoPackage file are 100 times greater than the one of the GeoTIFF File.

I have tried to change this using the -a_scale and -a_offset options of gdal_translate, but this was ineffective.

If I want them to be correctly interpreted I have to manually modify the coverage ancillary table:

UPDATE 
gpkg_2d_gridded_coverage_ancillary SET 
scale = 0.01, 
offset=-327.68, 
precision=0.01 
WHERE tile_matrix_set_name = 'Arctic_DEM';

Then the values are correctly rendered in QGiS.

What I consider as a bug or an incomplete feature of gdal_translate is that it doesn't offer a way to completely fill the gpkg_2d_gridded_coverage_ancillary table regarding the scale factor applied to the values.

Steps to reproduce the issue

The steps to reproduce the behavior are documented in this short report of the experimentation (with link to original DEM file):

GeoPackage-Integer-encoding-experiment.pdf

Versions and provenance

Windows 11, tested with gdal versions 3.10 and 3.12

Additional context

No response

[rouault]

What I consider as a bug or an incomplete feature of gdal_translate is that it doesn't offer a way to completely fill the gpkg_2d_gridded_coverage_ancillary table regarding the scale factor applied to the values.

Just rewording your observations:

• The scaling & offseting you apply at the gdal_translate level is intended to be user visible, so GDAL behaves as designed here. Said otherwise "-scale 0 1 0 100" means multiply input values by 100 and make sure those 100x values are the ones you get when reading the output raster.
• The GPKG driver indeed hardcodes on creation the scale/offset using for Float32->Int16 conversion and doesn't offer control to the user for customizing that. What you express in that ticket would rather be the addition of GPKG creation options to control those values