Функция Numpy .where в трехмерном массиве, выдающая неожиданный вывод

Question

Я пытаюсь выполнить радиометрическую и атмосферную коррекцию на некоторых спутниковых изображениях Landsat 5 TM - 6-полосный стек, введенный в 3-мерный массив кусков, расположенный следующим образом: band, Rows (Lat), Columns(Lon),

Проблема заключается в эмпирическом методе, который я пытаюсь использовать - вычитание темных пикселей. Расчет требует идентификации в самом темном пикселе изображения.

Проблема состоит в том, что способ, которым изображение упорядочено в массиве, включает в себя пиксели с 0 значениями, которые являются, так сказать, кадром изображения. Я хочу исключить эти пиксели из изображения, используя 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()

Код выполняется, но вывод не так, как должно быть. Вместо спутникового изображения я вижу это изображение:

Что-то в моем .where утверждение неверно, так как кажется, что все пиксели выбраны и даны NaN значение, потому что, когда я пытаюсь проверить значение пикселя, наведя курсор мыши на изображение с помощью курсора мыши, цифры не отображаются.

Может кто-нибудь помочь определить, что не так с тем, как я использую np.where?

Answer 1

Я думаю, что ваша проблема на самом деле здесь:

Lmini = L.min()

ndarray.min распростран NaN значения, поэтому, когда вы умножаете на это позже, вы превращаете весь массив в NaN, Ты хочешь

Lmini = np.nanmin(L)

Что даст вам минимум без NaN ценности