2

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
2

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 = gdo.ReadAsArray(lon_offset_px,lat_offset_px,1000,1000)
    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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.