python - 从光栅图像创建 numPy 数组

Question

我正在尝试将 4 波段(RGB 和 nr 红外线)栅格图像转换为 ArcMap 中的 numPy 数组。成功转换为 numpy 数组后，我想计算图像上没有数据的像素数。在 ArcMap 中检查时，这些像素颜色被标记为“无”，它们显示为黑色，但它们缺少来自波段 1、2 或 3 的红色、绿色和/或蓝色 channel 数据。我需要找到它们。

这是我目前所拥有的:

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))

Runtime error Traceback (most recent call last): File "", line 14, in File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy__init__.py", line 1814, in RasterToNumPyArray return _RasterToNumPyArray(*args, **kwargs) RuntimeError: ERROR 999998: Unexpected Error.

# 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 的有效方式。

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

这可能意味着三件事:1) 只有一个波段，2) 某些波段中的数据已设置为 0(或其他值)，3) 图像为灰度，但存储为 RGB .

您可以通过查看 numpy 数组来检查波段数:

nbands = data.shape[0]

或者直接使用arcpy:

nbands = raster.bandCount

这处理了第一种情况，但是，看起来您正在尝试检测 strip 何时没有信息，而不是它们是否存在。

如果您总是期望至少有红色、绿色和蓝色(有时是 alpha，有时不是)，那么最简单的解压缩波段有点类似于:

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

这样一来，如果有 alpha 波段，我们将忽略它，如果不存在，也无所谓。同样，这假设您的数据的形状是 nbands x nrows x ncolumns(而不是 nrows x ncolumns x nbands)。

接下来，如果我们想检查波段中的所有像素值是否都为零，请不要迭代。而是使用 numpy bool 比较。它们将(>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 数组中的每个像素。请改用矢量化表达式。