[gdal-dev] Issues with SetColorTable() and multi-sample TIFF

[afernandez]

Hello Even,
Thank you for the quick reply. I'll try using GDAL subdatasets later this afternoon (although I'm not sure if this solution would become a bit messy for hundreds of files). The format KEA is completely new to me. I was able to install KEALib but found only one source of information online (https://gdal.org/drivers/raster/kea.html). To start somewhere, I simply modified the first 2 lines of the script to:
driver_name = 'KEA'
driver = gdal.GetDriverByName(driver_name)
However, it doesn't pick it up (driver returns 'None'). Ideally there would be some example that I could follow but I might have to do this piecewise in the meantime. 
Even Rouault wrote:
Hi,
The TIFF format only supports only single color table per TIFF image, and it only makes sense for a single-sample TIFF. An alternative would be for you to create multiple TIFF image file directories in the same file (mapped as GDAL subdatasets), with one band per subdataset, by using the APPEND_SUBDATASET=YES creation option.
Another option is to use the KEA format (specialization of HDF5) which supports all the possibilities of the GDAL abstract data model.
Even
Le 14/02/2023 à 17:25, afernandez a écrit :
Hello,
I'm using the GDAL Tools plugin within QGIS to plot some raster layers. The snippet reads:
driver_name = 'GTIFF'
driver = gdal.GetDriverByName(driver_name) 
...
gdal_ds = driver.Create("my_arguments") 
gdal_ds.SetProjection("my_projection")
gdal_ds.SetGeoTransform("my_geo_transform")
for band_idx in range(1, times + 1):
band = gdal_ds.GetRasterBand(band_idx) 
band.SetNoDataValue(no_data)
band.SetDescription(time_step)
data = var[band_idx - 1]
data = ma.getdata(data)
band.WriteArray(data.astype(np_dtype, copy=False))
gdal_ds.FlushCache()
This works without any problem and plots the raster in greys. However, I need the plot to be colored so I modified the code based on some posts, which now reads: 
gdal_ds = driver.Create(out_path, cols, rows, times, gdal.GDT_UInt16) 
gdal_ds.SetProjection("my_projection")
gdal_ds.SetGeoTransform("my_geo_transform")
for band_idx in range(1, times + 1):
band = gdal_ds.GetRasterBand(band_idx) 
band.SetNoDataValue(no_data)
band.SetDescription(time_step)
data = var[band_idx - 1]
data = ma.getdata(data)
data_scaled = np.interp(data, (data.min(), data.max()), (0, 255))
data_scaled2 = data_scaled.astype(int)
colors = gdal.ColorTable()
colors.CreateColorRamp(0, (229, 245, 249), 127, (153, 216, 201))
colors.CreateColorRamp(127, (153, 216, 201), 254, (44, 162, 95))
band.SetRasterColorTable(colors)
band.SetRasterColorInterpretation(gdal.GCI_PaletteIndex)
band.WriteArray(data_scaled2)
gdal_ds.FlushCache()
The problem is that this produces a RunTimeError: SetColorTable() not supported for multi-sample TIFF. I simplified this for a single band and QGIS did plot the colored raster. However, this is not a solution because production runs will generate dozens (or even hundreds) of time step bands. I'm pretty new to GDAL so I'm unsure if this is a feature under current development or if there are other alternatives so any suggestion on how to proceed would be appreciated.
Thank you in advance.
_______________________________________________ gdal-dev mailing list gdal-dev at lists.osgeo.org <mailto:gdal-dev at lists.osgeo.org> https://lists.osgeo.org/mailman/listinfo/gdal-dev <https://lists.osgeo.org/mailman/listinfo/gdal-dev> -- http://www.spatialys.com <http://www.spatialys.com> My software is free, but my time generally not. 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20230214/10928f6e/attachment-0001.htm>