Numpy .where function in a 3-D array producing unexpected output

Question

I am trying to execute an radiometric and atmoshperic correction on some Landsat 5 TM satellite images - a 6-band stack, entered into a 3 dimensional numpy array arranged as follows: band, Rows (Lat), Columns(Lon).

The issue comes from the empirical method I am trying to use - dark pixel subtraction. The calculation requires the identification in the darkest pixel of the image.

The problem is that in the way the image is arranged in the array it includes the 0 value pixels, which are the, so to say, frame of the image. I wish to exclude these pixels from the image using numpy.where:

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


#choose the path of the directory with your data
wdir = 'E:/MOMUT1/Bansko_Sat_Img/OutputII/Ipython/stack/'
filename1 = 'Stack_1986.tif'

## open tiff files
data = gdal.Open(wdir+filename1)
tiff_image = data.ReadAsArray()
print(tiff_image.shape)

#exclude one of the bands (thermal band) out of the array
tiff_new = np.zeros((6,7051,7791))
tiff_new[0:4,:,:] = tiff_image[0:4,:,:]
tiff_new[5,:,:] = tiff_image[6,:,:]
del tiff_image


###########
#Convert to TOA (top-of-atmosphere) and surface reflectance

LMAX = [169, 333, 264, 221, 30.2, 16.5] # These are the Landsat calibration coefficients per band
LMIN = [-1.52, -2.84, -1.17, -1.51, -0.37, -0.15];
ESUN = [1983, 1796, 1536, 1031, 220, 83.44]; # These are the solar irradiances per waveband

QCALMAX = 255
QCALMIN = 0

theta=59.40 # Solar zenith angle
DOY=128; # Day of acquisition
d = 1-0.01674*cos(0.9856*(DOY-4)); # This is the formula for the Sun-Earth distance in astronomical units

#create empty matrices (with zeros) to fill in in the loop bellow:
refl_TOA = np.zeros (tiff_new.shape, single)  #same shape as the tiff_new
refl_surf = np.zeros (tiff_new.shape, single)

    for i in range(5):

       im = np.squeeze(tiff_new[i,:,:])#squeezing out the bands out of the array
       im = np.where(im == 0, (NaN) , im)#excluding 0 value pixels
       L = ((LMAX[i] - LMIN[i])/(QCALMAX - QCALMIN))* im + LMIN[i] #This formula calculates radiance from the pixel values
       L = np.single(L) # For memory reason we convert this to single precision
       L1percent = (0.01 * ESUN[i] * np.cos(np.rad2deg(theta))) / (d**2 * pi) # This calculates the theoretical radiance of a dark object as 1 % of the maximum possible radiance 
       Lmini = L.min() # This command finds the darkest pixel in the image
       Lhaze = Lmini - L1percent # The difference between the theoretical 1 % radiance of a dark object and the radiance of the darkest image pixel is due to the atmosphere (this is a simplified empirical method!)
       refl_TOA[i,:,:]=(pi * L  * d**2) / (ESUN[i] * np.cos(np.rad2deg(theta))) # This is the formula for TOA reflectance
       refl_surf[i,:,:]=(pi * (L - Lhaze) * d**2) / (ESUN[i] * np.cos(np.rad2deg(theta))) # This is the formula for surface reflectance in which Lhaze is subtracted from all radiance values

    imshow(refl_surf[1,:,:])
    colorbar()
    show()

The code runs, but the output is not as it should be. Instead of the satellite image I see this image:

Something in my .where statement is not correct, as it seems to be selecting all of the pixels and giving them a NaN value, because when I try to check for pixel value by hovering over the image with my mouse courser, no numbers are shown.

Can someone help identifying what is wrong with the way I am using np.where?

Answer 1

I think your problem is actually here:

Lmini = L.min()

ndarray.min propagates NaN values, so when you multiply through by it later you turn your whole array into NaN. You want

Lmini = np.nanmin(L)

Which will give you the minimum without NaN values