从光栅图像创建 numPy 数组
Create numPy array from raster image
我正在尝试将 4 波段（RGB 和 nr 红外线）栅格图像转换为 ArcMap 中的 numPy 数组。成功转换为 numpy 数组后，我想计算图像上没有数据的像素数。在 ArcMap 中检查时，这些像素颜色标记为 "None"，它们显示为黑色，但它们缺少来自波段 1,2 或 3 的红色、绿色 and/or 蓝色通道数据。我需要找到他们。
这是我目前的情况：
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 时就会下降。 :(
假设您已经知道图像大小 (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
你好像在问 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 的有效方法。
您提到您正在尝试检测条带何时被消除and/or设置为常数值。
这可能意味着三件事：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 数组中的每个像素。请改用矢量化表达式。
我正在尝试将 4 波段（RGB 和 nr 红外线）栅格图像转换为 ArcMap 中的 numPy 数组。成功转换为 numpy 数组后，我想计算图像上没有数据的像素数。在 ArcMap 中检查时，这些像素颜色标记为 "None"，它们显示为黑色，但它们缺少来自波段 1,2 或 3 的红色、绿色 and/or 蓝色通道数据。我需要找到他们。
这是我目前的情况：
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 时就会下降。 :(
假设您已经知道图像大小 (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
你好像在问 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 的有效方法。
您提到您正在尝试检测条带何时被消除and/or设置为常数值。
这可能意味着三件事：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 数组中的每个像素。请改用矢量化表达式。