# Merging multiple 16-bit image bands to create a True Color Tiff

I'm trying to merge three Tiff band images in python, red, green and blue into a True color image. PIL does not support the file depth as they are 16-bit files with 11-bit values. I have seen it mentioned to use a sequence of Image objects. However after attempting to research this I'm unsure how to approach it.

• Will QGIS read the images correctly? Are they 16bits per band with only 11 value bits or 16bit colour (5bits per channel + 1)? Is each band a separate 1 band image? – Michael Stimson Nov 5 '14 at 1:12
• They are 16-bits per band with only 11-bit values with individual bands. The answer below does work with 8bit values but GDAL was able to import and convert them to the 11-bit values needed. – ArcAngel Nov 5 '14 at 13:38

Hopefully this is useful. Take note, it is written for python3 and needs specific libraries.

First you would need to read each tiff band as a numpy array using something like this:

``````import numpy as np
from osgeo import gdal
import os.path

def file_test(input_file):
return os.path.isfile(input_file)

def open_tiff(tiff_file):
print(tiff_file)
if file_test(tiff_file):
return gdal.Open(tiff_file)
else:
print("Error!")

def tiff_to_array(tiff_file, lon_offset_px=0, lat_offset_px=0):
gdo = open_tiff(tiff_file)
tiff_array = np.array(bytescale(tiff_array, low=0, high=255), dtype=np.int16)
gdo = None
return tiff_array
``````

This is used to create image array. I use the scikit-image library to do some enhancements.

``````import numpy as np
from scipy.misc import bytescale
from skimage import exposure
from matplotlib import pyplot as plt

def create_composite(red_band, green_band, blue_band):
img_dim = red.shape
img = np.zeros((img_dim[0], img_dim[1], 3), dtype=np.float)
img[:,:,0] = red
img[:,:,1] = green
img[:,:,2] = blue
p2, p98 = np.percentile(img, (2, 98))
img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))
return bytescale(img_rescale)

image = create_composite(red, green, blue)
plt.imshow(image)
``````

You could convert each band to 8-bit using scipy.misc.bytescale first before passing it to the function.