3-D 数组中的 Numpy .where 函数产生意外的输出

Question

我正在尝试对一些 Landsat 5 TM 卫星图像执行辐射测量和大气校正 - 一个 6 波段堆栈，输入一个 3 维 numpy 数组，排列如下：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 值，这将为您提供最小值