How to add color interpretation for raster bands using GDAL?

Question

I have 3 separate 1-band GeoTIFF files. The color interpretation for band is gray. I want a 3-band RGB file. I've used gdal_merge.py

gdal_merge.py -separate file1.tif file2.tif file3.tif -o output_file.tif

But the color interpretation for the 3 bands of the resulting output_file.tif is gray, undefined, undefined. Apart from that, all pixels are 0.

If I execute pct2rgb.py, I get a RGB file:

pct2rgb.py output_file.tif output_file_rgb.tif

But of course, the pixels are 0 too. So, I have 3 questions:

• Is gdal_merge.py the right tool to combine 3 1-band files in one 3-band RGB file?
• Why am I getting undefined color interpretation for bands?
• Is pct2rgb the right tool to transform 3-band files with this color interpretation to RGB files?

UPDATE: The rasters does not have a color table. Just color interpretation: Gray.

On other hand, pixel values goes from 0 to 1023 (this is deliberate)

More data: they're rotated rasters (no north up), but all of them have the same geotransform.

UPDATE 2: I can warp the images to make them north up, construct a VRT and add ColorInterp for each band, but I still get color interpretation as gray, undefined, undefined in the output result.

The problem is I need to create a color table in, at least, the first band. I know the way to create them, but I don't know how many entries should my table have. Why are there 13 entries in the example of GDAL Raster FAQ? All the pixels have values between 0 and 1023, if helps.

UPDATE 3: Apparently, there's no way in the TIFF format to really specify the color interpretation of each band. The way GDAL builds the color interpretation when reading a TIFF file is a combination of the value of the PHOTOMETRIC and EXTRASAMPLES tag.

Reading about these tags:

• PHOTOMETRIC represents the color space of the image data. A value of 2 means that the components of a pixel value are RGB, but it assumes Byte pixels, and I have UInt16 pixels (I tried -co "PHOTOMETRIC=rgb", and got an error). So, I can't specify PHOTOMETRIC tag for the output file.

• EXTRASAMPLES specifies that each pixel has N extra components. I'm not sure about how to use this tag to create my merged file. Or if a I need it.

So, in update 2 I suggest the creation of a ColorTable, but how? In my 3 input files, pixel values go from 0 to 1023. Do I have to match them with colors? Do I have to create a ColorTable with 1024 inputs? How?

In update 3, seems that I could use some GeoTIFF tags when creating the merged file, but I'm not sure if I really can use them, or how.

Answer 1

gdal_merge.py is the correct tool to 'stack' your input images.

Assuming that your first band has a valid color table you could use:

gdal_merge.py -separate -pct -o output_file.tif file1.tif file2.tif file3.tif

Note: The command has been reformatted with -o output_file.tif before the list of inputs.

From the docs:

-pct: Grab a pseudocolor table from the first input image, and use it for the output. Merging pseudocolored images this way assumes that all input files use the same color table.

I would test your output with gdalinfo -stats to make sure it is being stacked properly.

Updated for OP

From the osgeo list, it looks like you might try a different format to check the results:

There's no way in the TIFF format to really specify the color interpretation of each band. The way GDAL builds the color interpretation when reading a TIFF file is a combination of the value of the PHOTOMETRIC and EXTRASAMPLES tag.

-Evan (the poster) knows GDAL inside and out.

Answer 2

The first way I can think of is to build a vrt, edit and translate:

gdalbuildvrt -separate output.vrt file1.tif file2.tif file3.tif

add the color interp tag where needed:

...
<VRTRasterBand dataType="Byte" band="1">
<ColorInterp>Red</ColorInterp>
<NoDataValue>255</NoDataValue>
<ComplexSource>
  <SourceFilename relativeToVRT="1">test.tif</SourceFilename>
  <SourceBand>1</SourceBand>
  <SourceProperties RasterXSize="42" RasterYSize="22" DataType="Byte" BlockXSize="42" BlockYSize="22" />
  <SrcRect xOff="0" yOff="0" xSize="42" ySize="22" />
  <DstRect xOff="0" yOff="0" xSize="42" ySize="22" />
  <NODATA>-32768</NODATA>
</ComplexSource>

and then:

gdal_translate output.vrt output.tif

Answer 3

Your input grayscale band has values from 0 to 1024, but RGB only allows values from 0 to 255 for each band.

You can use gdal_translate -scale [src_min src_max [dst_min dst_max]] to rescale your source file.

In a second run, you can stack your three bands with gdal_merge.py and apply -co "PHOTOMETRIC=rgb".