2

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: visualization of one band Visualisation of 3 bands

  • One Solution is to use imshow function from the spectral module. add from spectral import * to the preamble... – user1771806 Mar 31 '14 at 15:27
3

This happens because you are crossing the matplotlib's color limits many times. Matplotlib handles 256 shades of a channel, if your image has 257 shades it will start over. In other words, the value 257 will be the same shade as the 0.

So, for example, if your image value is from 0 to 8192, you should divide the value so that it doesn't cross the 256 shade limit. For simplicity, just divide all the values by the maximum value of the image.

In your case, suppose your maximum value is 8192. Then you should do:

img = np.dstack((img1/8192,img2/8192,img3/8192))
2

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.

Output of above code

0

A solution here can be using

imread

It reads the image as such, with the bands included into a 3D array. You can, in return, call each band. More information can be found here.

You can also work with Spectral Python - SPy and/or HyperSpy. The issues with HyperSpy could be its stable versions and errors, but SPy works. Direct functions to read the hyperspectral dataset will avoid the need of excess lines of code. They provide direct functions for many applications too.

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.