从栅格图像创建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))

运行时错误Traceback(最近一次调用最后一次):文件"",第14行,在文件"c:\ program files(x86)\ arcgis\desktop10.2\arcpy\arcpy__init __.py",第1814行,在RasterToNumPyArray中返回_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
# 

我能够从这里的代码将光栅图像更改为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

Vectorize,不迭代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数组中的每个像素.改为使用矢量化表达式.