[gdal-dev] gdal_calc.py looses the color lookup table

[Stephen Woodbridge]

Hi All,

I have a geotiff file like:

[] ~/work/oceandata/test$ gdalinfo A2015173174000.L0_LAC.L2_OC.tif -noct
Driver: GTiff/GeoTIFF
Files: A2015173174000.L0_LAC.L2_OC.tif
        A2015173174000.L0_LAC.L2_OC.tif.aux.xml
Size is 2433, 1727
Coordinate System is:
GEOGCS["WGS 84",
     DATUM["WGS_1984",
         SPHEROID["WGS 84",6378137,298.257223563,
             AUTHORITY["EPSG","7030"]],
         AUTHORITY["EPSG","6326"]],
     PRIMEM["Greenwich",0],
     UNIT["degree",0.0174532925199433],
     AUTHORITY["EPSG","4326"]]
Origin = (-84.000000000000000,45.000000000000000)
Pixel Size = (0.012741471640766,-0.012738853693008)
Metadata:
   AREA_OR_POINT=Area
Image Structure Metadata:
   INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -84.0000000,  45.0000000) ( 84d 0' 0.00"W, 45d 0' 0.00"N)
Lower Left  ( -84.0000000,  22.9999997) ( 84d 0' 0.00"W, 23d 0' 0.00"N)
Upper Right ( -52.9999995,  45.0000000) ( 53d 0' 0.00"W, 45d 0' 0.00"N)
Lower Right ( -52.9999995,  22.9999997) ( 53d 0' 0.00"W, 23d 0' 0.00"N)
Center      ( -68.4999997,  33.9999998) ( 68d30' 0.00"W, 34d 0' 0.00"N)
Band 1 Block=2433x3 Type=Byte, ColorInterp=Palette
   Min=0.000 Max=250.000
   Minimum=0.000, Maximum=250.000, Mean=0.850, StdDev=9.542
   Metadata:
     STATISTICS_MAXIMUM=250
     STATISTICS_MEAN=0.85033548789076
     STATISTICS_MINIMUM=0
     STATISTICS_STDDEV=9.542216344358
   Color Table (RGB with 256 entries)

It has multiple entries to the color lookup table that should be treated 
as NODATA for example entries 251-255.

So I run gdal_calc.py like:

[] ~/work/oceandata/test$ gdal_calc.py -A 
A2015173174000.L0_LAC.L2_OC.tif --outfile=test.tif --calc="A*(A<251)" 
--NoDataValue=0 --co="TILED=YES"
0 .. 10 .. 20 .. 30 .. 40 .. 50 .. 60 .. 70 .. 80 .. 90 .. 100 - Done

and it create a gray-scale image with no color lut.

[] ~/work/oceandata/test$ gdalinfo test.tif
Driver: GTiff/GeoTIFF
Files: test.tif
Size is 2433, 1727
Coordinate System is:
GEOGCS["WGS 84",
     DATUM["WGS_1984",
         SPHEROID["WGS 84",6378137,298.257223563,
             AUTHORITY["EPSG","7030"]],
         AUTHORITY["EPSG","6326"]],
     PRIMEM["Greenwich",0],
     UNIT["degree",0.0174532925199433],
     AUTHORITY["EPSG","4326"]]
Origin = (-84.000000000000000,45.000000000000000)
Pixel Size = (0.012741471640766,-0.012738853693008)
Metadata:
   AREA_OR_POINT=Area
Image Structure Metadata:
   INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -84.0000000,  45.0000000) ( 84d 0' 0.00"W, 45d 0' 0.00"N)
Lower Left  ( -84.0000000,  22.9999997) ( 84d 0' 0.00"W, 23d 0' 0.00"N)
Upper Right ( -52.9999995,  45.0000000) ( 53d 0' 0.00"W, 45d 0' 0.00"N)
Lower Right ( -52.9999995,  22.9999997) ( 53d 0' 0.00"W, 23d 0' 0.00"N)
Center      ( -68.4999997,  33.9999998) ( 68d30' 0.00"W, 34d 0' 0.00"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Gray
   NoData Value=0


[] ~/work/oceandata/test$ gdalinfo --version
GDAL 1.10.1, released 2013/08/26

So is there away to get gdal_calc.py to copy the LUT to the output file?
Is there a better way/different way to set these pixels to be NODATA?

Thanks,
   -Steve