[gdal-dev] Trouble with GPKG and colors

[Ben Avrahami]

Hello, I'm trying to save a byte raster to a geopackage. I ran into a
problem that the geopackage raster driver does not behave as expected when
working with bytes and color tables. I have written a python script that
explains the issue (attached):

The script creates 4 rasters in similar ways but prints different outputs:

GTIFF(no palette):  Minimum=0, Maximum=9, Color Table=False, band count=1
GTIFF(palette):  Minimum=0, Maximum=9, Color Table=True, band count=1
GPKG(no palette):  Minimum=0, Maximum=9, Color Table=False, band count=4
GPKG(palette):  Minimum=0, Maximum=255, Color Table=False, band count=4

The problem is that even though we set a color table in the geopackage when
writing, it reads as having no color table when reading. Worse yet, the
original raster values are lost when written with a color table.

Is there any workaround or fix for this? We want to be able to keep
non-image byte data with a palette to geopackage.

Regards, Ben
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20181029/8ce220cd/attachment.html>
-------------- next part --------------
import gdal
import numpy as np

template_array = np.array([[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]])

wgs84_geo_wkt = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]'

color_table = gdal.ColorTable()
for i in range(10):
    r = int((9 - i) * 255 / 9)
    color_table.SetColorEntry(i, (r, 127, 127, 255))


def write(driver, path, palette):
    ds = gdal.GetDriverByName(driver).Create(path, *template_array.shape[::-1], 1, gdal.GDT_Byte)
    ds.SetGeoTransform((30, 0.01, 0.0, 30, 0.0, -0.01))
    ds.SetProjection(wgs84_geo_wkt)
    band = ds.GetRasterBand(1)
    band.WriteArray(template_array)
    if palette:
        band.SetColorTable(color_table)
    band.FlushCache()
    del band
    del ds


def read(name, palette, path):
    ds = gdal.OpenEx(path, gdal.OF_READONLY)
    band = ds.GetRasterBand(1)
    stats = band.GetStatistics(True, True)
    ct = band.GetColorTable()
    band_count = ds.RasterCount
    print(
        f"{name}({'palette' if palette else 'no palette'}):  Minimum={stats[0]:.0f}, Maximum={stats[1]:.0f},"
        f" Color Table={ct is not None}, band count={band_count}")


for name, path, pal in [
    ('GTIFF', r'D:\temp\demo_tif.tif', False),
    ('GTIFF', r'D:\temp\demo_tif_palette.tif', True),
    ('GPKG', r'D:\temp\demo_gpkg.gpkg', False),
    ('GPKG', r'D:\temp\demo_gpkg_palette.gpkg', True),
]:
    write(name, path, pal)
    read(name, pal, path)