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.
You can try to use gdalbuildvrt (http://www.gdal.org/gdalbuildvrt.html) instead of gdal_merge.py. As the result the xml-like file describing virtual raster created. GDAL and QGIS works with this file.
gdalbuildvrt qqv.vrt -separate qq1.tif qq3.tif
The only problem is, that only first band from all files will be added. This simply fixed by editing vrt file.
$ gdalinfo qqv.vrt
Driver: VRT/Virtual Raster
Files: qqv.vrt
/home/bishop/tmp/qq1.tif
/home/bishop/tmp/qq3.tif
Size is 735, 547
Coordinate System is:
GEOGCS["Unknown datum based upon the WGS 84 ellipsoid",
DATUM["Not_specified_based_on_WGS_84_ellipsoid",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6030"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4030"]]
Origin = (513165.703783421427943,5016826.563432835973799)
Pixel Size = (59.969933184855272,-60.001568937759295)
Corner Coordinates:
Upper Left ( 513165.704, 5016826.563) (Invalid angle,Invalid angle)
Lower Left ( 513165.704, 4984005.705) (Invalid angle,Invalid angle)
Upper Right ( 557243.605, 5016826.563) (Invalid angle,Invalid angle)
Lower Right ( 557243.605, 4984005.705) (Invalid angle,Invalid angle)
Center ( 535204.654, 5000416.134) (Invalid angle,Invalid angle)
Band 1 Block=128x128 Type=Byte, ColorInterp=Undefined
Min=0.000 Max=255.000
NoData Value=nan
Band 2 Block=128x128 Type=Byte, ColorInterp=Undefined
Min=0.000 Max=255.000
NoData Value=nan
Band 3 Block=128x128 Type=Byte, ColorInterp=Undefined
Min=0.000 Max=255.000
NoData Value=nan
Band 4 Block=128x128 Type=Int16, ColorInterp=Undefined
Raster VRT format tutorial: http://www.gdal.org/gdal_vrttut.html
Best Answer
So I figured it out!
gdal_translate -b 2 input_rgb.tiff output_green.tiffgdaldem color-relief output_green.tiff green_colortable.txt output_green_rgb.tifftaking the colortable from QGIS.I did this for both red and green bands. After that I combined them by using Mapnik!