Loading satellite spectral bands in Ipython error

Question

I am quite new to working in the prompt for image analysis and data handling, so I realize this might very well be a very silly problem. In any occasion, I am using my tutor's code in order to radiometrically, atmospherically and topographically correct my images, which are in .tif format. I want to correct each band separately. So far the code looks like this:

    # import libraries

import numpy as np                               #numerical python

import matplotlib.pyplot as plt                  #plotting libraries
from pylab import *   
from matplotlib import cm  
from matplotlib.colors import ListedColormap

from spectral.io import envi                     #envi files library

import gdal                                      #Geospatial Data Abstraction 
                  Library

wdir = 'E:/MOMUT1/Bansko_Sat_Img/OutputII/Ipython/copies/'

def readtiff ():
   filename1 = 'band1'
   info = envi.open(wdir+'band1')
   data = info.load()
   return data, info

[Wallow_TM,infoa]=readtiff(band1.tif)  # This load a Landsat 5 Thematic Mapper (TM) image

It is at this particular point that it returns a name error, saying "band1 is not define".

Answer 1

You seem to be using spectral python's class for ENVI format files, whilst it looks like you are trying to read a single band of a Landsat image as a TIFF file.

A TIFF can easily be loaded into a numpy array using matplotlib, so your function would become

import matplotlib.pyplot as plt
def readtiff(tiff_file):
    return, plt.imread(tiff_file)

# load the data from the file to a numpy array
filename = 'band1.tif'
data = readtiff(filename)

Note that matplotlib here will just be returning an ndarray of the image data and will not include any of the image metadata.

For this sort of image processing, you ought to be using the Geospatial Data Abstraction Library (GDAL) directly. There is an excellent tutorial on loading and processing raster imagery in python.

In this case, with GDAL you can get your image as a numpy array with:

from osgeo import gdal
import numpy as np
dataset = gdal.Open('filename.tif')
band = dataset.GetRasterBand(1)
data = band.ReadAsArray()

Answer 2

Well, the problem got resolved with due help from my tutor. Turns out the array in the code I was using was arranged: "Lon-Lat-Band", whereas the images I was using were arranged like: "Band-Lon-Lat", which produced the error. Switching the places where the column was in the list fixed the problem! Thanks to the whole community for the communal effort!