可以使用 Raster Cell Iterator (RCI) 来执行自定义栅格分析。它允许您访问单个栅格像元位置,并提供迭代控制来查询和修改栅格像元值,所有这些操作均在 Python 环境中执行。
使用栅格像元迭代器的基本操作
在本部分中,您将对可使用 RCI 执行的基本操作有一个概念性了解。您将学习如何执行以下操作:创建隐式 RCI、查询值并将值分配给栅格像元、使用元数据创建空栅格、为栅格对象创建索引以及处理 NoData。每个部分均包含代码示例,用于演示如何在 Python 窗口中使用 RCI。建议您自上而下阅读本文档以理解这些代码示例。
隐式栅格像元迭代器
可以在栅格对象上隐式调用 RCI,以迭代栅格的行索引和列索引。以下显示的代码示例将根据现有栅格数据集 myRas 创建一个名为 "myras" 的栅格对象。将在栅格对象上定义一个隐式栅格迭代器,用于在循环中对整个栅格的行和列索引对进行枚举。在迭代器内,将打印每个像元位置处的行索引 i、列索引 j 和像元值 myRas[i,j]。本示例将使用索引记法来查询给定像元位置处的像元值。
from arcpy.sa import *
myRas = Raster("myras")
for i,j in myRas:
print(i, j, myRas[i,j])
为像元分配值
借助 RCI,可以为栅格像元分配值,由此修改栅格数据集。默认情况下,根据现有栅格数据集创建的栅格对象具有 readOnly 属性,该属性设置为 True。要修改栅格数据集,应首先将 readOnly 属性设置为 False。修改栅格像元后,应在栅格对象上调用 save() 以保留这些更改。有关 readOnly 属性和保存栅格的详细信息,请参阅 Raster 对象帮助主题。
myRas.readOnly = False
for i,j in myRas:
if myRas[i,j] == 0:
myRas[i,j] = 128
myRas.save()
创建空栅格
借助 getRasterInfo() 方法,可以轻松创建一个空栅格对象,并从现有栅格对象复制元数据。将在现有栅格对象上调用 getRasterInfo() 方法,该方法将针对此栅格返回一个 rasterInfo() 对象。然后将 rasterInfo() 对象作为输入参数传递给 Raster() 类以实例化一个空栅格。空栅格对象的元数据与根据其创建空栅格对象的栅格对象的元数据相同,但所有栅格像元中均为 NoData。栅格元数据包括 bandCount、像元大小、extent、spatialReference、pixelType 和 noDataValues,以及其他属性。使用 getRasterInfo() 方法创建栅格对象时,readOnly 属性将默认设置为 False,由此可将值分配给此空栅格对象。
outRas = Raster(myRas.getRasterInfo())
您还可以创建一个空的 rasterInfo 类对象,然后在 FromJSONString() 对象上使用 rasterInfo() 方法来填充栅格元数据。随后,可以使用该对象来创建一个空栅格,如下所示。
ri = arcpy.RasterInfo()
ri_as_jsonstring = '{"bandCount":1,' \
'"extent":{"xmin":0,"ymin":0,"xmax":10,"ymax":10,"spatialReference":{"wkid":26912}},' \
'"pixelSizeX":1,' \
'"pixelSizeY":1,' \
'"pixelType":"F32"}'
ri.fromJSONString(ri_as_jsonstring)
outRas = Raster(ri)
为栅格对象创建索引
为栅格对象创建索引允许对栅格数据集中的特定像元进行读写访问。可以在栅格像元迭代器内外使用索引。
对于单波段栅格,栅格对象具有两个索引 [i, j],分别为行索引和列索引。查询和分配值至少需要两个索引。所有索引均以零为基础,从零开始。行索引的值范围定义为 0 到 n(rows) - 1,其中 n(rows) 为该栅格中的行数。同样,列索引的值范围定义为 0 到 n(cols) - 1,其中 n(cols) 为该栅格中的列数。
索引将从栅格数据集的左上角开始,其中最左上角像元的索引为 [0, 0],最右下角像元的索引为 [n(rows) - 1, n(cols) - 1]。
处理 NoData 像元
在 RCI 上下文内外使用索引记法来查询具有 NoData 的像元时,在 Python 中返回的值将为 NaN。例如,如果已知 myRas 对于行和列索引 [2, 2] 具有 NoData 值,则查询此像元值将返回 NaN。可以使用 math.isnan() 方法对其进行验证,该方法针对此像元将返回布尔值 True。
import math
math.isnan(myRas[2,2])
out: True
如果在 Python 中为像元分配 NaN,则栅格像元将变为 NoData,不考虑该栅格数据集的格式或 NoDataValues。这是将像元转换为 NoData 的正确方法。例如,在行和列索引 NaN 处将 myRas 分配给 [3, 3] 会将该像元修改为 NoData。
myRas[3,3] = math.nan
建议您不要通过将该栅格数据集的 noDataValues 直接分配给像元以将像元转换为 NoData,而应该按照上述方法执行此操作。
使用栅格像元迭代器的高级操作
在本部分中,您将对可使用 RCI 执行的高级操作有一个概念性了解。您将学习如何创建一个显式栅格像元迭代器以迭代多个栅格,以及如何使用 padding 和 skipNoData 选项。您将学习创建多波段栅格的索引并在 RCI 中使用该索引。您将学习如何处理不同的栅格像素类型以及如何处理超出范围的值。
显式栅格像元迭代器
调用 RCI 的另一种方法是使用显式 RasterCellIterator 类对象。隐式迭代器适用于对单个栅格数据集进行迭代,而显式栅格迭代器用于对多个栅格进行迭代,并针对此目的对其进行了优化。在迭代多个栅格时,可以使用输入栅格和输出栅格来执行分析,在计算中将使用输入栅格的像元值,并且在计算过程中将写入输出栅格的像元值。
以下示例演示了如何使用显式 RCI 对多个输入和输出栅格进行迭代。输入栅格 myRas1 和 myRas2 具有相同数量的行和列以及其他属性,例如像元大小、范围和空间参考。使用第一个输入栅格的 outRas1 类对象创建输出栅格 outRas2 和 rasterInfo。对于 myRas1 和 myRas2,将在每个栅格像元处执行焦点运算,以计算 3 x 3 像元邻域的平均值。在此像元位置,将为 outRas1 分配一个像元值作为平均值的最小值,并为 outRas2 分配一个像元值作为平均值的最大值。
myRas1 = Raster("myras1")
myRas2 = Raster("myras2")
outRas1 = Raster(myRas1.getRasterInfo())
outRas2 = Raster(myRas1.getRasterInfo())
with RasterCellIterator({'rasters':[myRas1, myRas2, outRas1, outRas2]}) as rci:
for i,j in rci:
meanMyRas1 = (myRas1[i-1,j-1] + myRas1[i-1, j] + myRas1[i-1, j+1] + \
myRas1[i,j-1] + myRas1[i, j] + myRas1[i, j+1] + \
myRas1[i+1,j-1] + myRas1[i+1, j] + myRas1[i+1, j+1]) / 9
meanMyRas2 = (myRas2[i-1,j-1] + myRas2[i-1, j] + myRas2[i-1, j+1] + \
myRas2[i,j-1] + myRas2[i, j] + myRas2[i, j+1] + \
myRas2[i+1,j-1] + myRas2[i+1, j] + myRas2[i+1, j+1]) / 9
outRas1[i,j] = min(meanMyRas1, meanMyRas2)
outRas2[i,j] = max(meanMyRas1, meanMyRas2)
outRas1.save()
outRas2.save()
在迭代器中使用多个栅格时,列表中的第一个栅格将定义迭代器的栅格分析环境。列表中的所有其他栅格都应在空间参考、像元大小和范围方面与第一个栅格相匹配。如果此列表中第二个栅格及后续栅格的这些属性与第一个栅格不匹配,则系统会将第一个栅格的分析环境应用于它们。建议您使用第一个栅格的 rasterInfo() 对象来创建所有输出栅格,以便这些输出栅格继承其空间参考、像元大小和范围。
此外,如果您希望使用环境设置,则在迭代器中使用 ApplyEnvironment 函数之前,应使用该函数将环境设置应用于此列表中的第一个栅格。以下示例中演示了此操作。
myRas1 = Raster("myras1")
myRas2 = Raster("myras2")
#Apply environment settings to myRas1
myRas1_env = arcpy.sa.ApplyEnvironment(myRas1)
outRas1_env = Raster(myRas1_env.getRasterInfo())
with RasterCellIterator({'rasters':[myRas1_env, myRas2, outRas1_env]}) as rci:
for i,j in rci:
meanMyRas1 = (myRas1_env[i-1,j-1] + myRas1_env[i-1, j] + myRas1_env[i-1, j+1] + \
myRas1_env[i,j-1] + myRas1_env[i, j] + myRas1_env[i, j+1] + \
myRas1_env[i+1,j-1] + myRas1_env[i+1, j] + myRas1_env[i+1, j+1]) / 9
meanMyRas2 = (myRas2[i-1,j-1] + myRas2[i-1, j] + myRas2[i-1, j+1] + \
myRas2[i,j-1] + myRas2[i, j] + myRas2[i, j+1] + \
myRas2[i+1,j-1] + myRas2[i+1, j] + myRas2[i+1, j+1]) / 9
outRas1_env[i,j] = min(meanMyRas1, meanMyRas2)
outRas1_env.save()
填充选项
显式迭代器包含用于提供填充的选项。填充不会以任何方式更改输出,但是如果在迭代像元位置时访问邻域像元值,则填充可提高性能。例如,如果邻域操作具有一个内核,该内核在迭代过程中会访问距离像元位置最多两行像元或两列像元的像元值,则将提供填充 2。
with RasterCellIterator({'rasters':[myRas1, outRas1], 'padding': 2}) as rci_padded:
for i,j in rci_padded:
outRas1[i,j] = (myRas1[i-2,j-2] + myRas1[i-2, j] + myRas1[i-2, j+2] + \
myRas1[i,j-2] + myRas1[i, j] + myRas1[i, j+2] + \
myRas1[i+2,j-2] + myRas1[i+2, j] + myRas1[i+2, j+2]) / 9
outRas1.save()
跳过 NoData 选项
显式迭代器具有完全跳过 NoData 像元的功能,而非逐个单元地对其进行处理。可以使用 skipNoData 键以提供一个栅格列表,这些栅格可确定迭代器将跳过的栅格像元。对于任何给定的像元,如果此列表中的任何栅格具有 NoData 值,则将跳过该像元。由于这些栅格中 NoData 像元的状态可能会在迭代过程中发生变化,因此,是否应用跳过将取决于迭代器正在访问的栅格中每个像元的状态。此选项的主要好处是可以提高在稀疏栅格数据集上进行迭代时的性能,因为不再需要访问许多 NoData 像元。
在以下示例中,稀疏栅格 streamRas 仅在表示河流的像元中具有值,在其他像元中均为 NoData。流向栅格 flowDirRas 在栅格范围内的任何像元中都具有值,而没有任何 NoData 像元。在使用 RCI 提取沿流像元的流向栅格值时,可以使用 skipNoData 键以使迭代器仅访问具有值的流像元,从而将访问的总像元数降至最低。
streamRas = Raster("streams")
flowDirRas = Raster("flowdir")
outRas = Raster(flowDirRas.getRasterInfo())
with RasterCellIterator({'rasters':[streamRas, flowDirRas, outRas],'skipNoData':[flowDirRas, streamRas]}) as rci_skip:
for i,j in rci_skip:
outRas[i,j] = flowDirRas[i,j]
outRas.save()
创建多波段栅格的索引
对于多波段栅格,可以提供两个或三个索引。两者之间的行为差异将在下文中进行讨论。
如果提供两个索引,则其将解释为行和列索引,并且将针对指定的 [row, column] 索引从每个波段返回一组值。例如,可以使用两个索引 [2,2] 查询具有三个波段的多波段栅格“mbras”。这些索引将解释为行和列索引,并且将返回三个波段中每个波段的值作为一组值。
multibandRas = Raster("mbras")
multibandRas[2,2]
out: (4,5,6)
如果指定三个索引,则其将解释为 [band, row, column] 索引,并且将返回该波段、行像元位置和列像元位置的像元值。波段索引从零开始,并且其多波段栅格的值范围由 0 到 n(bands) - 1 确定,其中 n(bands) 为波段的总数。
multibandRas[0,2,2]
out: 4
处理超出范围的值
所有栅格都具有像素类型和相应的位深度。有关位深度的完整列表和每个像元可包含的值范围,请参阅栅格数据集像素的位深度容量。如果为栅格像元分配的值超出其位深度的可接受值范围,则将出现超出范围的情况。在上述情况下,将为像元分配 NoData。例如,栅格对象 myRas1 的像素类型为 'U8'(无符号 8 位整数)。此栅格可以包含 0-255 范围内的整数值。如果将值 260 分配给 myRas1 中的像元,则其将分配 NoData。完成此分配后,查询此像元的值将返回 NaN。
myRas1 = Raster("myras1")
myRas1Info = myRas1.getRasterInfo()
pixelType_MyRas1 = myRas1Info.getPixelType()
print(pixelType_MyRas1)
out: 'U8'
myRas1[3,3] = 260
math.isnan(myRas1[3,3])
out: True