[GIS] How to add color interpretation for raster bands using GDAL

gdalgdal-merge

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.

Best Answer

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.