[GIS] Visualizing.tif image using Python

Question

I am trying to read a hyperspectral image in Uint16 .tif format and visualize it.
I use GDAL in Python to read the data and everything works perfectly when I use matplotlib to show one band of the data (The first image). The image however becomes completely corrupted when I want to concatenate(stack) 3 different selected bands and visualize them. Here is the code:

import numpy as np
import struct
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
from osgeo import gdal_array
from osgeo.gdalconst import *
import matplotlib.pyplot as plt

cube = gdal.Open(filename)
bnd1 = cube.GetRasterBand(29)
bnd2 = cube.GetRasterBand(19)
bnd3 = cube.GetRasterBand(9)

Now I turn each band into a ndarray:

img1 = bnd1.ReadAsArray(0,0,cube.RasterXSize, cube.RasterYSize)
img2 = bnd2.ReadAsArray(0,0,cube.RasterXSize, cube.RasterYSize)
img3 = bnd3.ReadAsArray(0,0,cube.RasterXSize, cube.RasterYSize)    

and stack them to have a 3 band image

img = np.dstack((img1,img2,img3))

f = plt.figure()
plt.imshow(img)    
plt.show()

here are the results:

Answer 1

I don't know your input data, but the reason is most likely that you haven't normalized the DN or radiance values. matplotlib.pyplot.imshow() parses RGB data only if all channels are normalized to values between 0 and 1. Also, you should convert the data to float32 or uInt8 for matplotlib.

This happened to me before, so here's a (very verbose) example to visualize what happens if your bands are not normalized -- for anyone who comes across it in the future. You'll see that the modified array plots differently than the original or the normalized array.

import matplotlib.pyplot as plt
import numpy as np

def normalize(arr):
    ''' Function to normalize an input array to 0-1 '''
    arr_min = arr.min()
    arr_max = arr.max()
    return (arr - arr_min) / (arr_max - arr_min)


# Make up some random data (5*5 arrays) between 0 and 1
r = normalize(np.random.rand(5,5))
g = normalize(np.random.rand(5,5))
b = normalize(np.random.rand(5,5))

# Modify it with different offsets and factors
rtwo = r * 200 + 12.5
gtwo = g * 100 + 5.4
btwo = b * 300 + 1.1

# Normalize it
rtwo_n = normalize(rtwo)
gtwo_n = normalize(gtwo)
btwo_n = normalize(btwo)

# Plot them
f, (ax1, ax2, ax3) = plt.subplots(1, 3)
ax1.imshow(np.dstack((r,g,b)), interpolation='nearest')
ax2.imshow(np.dstack((rtwo,gtwo,btwo)), interpolation='nearest')
ax3.imshow(np.dstack((rtwo_n, gtwo_n, btwo_n)), interpolation='nearest')

ax1.set_title("Original")
ax2.set_title("Modified")
ax3.set_title("Modified and normalized")
plt.show()

This is what the plot looks like. Note that when dealing with sattelite images, many people tend to use raw DN or radiance values, which is not normalized to any common factor. You will have to normalize your input or convert it to surface reflectance.