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
3 Answers 3
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
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))
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.
imshow
function from thespectral
module. addfrom spectral import *
to the preamble...