从光栅图像创建 numPy 数组

Question

我正在尝试将 4 波段（RGB 和 nr 红外）光栅图像转换为 ArcMap 中的 numPy 数组。成功转换为 numpy 数组后，我想计算图像上没有数据的像素数。在 ArcMap 中检查时，这些像素颜色被标记为“无”，它们显示为黑色，但它们缺少带 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 return _RasterToNumPyArray(*args, **kwargs) RuntimeError: ERROR 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
# 

我能够将光栅图像从代码 here 更改为 numPY 数组。

不确定 numPY 数组是如何存储的，但是当迭代它时，数据会从 y 轴开始打印出来（逐列）图像而不是 x（逐行）。

我需要切换它，这样我才能从左上角到右下角逐个像素 (RGBA) 读取数据。但是，我对 numPy 知之甚少，无法做到这一点。

我认为有问题的错误可能是由于有问题的 tiff 的大小：它在 2.5MB 时工作正常，但在 4GB 时会下降。 :(

Answer 1

您似乎在询问np.nditer。

你不想使用nditer，除非你需要低级控制。但是，您几乎不需要那种级别的控制。除非您确切知道为什么需要它，否则最好不要使用 nditer。

您拥有的是一个 3D numpy 数组。您当前正在迭代数组中的每个元素。相反，您只想遍历数组的前两个维度（宽度和高度）。

---

遍历 3D 数组

作为一个简单的例子来重现您在没有 ArcMap 的情况下看到的内容：

import numpy as np

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

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

(快速说明：我在这里使用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 布尔比较。它们会快很多 (>100x)：

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 数组中的每个像素。请改用矢量化表达式。

Answer 2

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

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

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

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

但是给了你 A = 数组（[1, 3, 5, 2, 4, 6]） 你想要的输出是

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