从光栅图像创建NumPy数组

Question

我正在尝试将一张四波段（RGB和nr红外）的栅格图像转换成ArcMap中的numPy数组。成功转换为numpy数组后，我想计算图像中没有数据的像素数。在ArcMap中检查时，这些像素的颜色标记为“None”，它们看起来是黑色的，但它们缺少第1、2或3波段的红、绿和/或蓝通道数据。我需要找到它们。

以下是我目前的进展：

import numpy
import os

myDir = "C:\\Temp\\temp"
# myFile = "4_pixel_test.tif"
myFile = "4band.tif"

# import 4band (R,G,B & nr Infrared) image
fName = os.path.join(myDir, myFile)
head, tail = os.path.split(fName)


# Convert Raster to Array, Default using LowerLeft as Origin
rasArray = arcpy.RasterToNumPyArray(fName)

# find out the number of bands in the image
nbands = rasArray.shape[0] # int
# print nbands (int)

blackCount = 0 # count the black pixels per image
th = 0 # Threhold value

# print rasArray

r, g, b, a = rasArray # not working

rCheck = numpy.any(r <= th)
gCheck = numpy.any(g <= th)
bCheck = numpy.any(b <= th)
aCheck = numpy.any(a == 0)

print rCheck
print gCheck
print bCheck
print aCheck


# show the results
if rCheck:
  print ("Black pixel (red): %s" % (tail))

elif gCheck:
  print ("Black pixel (green): %s" % (tail))

elif bCheck:
  print ("Black pixel (blue): %s" % (tail))

else:
  print ("%s okay" % (tail))

if aCheck:
  print ("Transparent pixel: %s" % (tail))

运行时错误 跟踪（最近的调用在最上面） 文件“”，第14行 在c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy__init__.py中的文件，第1814行，RasterToNumPyArray函数 返回_RasterToNumPyArray(*args, **kwargs) 运行时错误：错误999998：意外错误。

# previous code which might have incorrect numpy import
# options so I'm going with default options until I know better
# import numpy
# import os
# 
# myDir = "C:\\Temp\\temp"
# myFile = "4_pixel_test.tif"
# fName = os.path.join(myDir, myFile)
# 
# Convert Raster to Array
# rasArray = arcpy.RasterToNumPyArray(fName)
# maxVal = rasArray.max()
# minVal = rasArray.min()
# maxValpos = numpy.unravel_index(rasArray.argmax(),rasArray.shape) 
# minValpos = numpy.unravel_index(rasArray.argmin(),rasArray.shape)
# 
# desc = arcpy.Describe(fName)
# utmX = desc.extent.upperLeft.X + maxValpos[0]  
# utmY = desc.extent.upperLeft.Y - maxValpos[1]  
# 
# for pixel in numpy.nditer(rasArray):
#   # r,g,b = pixel # doesn't work  - single dimension array
#   print pixel
# 

我能够从这里的代码中将光栅图像转换为numPY数组。

不确定numPY数组如何存储，但是在迭代它时，数据从y轴开始打印，并按列向下（逐列）打印图像，而不是按x轴（逐行）打印。

我需要更改此设置，以便可以从左上到右下按像素（RGBA）读取数据。然而，我对numPy的了解不足以做到这一点。

我认为问题可能是所讨论的tiff文件大小：它在2.5MB时运行良好，但在4GB时会出现问题。 :(

- Mr Mystery Guest

---

这很可能是ArcMap的限制，如果它在>2GB时崩溃。ArcMap对于大多数版本来说都是32位应用程序（不确定最新版本...）。因此，它无法访问超过2GB的内存。使用32位版本的Python无法将>2GB的文件加载到内存中。如果您愿意，可以改用GDAL和64位版本的Python来完成此操作。 - Joe Kington

Answer 1

看起来你在问关于np.nditer的问题。除非你需要底层控制，否则不要使用nditer。然而，你几乎永远不需要那种级别的控制。最好不要使用nditer，除非你确切知道为什么需要它。

你有一个三维numpy数组。你目前正在遍历数组中的每个元素。相反，你想要遍历数组的前两个维度（宽度和高度）。

---

遍历3D数组

作为一个快速的例子，重现你在ArcMap中看到的情况：

import numpy as np

data = np.random.random((3, 10, 10))

for value in np.nditer(data):
    print value

< p >（快速提醒：我在这里使用arcpy的形状约定nbands x nrows x ncolumns。看到nrows x ncolumns x nbands也非常普遍。在这种情况下，后面部分的索引表达式将不同。）

再次强调，nditer不是您想要的东西，因此，如果您确实想要做到这一点（而不是每个r,g,b像素中的每个值），那么做法会更加易读：

import numpy as np

data = np.random.random((3, 10, 10))

for value in data.flat:
    print value

在这种情况下，两者是相同的。

---

遍历像素

接下来，您想要遍历每个像素。在这种情况下，您可以执行以下操作：

import numpy as np

data = np.random.random((3, 10, 10))

for pixel in data.reshape(3, -1).T:
    r, g, b = pixel
    print r, g, b

在这种情况下，我们将10x10x3数组临时视为100x3数组。因为numpy数组默认按照第一个轴迭代，所以这将迭代每个r,g,b元素。

如果您更喜欢，您也可以直接使用索引，但速度会慢一些：

import numpy as np

data = np.random.random((3, 10, 10))

for i, j in np.ndindex(data.shape[:-2]):
    r, g, b = data[:, i, j]
    print r, g, b

向量化，不要遍历numpy数组

总的来说，像这样逐个元素地遍历数组并不是使用numpy的有效方法。

您提到您正在尝试检测带宽何时被消除和/或设置为常数值。

有三件事可能是您所指的：1）只有一个带宽，2）某些带宽中的数据已设置为0（或另一个值），3）图像是灰度的，但存储为RGB。

您可以通过查看numpy数组来检查带宽数量：

nbands = data.shape[0]

或者直接使用arcpy：

nbands = raster.bandCount

这处理了第一种情况，然而，看起来你试图检测乐队没有信息的情况，而不是它们是否存在。

如果你总是期望至少有红色、绿色和蓝色（有时有 alpha，有时没有），那么最容易的方法是类似于解包乐队：

r, g, b = data[:3, :, :]

那样做，如果有alpha通道，我们将忽略它，如果没有，也不会有影响。再次强调，这假设你的数据形状为nbands x nrows x ncolumns（而不是nrows x ncolumns x nbands）。

接下来，如果我们想检查一个波段中所有像素值是否都为零，不要使用迭代。相反，使用numpy布尔比较，速度会快得多（>100倍）：

r, g, b = data[:3, :, :]
print np.all(r == 0) # Are all red values zero?

然而，我猜测你最常想检测的是已存储为RGB格式的灰度图像。在这种情况下，每个像素的红色、绿色和蓝色值将相等，但像素不会完全相同。你可以通过以下方式进行检查：

gray = (r == b) & (b == g)
print np.all(gray)

一般来说，你真的不想遍历numpy数组中的每个像素。相反，使用向量化表达式。

- Joe Kington

---

+1 非常好的解释。我想做什么？之前的SO问题 获取大型（2GB + TIFF）检查所述TIFF中黑色，透明像素和/或缺少的RGB通道数据的存在。然后拯救公主。 - Mr Mystery Guest

我试图将您上面的代码与我的代码结合起来，用于 i、j in numpy.ndindex(rasArray.shape[:2]): r、g、b = rasArray[i, j, :] 但是我收到了“ValueError：要展开的值太多”。 - Mr Mystery Guest

@MrMysteryGuest - 在这种情况下，我猜测了arcpy用于多波段数据的形状（目前没有Arc，所以无法测试）。相应的随机数据将是data = np.random.random((3, 10, 10))。我会更新答案以反映Arcpy用于多波段数据的形状约定。 - Joe Kington

如果你想找出带有全零的波段，不要使用循环。使用布尔比较。例如：r, g, b = data，然后 np.all(r == 0) 等。 - Joe Kington

Answer 2

假设您已经知道图像大小（n x m），并且您的1d numpy数组是A，那么这将起作用。

img2D = np.reshape(A, (m,n)).T

例如：假设您的图像数组是

img2D = array([[1, 2],
               [3, 4],
               [5, 6]])

但是你已经有了 A = array([1, 3, 5, 2, 4, 6]) 你想要的输出是

 img2D = np.reshape(A, (2, 3)).T

- Hun