10 栅格的使用（3）

最新推荐文章于 2022-07-25 16:14:52 发布

GISideas

最新推荐文章于 2022-07-25 16:14:52 发布

阅读量1.2k

点赞数

分类专栏： Arcgis Pro脚本学习文章标签： python 开发语言

原文链接：https://learngis.maps.arcgis.com/home/group.html?id=8ad8348a2ce4477c8dbb4d4bd8aa18c8&view=list#content

版权

Arcgis Pro脚本学习专栏收录该内容

55 篇文章 71 订阅

订阅专栏

10.9 使用 arcpy.sa 和 arcpy.ia 模块的类

arcpy.sa 和 arcpy.ia 模块包含几个主要用于定义栅格工具参数的类。通常，这些类用作工具参数的快捷方式，否则需要更复杂的字符串值。
考虑重分类工具的示例。使用此工具，栅格像元会根据重映射表（也称为“重分类表”）赋予新值。重映射表显示每个旧值如何映射到新值。图中的工具对话框显示了土地利用栅格作为适宜性模型的一部分被重新分类为几个新值的示例。

在这里插入图片描述
在此示例中，输入值由描述土地利用类型的字符串组成，但这些值也可以是数字。
重分类工具的语法是 Reclassify(in_raster, reclass_field, remap, {missing_values}) 将重映射表的所有值直接作为参数键入会很复杂，因为该表可以有许多不同的条目。而是使用 remap 参数，以便可以将重映射表创建为单独的对象。 remap 参数可以使用两个不同的类，具体取决于重新分类的性质：

RemapValue - 要重新分类的单个输入值的列表
RemapRange - 标识要重新分类的输入值范围的列表

RemapValue 对象的语法是:
RemapValue(remapTable)

使用 Python 列表定义 remapTable 对象，其中每个列表包含旧值和新值，类似于工具对话框中的重映射表。用于 RemapValue 对象的重映射表的语法是
[[oldValue1, newValue1], [oldValue2, newValue2], …]

以下代码说明了如何使用重映射对象对表示土地利用的栅格进行重分类：

import arcpy
from arcpy.sa import *
arcpy.env.workspace = "C:/Raster"
myremap = RemapValue([["Brush/transitional", 2],
                                         ["Water", 0], ["Barren land", 1],
                                        ["Built up", 1], ["Agriculture", 3],
                                        ["Forest", 5], ["Wetlands", 4]])
outreclass = Reclassify("landuse", "LANDUSE", myremap)
outreclass.save("lu_reclass")

除非需要将输入中的 NoData 值分配给不同的值，否则无需为 NoData 添加条目。使用重分类工具对话框时，您可以选择将在工具对话框中创建的重映射表保存到文本文件，或将先前创建的重映射表作为文本文件加载到工具对话框中。对于脚本，重映射表必须在脚本中创建为 Python 列表。

RemapRange 对象以类似的方式工作，但使用范围而不是单个值。用于 RemapRange 对象的重映射表的语法是

[[startValue, endValue, newValue], ...]

以下代码说明了如何使用重映射对象对高程栅格进行重新分类：

import arcpy
from arcpy.sa import *
arcpy.env.workspace = ("C:/Raster")
myremap = RemapRange([[1, 1000, 0], [1000, 2000, 1],
                                            [2000, 3000, 2], [3000, 4000, 3]])
outreclass = Reclassify("elevation", "VALUE", myremap)
outreclass.save("elev_recl")

请注意，第一个范围 (1–1000) 的结束值与第二个范围 (1000–2000) 的起始值相同，依此类推。当数据是连续的时，这种类型的重映射表很常见，例如高程栅格。除重分类工具外，重映射表也用于加权叠加工具。
arcpy.sa 模块中还有许多其他类。它们可以根据逻辑功能分为几个类别。表 10.2 列出了这些类别。
在这里插入图片描述
arcpy.ia 模块还包含三个用于处理像素块和栅格集合的类：PixelBlock、PixelBlockCollection 和 RasterCollection。
在 arcpy.sa 模块中使用更广泛的类中，除了已经讨论过的重映射类之外，还有邻域类，它定义了不同形状和大小的邻域。
例如，考虑焦点统计工具。该工具以及 Neighborhood 工具箱中的其他工具都需要定义特定的邻域。

在这里插入图片描述
邻里设置因邻里类型而异。例如，对于默认矩形邻域，设置包括单元格或地图单位的高度和宽度。但是，对于楔形邻域，参数包括开始角度、结束角度和单元格或地图单位的半径。

由于这些参数的可变性，邻域函数包括邻域对象。例如，焦点统计工具的语法如下：

FocalStatistics(in_raster, {neighborhood}, {statistics_type},
                        {ignore_nodata}, {percentile_value}))

有六个不同的邻域对象，每个对象都有自己独特的参数集。它们如下：

NbrAnnulus——定义一个环形或环状邻域，它是通过以地图单位或像元数指定内圆和外圆的半径来创建的。
NbrCircle - 定义一个圆形邻域，它是通过以地图单位或像元数指定半径来创建的。
NbrIrregular—定义一个不规则的邻域，由内核文件创建。
NbrRectangle—定义一个矩形邻域，它是通过以地图单位或像元数指定高度和宽度来创建的。
NbrWedge—定义一个楔形邻域，它是通过以地图单位或像元数指定半径和两个角度来创建的。
NbrWeight - 定义权重邻域，用于标识要包含在邻域内的像元位置以及用于乘以输入栅格像元值的权重。

例如，NbrRectangle 对象的语法是

NbrRectangle({width}, {height}, {units})

以下代码定义了一个邻域对象并在 FocalStatistics 函数中使用它

import arcpy
from arcpy.sa import *
arcpy.env.workspace = "C:/Raster"
mynbr = NbrRectangle(5, 5, "CELL")
outraster = FocalStatistics("landuse", mynbr, "VARIETY")
outraster.save("lu_var")

在此示例中，输出是基于 5 个像元 x 5 个像元的矩形邻域的土地覆盖变化栅格。
许多其他类别的类是专门的，仅用于一种或少数工具。例如，地形类包括九个不同的类，但这些类仅在地形转栅格工具中使用。

10.10 使用光栅单元迭代器

arcpy.sa 模块包括另一种处理栅格数据的方法，称为栅格像元迭代器 (RCI)。栅格像元迭代器，也称为栅格迭代器，可以迭代一个或多个栅格的所有像元。您可以读取和写入单个像元值，并且使用由行和列对组成的索引系统通过它们在栅格中的位置来引用像元。这种处理栅格的方法与地理处理工具不同，后者通常处理整个栅格而不显式检查每个位置的单个像元值。
使用 RCI 执行的典型任务包括为特定像元分配值、创建空栅格和处理 NoData 值。 RCI 作为 ArcPy 的 arcpy.ca 模块的一部分提供。 ArcGIS Pro 应用程序中没有直接等效功能，因此地理处理工具无法复制相同的功能。

有两种方法可以使用 RCI 的功能：(1) 使用在栅格对象上定义的隐式迭代器，以及 (2) 使用显式 RasterCellIterator 类。主要区别在于隐式迭代器仅适用于单个栅格数据集，而显式迭代器旨在用于多个栅格数据集。本节只讨论隐式迭代器。有关显式迭代器的详细信息，请参阅 ArcGIS Pro 帮助页面中的 RasterCellIterator 类文档。

使用 RCI 需要使用 Raster 类创建的栅格对象。以下代码创建了一个名为 vegras 的栅格对象：

from arcpy.sa import *
vegras = Raster("C:/Raster/vegetation.tif")

在此栅格对象上定义了一个隐式迭代器，它枚举行和列索引对。开始使用迭代器不需要额外的代码，因为它在创建栅格对象时可用。按照惯例，行通常用小写字母 i 表示，列用字母 j 表示。
栅格中任何像元的位置都可以通过其索引对 [i, j] 来引用。可以使用这个索引对来获得一个单元格的单元格值——例如，vegras[i, j]。

以下代码遍历栅格对象中的所有像元，并使用此索引表示法打印每个像元的位置及其值：

for i, j in vegras:
     print(i, j, veg[i, j])

考虑一个只有五行五列的小栅格，具有五个不同的像元值，如图所示。
在这里插入图片描述
上述示例脚本的打印输出的顶部如下所示

完整的打印输出由 25 行组成，每个单元格一行。打印出整个栅格的像元值本身没有用，但它说明了使用隐式栅格迭代器检查所有像元的值以及使用索引表示法来确定各个像元的位置。

使用 RCI 的一项常见任务是为特定栅格像元分配值并修改栅格数据集。从现有栅格数据集创建栅格对象时，其 readOnly 属性设置为 True。要修改栅格数据集，必须将此属性设置为 False。为栅格像元分配新值后，必须使用 save() 方法保存栅格数据集。完整的脚本如下：

from arcpy.sa import *
vegras = Raster("C:/Raster/vegetation.tif")
vegras.readOnly = False
for i, j in vegras:
    if vegras[i, j] == 4:
       vegras[i, j] = 2
vegras.save()

在这个脚本中，栅格迭代器用于对所有像元进行迭代，如果像元值等于4，则将其替换为2。得到的栅格如图所示。
在这里插入图片描述
这种细粒度控制的另一个示例是根据特定位置更改单元格值。以下脚本更改栅格数据集单行的像元值：

from arcpy.sa import *
vegras = Raster("C:/Raster/vegetation.tif")
vegras.readOnly = False
for i, j in vegras:
      vegras[4, j] = 5
vegras.save()

在此脚本中，第 4 行中的所有像元值都替换为值 5。生成的栅格如图所示。
在这里插入图片描述
使用栅格迭代器的另一个典型示例是使用 NoData 值。 NoData 的像元值对应于 Python 中的值 NaN。 NaN 代表“不是数字”，用于缺失值。您可以使用数学模块或 NumPy 包处理值 NaN。以下脚本使用浮点 NaN 值 math.nan 将特定像元值更改为 NoData。

import math
from arcpy.sa import *
vegras = Raster("C:/Raster/vegetation.tif")
vegras.readOnly = False
for i, j in vegras:
    if vegras[i, j] == 3:
        vegras[i, j] = math.nan
vegras.save()

在此脚本中，所有像元值 3 都将替换为 NoData 的值。生成的栅格如图所示。
在这里插入图片描述
要查询像元值是否为 NoData，可以使用 math.isnan() 函数。此函数确定一个值是否为 NaN 并返回一个布尔值。

其中一些结果也可以使用其他方法来完成，例如，使用重新分类等地理处理工具或使用 Con 语句的地图代数表达式。但是 RCI 提供了这些方法无法实现的额外细粒度控制，并且对现有栅格数据集进行了更改，而不是生成新的栅格对象。光栅单元迭代器文档中包含更高级的示例。

10.11 使用空间分析函数

arcpy.sa 和 arcpy.ia 模块包括许多用于处理栅格数据的附加功能。

这些函数中的大多数都基于 ArcGIS Pro 中的栅格函数。 ArcGIS Pro 中的栅格函数是直接将处理应用到栅格数据集的像元而不将新栅格写入磁盘的操作。此处理与地理处理工具形成对比，后者通常以新栅格数据集的形式将结果写入磁盘。使用栅格函数时，会在显示栅格时将计算应用于单元格，因此仅处理屏幕上可见的单元格。栅格函数的输出是虚拟栅格图层。
在 ArcGIS Pro 应用程序中，您可以通过单击栅格函数从影像选项卡中访问栅格函数，这将打开栅格函数窗格。
在这里插入图片描述
有许多不同的功能。一些功能等同于地理处理工具（例如，加权叠加），而另一些则是独特的（例如，Wind Chill）。每个函数都有一个对话框，看起来像一个地理处理工具对话框。

栅格函数和地理处理工具的对话框之间存在一些差异，但主要区别在于栅格函数的输出是新的虚拟栅格图层，而不是保存在磁盘上的新栅格数据集。这就是为什么栅格函数对话框显示“创建新图层”（而不是“运行”），并且没有输出参数的原因。

注意：光栅函数的详细讨论超出了本书的范围。可在 ArcGIS Pro 帮助页面上的“栅格函数”主题下找到详细信息。

栅格函数存储为栅格函数模板 (RFT)。您可以使用函数编辑器创建包含一系列栅格函数的工作流。函数编辑器是一种可视化编程语言，类似于模型构建器，但专为影像和栅格分析工作流而设计。这些工作流程可以保存为新的 RFT。栅格函数是用 Python 编写的，但无法从 ArcGIS Pro 中轻松访问代码。您可以在 ArcGIS Pro 安装文件夹中找到一些关联的 .py 文件，该文件夹通常位于 C:\Program Files\ArcGIS\Pro\Resources\Raster\Functions\System。脚本 Windchill.py 作为示例显示在图中。
在这里插入图片描述
您还可以创建自己的自定义栅格函数，这些函数是用 Python 编写的。尽管栅格函数使用 Python，但它们不使用 ArcPy。例如，前面的 Windchill.py 脚本不使用 import arcpy，而是严重依赖 NumPy 进行栅格数据处理。
尽管栅格函数不使用 ArcPy，但许多栅格函数也是 arcpy.sa 和 arcpy.ia 模块的函数。这些函数有效地复制了功能，但可以在 ArcPy 中使用选定的栅格函数。其中许多功能不需要 Spatial Analyst 或 Image Analyst 许可，但它们是 arcpy.sa 和 arcpy.ia 模块的一部分。
考虑 arcpy.sa 的 WindChill() 函数。该函数的语法如下：

WindChill(temperature_raster, wind_speed_raster, {temperature_units},
                        {wind_speed_units}, {wind_chill_units})

此函数使用温度和风速的输入栅格来使用众所周知的公式计算风寒系数。温度单位的默认值为华氏度，风速单位的默认值为英里/小时 (mph)。以下示例脚本使用默认单位创建新的 Windchill 栅格，如下所示：

from arcpy.sa import *
arcpy.env.workspace = "C:/Raster"
t = "temperature.tif"
ws = "windspeed.tif"
chill = WindChill(t, ws)
chill.save("windchill.tif")

由于 WindChill() 函数的输出是临时栅格对象，因此使用 save() 方法使其永久化。
重要的是要认识到 WindChill() 函数是 arcpy.sa 模块的函数，因此它可以在采用其他 ArcPy 功能的脚本中使用。该功能与 ArcGIS Pro 栅格函数窗格中名为 Wind Chill 的栅格函数相同。但是，您不会导入此功能，因为与地理处理工具不同，无法使用 ArcPy 调用 ArcGIS Pro 中的栅格函数。此外，脚本 WindChill.py 也没有直接使用。如果您有兴趣了解如何准确地执行计算，可以查阅此脚本，但在导入 arcpy.sa 模块时可以使用执行计算的实际代码。复制底层代码以进一步增强 ArcPy 的功能。
几十个栅格函数的功能可作为 arcpy.sa 和 arcpy.ia 模块的函数使用。这些功能分为几个类别，包括分析、外观、分类、转换、校正、数据管理、距离、水文、数学、重新分类、统计和表面。其中一些功能需要 Spatial Analyst 或 Image Analyst 的许可，但很多不需要。尽管与地理处理工具有一些重叠，但许多功能提供了独特的分析和处理功能。

10.12 使用 NumPy 数组

应该提到另外两个函数：NumPyArrayToRaster() 和 RasterToNumPyArray()。这些是常规 ArcPy 函数，而不是 arcpy.sa 或 arcpy.ia 模块的函数。这两个函数允许在栅格和 NumPy 数组之间进行转换。
注意：arcpy.da 模块中有几个函数可用于处理 NumPy 数组，但这些函数适用于要素类和表。但是，在栅格和 NumPy 之间进行转换的两个函数是常规 ArcPy 函数。

NumPy 数组旨在处理非常大的数据集。 NumPy 是一个 Python 包，用于使用 Python 进行科学计算。 NumPy 是“numerical Python”的缩写，通常发音为 NUM-pie，但有时也发音为 NUM-pee。除此之外，NumPy 提供了一个强大的 n 维数组对象。这种类型的对象可以在数据库之间移动数据。例如，SciPy 包包含许多可能对应用程序有用的算法，例如傅里叶变换、最大熵模型和多维图像处理，仅举几例。您可以编写一个脚本，将栅格转换为 NumPy 数组，然后调用 SciPy 包中的专用函数，而不是尝试在 ArcGIS Pro 中构建执行这些专用功能的工具或模型。通用脚本如下：

import arcpy
import scipy
from arcpy.sa import *
inRaster = arcpy.Raster("C:/Raster/myraster")
my_array = RasterToNumPyArray(inRaster)
outarray = scipy.<module>.<function>(my_array)
outraster = NumPyArrayToRaster(outarray)
outraster.save("C:/Raster/result")

这是一个简化的示例并引用了一个通用 SciPy 函数，但它说明了如何使用 NumPy 数组函数导出数据以使用另一个包进行处理，并将结果导入回 ArcGIS Pro 兼容格式——所有这些都在同一个 Python 脚本中.

除了支持专门的功能外，NumPy 数组还用于更简单的任务以提高性能。 NumPy 数组被读入内存，这使得处理它们的速度非常快。例如，以下脚本在将栅格转换为 NumPy 数组后确定栅格的基本描述性统计数据：

import arcpy
import numpy
raster = "C:/Raster/elevation"
array = arcpy.RasterToNumPyArray(raster)
print(array.min())
print(array.max())
print(array.mean())
print(array.std())

同样的基本统计数据也可以直接作为栅格对象的属性获得，如下

import arcpy
dem = arcpy.Raster("C:/Raster/elevation")
print(dem.minimum)
print(dem.maximum)
print(dem.mean)
print(dem.standardDeviation)

尽管两种解决方案都产生相同的统计数据，但使用 NumPy 数组通常更快。此外，NumPy 包括许多其他在 ArcPy 中没有直接等效的函数，例如中值和方差。
在某些情况下，您可以将栅格和其他数据集直接读入 NumPy，而无需依赖 ArcGIS 转换函数。考虑 ASCII 格式的高程数据集示例。此文本文件遵循特定类型的格式。 ASCII 是在不同应用程序和平台之间传输数据的常用格式。可以使用 ArcGIS Pro 中的 ASCII 转栅格工具将此类文件转换为栅格数据集，但也可以使用 numpy.loadtext() 函数将其直接读取为 NumPy 数组。在文本编辑器中打开时，代表栅格数据的典型 ASCII 文件如下图所示。
在这里插入图片描述
前六行代表标题信息，包括列数和行数、左下角像元左下角的 x,y 坐标、像元大小和 NoData 值。
读取 ASCII 文件时，前六行不应是实际高程值的一部分。您可以使用 numpy.loadtxt() 函数的 skiprows 参数跳过前六行，如下所示：

import numpy
infile = "C:/Project/dem.asc"
array = numpy.loadtxt(infile, skiprows=6)
print array.min()
print array.max()
print array.mean()
print array.std()

注意：这是将 ASCII 光栅读入 NumPy 的脚本的简化版本。如果 NoData 值出现在数据中，则需要附加代码来读取标头信息并说明它们。
能够使用 NumPy 和 SciPy 处理栅格数据集可以访问许多其他分析。

Points to remember

ArcPy includes the arcpy.sa (Spatial Analyst) and arcpy.ia (Image Analyst) modules, which provide access to many geoprocessing tools to work with imagery and raster data.

The Raster class of ArcPy is used to reference a raster dataset. Raster objects are important when working with raster datasets in Python scripting because they are widely used as the inputs and outputs of geoprocessing tools and map algebra statements.

ArcPy includes several functions to work with raster data. The Describe() and da.Describe() functions describe raster datasets and raster bands. The ListRasters() function lists rasters in a workspace.

The arcpy.sa and arcpy.ia modules integrate raster analysis and map algebra into the Python environment. In addition to providing access to all the tools in the Spatial Analyst and Image Analyst toolboxes, the modules contain several map algebra operators, which make scripting with rasters more efficient.

The arcpy.sa and arcpy.ia modules also contain several classes, which are used primarily for defining specific parameters unique to raster tools.

The Raster Cell Iterator makes it possible to iterate through all the cells of one or more rasters. Cells are referenced by their location in the raster using an indexing system consisting of row and column pairs. You can read and write raster cell values, which provide a detailed level of control not found in geoprocessing tools.

Raster functions in ArcGIS Pro are operations that apply processing directly to the cells of raster datasets without writing a new raster to disk. Many raster functions also are available as functions of arcpy.sa and arcpy.ia, which provide additional functionality for the processing and analyses of raster data. This functionality includes several specialized functions not available as geoprocessing tools.

Conversion functions are available to export a raster to a NumPy array, which makes it possible to use analysis functions from other Python packages.

ArcPy 包括 arcpy.sa (Spatial Analyst) 和 arcpy.ia (Image Analyst) 模块，它们提供对许多地理处理工具的访问以处理影像和栅格数据。

ArcPy 的 Raster 类用于引用栅格数据集。在 Python 脚本中处理栅格数据集时，栅格对象非常重要，因为它们被广泛用作地理处理工具和地图代数语句的输入和输出。

ArcPy 包含多个处理栅格数据的函数。 Describe() 和 da.Describe() 函数描述栅格数据集和栅格波段。 ListRasters() 函数列出工作空间中的栅格。arcpy.sa 和 arcpy.ia 模块将栅格分析和地图代数集成到 Python 环境中。除了提供对 Spatial Analyst 和 Image Analyst 工具箱中所有工具的访问之外，这些模块还包含多个地图代数运算符，这使得使用栅格的脚本编写更加高效。

arcpy.sa 和 arcpy.ia 模块还包含几个类，它们主要用于定义栅格工具特有的特定参数。
栅格像元迭代器可以遍历一个或多个栅格的所有像元。

使用由行和列对组成的索引系统，通过它们在栅格中的位置来引用像元。您可以读取和写入栅格像元值，从而提供地理处理工具中没有的详细控制级别。

ArcGIS Pro 中的栅格函数是直接将处理应用到栅格数据集的像元而不将新栅格写入磁盘的操作。许多栅格函数也可用作 arcpy.sa 和 arcpy.ia 的函数，它们为栅格数据的处理和分析提供附加功能。此功能包括一些不能作为地理处理工具使用的专用功能。

转换函数可用于将栅格导出到 NumPy 数组，从而可以使用其他 Python 包中的分析函数。