Re: [gdal-dev] Trouble with GPKG and colors

2018-10-29 Thread Idan Miara 
Hi Even,

Thanks for you prompt answer!

Could you please point out a valid example on how to read/write a single
byte raster dataset (preferably with a color table) in GPKG with any
supported setup (tiling, format, padding etc)?

GTiff also can store different band count, padding, etc, but there isn't a
similar problem with it. I'm not sure I understand the difference.
So I'm not sure I understand the RGBA expansion you were talking about, do
you have a link to more info on that matter?

Kind regards,
Idan

On Mon, 29 Oct 2018 at 14:29, Even Rouault 
wrote:

> On lundi 29 octobre 2018 14:19:12 CET Ben Avrahami wrote:
> > 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.
> >
>
> Change
> ds = gdal.OpenEx(path, gdal.OF_READONLY)
> to
> ds = gdal.OpenEx(path, gdal.OF_READONLY, open_options=['BAND_COUNT=1'])
>
> if you know that all the tiles have a color table. Note: that will not
> work if
> your raster is made of a single tile that has some padding (which seems to
> be
> the case in your demo code), because the driver has to expand to RGBA to
> be
> able to fill transparent pixels.
>
> The fact that you don't get excatly the characteristics of the input
> dataset
> comes from the fact that GPKG can store tiles of different formats, band
> count, etc, so the driver defaults to exposing at RGBA as a safe setting
>
> --
> Spatialys - Geospatial professional services
> http://www.spatialys.com
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Trouble with GPKG and colors

2018-10-29 Thread Even Rouault 
On lundi 29 octobre 2018 14:19:12 CET Ben Avrahami wrote:
> 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.
> 

Change
ds = gdal.OpenEx(path, gdal.OF_READONLY)
to
ds = gdal.OpenEx(path, gdal.OF_READONLY, open_options=['BAND_COUNT=1'])

if you know that all the tiles have a color table. Note: that will not work if 
your raster is made of a single tile that has some padding (which seems to be 
the case in your demo code), because the driver has to expand to RGBA to be 
able to fill transparent pixels.

The fact that you don't get excatly the characteristics of the input dataset 
comes from the fact that GPKG can store tiles of different formats, band 
count, etc, so the driver defaults to exposing at RGBA as a safe setting

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Trouble with GPKG and colors

2018-10-29 Thread 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
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)
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev