七、处理栅格数据【ArcGIS Python系列】

修改时间：2023年11月27日

栅格(rasters)是一种独特的空间数据类型，许多地图处理工具都是专门为栅格数据结构而设计的。可视化、处理、分析和管理影像和栅格数据是 GIS 的重要组成部分。 此处所说的栅格数据分为影像(imagery)和栅格(raster)数据。

1.影像和栅格数据简介

地理信息系统中的栅格数据将地理空间区域划分为规则的网格单元，每个单元称为一个像元(Pixel)，并为每个像元分配一个值来表示该区域的某种属性或特征。

栅格数据通常用于表示连续分布的地理现象，例如高程、温度、降雨量等。每个像元的值可以是数字或分类数据，用于描述该区域的属性或特征的某种度量。

以下是一个示例栅格数据的图示：

以栅格格式存储的数据通过以下方式表示实际现象：

• 连续数据表示光谱数据（例如卫星、航空和无人机影像）及物理和环境数据（例如高程和温度）。
• 专题数据（亦称离散数据）主要用于表示土地利用、土壤类型等属性。
• 图片包括扫描地图或工程图和建筑物照片。

栅格数据在GIS中具有广泛的应用，例如地形分析、遥感图像处理、环境建模等。它们可以用于分析、可视化和模拟地理现象，帮助我们理解和处理地理空间数据。

1）影像数据的地理属性

影像数据集的四个地理属性如下所示：

• 坐标系
• 参考坐标或 x、y 位置（通常在影像左上角或左下角）
• 像素大小
• 行计数和列计数

TODO：结合numpy处理影像数据。

2）栅格波段

一些影像具有单个波段或图层（单个特征的量度）的数据，而另一些影像具有多波段。 波段由单个像素值矩阵表示，而具有多个波段的栅格包含多个在空间上重合的表示同一空间区域的像素值矩阵。 数字高程模型 (DEM) 即是一个单波段栅格数据集的示例。 DEM 中的每个像素只包含一个表示表面高程的值。 还有一种称为全色影像或灰度影像的单波段正射影像。

卫星影像通常包含表示不同波长的多个波段，即从电磁波谱的紫外部分到可见光部分、红外部分和短波红外部分。例如，Landsat-9 影像的数据采集自电磁波谱的 11 个波段。 波段 1–7 表示来自可见光区、近红外区和中红外区的数据。 波段 6 从热红外区采集数据。 另一个多波段影像的示例是自然色正射影像，该影像包含分别表示红光、绿光和蓝光的三个波段。

可以使用栅格波段进行多种分析，以便观察和度量人眼可见和不可见的现象。 不同的材料和要素会反射和吸收电磁波谱不同部分的能量，这称为其光谱特征或光谱图。

ArcGIS Pro 中提供多种基于影像中栅格波段的预设波段组合。 “自然色彩”和“彩色红外”是常见的波段组合。 对于具有其他栅格波段（包括短波红外和热波段）的卫星传感器，可使用其他波段组合。ArcPy中影像模块同样也提供了波段合成CompositeBand (raster)函数来创建栅格对象，举例：

compband_raster = arcpy.ia.CompositeBand(["Band1.TIF", "Band2.TIF", "Band3.TIF"])

3）波段索引

传感器旨在记录电磁波谱特定部分的能量，这些部分对应于人们通常感兴趣的要素（例如水、人为要素和植被）的光谱图。 例如，许多传感器记录近红外能量 (750 - 1,000 nm)，这对于分析和监测植被类型、相对健康状况、环境压力和其他物理特性和现象非常重要。 其他波段和波段组合非常适合识别和量化各种应用的土地利用和土地覆被类别。

可使用不同的栅格波段组合创建可视化，从而能够对影像执行不同的分析。 每个可视化样式可以提供影像的不同视图，可以在其中观察和度量特定现象。 例如，在研究植被时，使用红外波段的栅格波段组合可显示健康植被，而短波红外波段可能更适合其他地质研究。 可以针对自定义可视化创建影像中任意三个波段的组合。 对于某些包含 alpha 波段的影像，第四个波段可用于表示透明度，但它在影像的图例中不可见。

遥感科学的基础是处理多光谱影像波段，以提取关于要素和现象的数据和信息。 ArcGIS Pro 中的影像处理工具、函数和功能以遥感和摄影测量概念为基础，理解和管理影像波段是影像分析和可视分析的核心。 多光谱波段可以使用算术运算进行关联、组合和处理，从而派生出特定类型的要素。 一些标准的算法和处理进程已为用户所熟知，称为指数。 这些指数按照应用进行分类，例如植被和土壤、水体、地质和景观。

常用的指数有归一化差值植被指数 (NDVI)，增强型植被指数 (EVI)，归一化差值含水指数(NDMI)。

完整的指数表格：

Arcgis Pro程序中的栅格函数

在Arcgis Pro程序中，可以通过在菜单栏的影像中找到栅格函数，单击并访问会打开栅格函数窗口。

关于栅格和影像的更多基础概念可以浏览ArcGIS 中的影像和遥感，本次重点聊聊ArcPy在处理栅格数据时的作用。

arcpy中的栅格函数则通过影像分析模块和空间分析模块提供：

Image Analyst、Spatial Analyst都提供计算指数的函数，我们拿计算 Landsat 8 影像的增强型植被指数 (EVI)为例：

增强型植被指数 (EVI) 方法是一个经过优化的植被指数，可以解释大气影响和植被背景信号。 它类似于 NDVI，但是对于背景和大气噪音不是很敏感，在查看绿色植被非常密集的区域时，颜色也不如 NDVI 那么深。

import arcpy

EVI_raster = arcpy.sa.EVI("Landsat8.tif", 5, 4, 2)

# 也可以使用
EVI_raster = arcpy.ia.EVI("Landsat8.tif", 5, 4, 2)

提示：栅格函数的输出是新的虚拟栅格图层，不是保存在磁盘上的新栅格数据集。

提示：栅格函数的使用Python脚本可以在C\Program Files\ArcGIS\Pro\Resources\Raster\Functions\System查看。

2. Arcpy中的栅格数据

ArcGIS Pro具有许多用于管理、分析、可视化和共享栅格数据的功能。通过ArcPy在Python脚本中使用这些功能有多种方法。这些办法包括：

1. 使用ArcPy函数可以完成许多栅格和影像的数据管理任务。例如，ArcGIS Pro中的**数据管理工具箱(Data Management Toolbox)包括一个栅格工具集(Raster toolset)，该工具集包含数十种地理处理工具，用于栅格数据的基本操作，包括处理栅格数据集、使用镶嵌(mosaics)**以及管理栅格数据集的属性。在其他工具箱和工具集中可以找到用于光栅数据的其他几个工具，例如，用于在栅格格式和其他投影格式之间进行数据转换。
2. 空间分析（Spatial Analyst）和影像分析（Image Analyst）扩展模块包括许多专用地理处理工具，用于处理影像和栅格数据。
3. 栅格数据集可以转换为NumPy数组以用于其他Python包。‘

Spatial Analyst和Image Analysis的功能差异

Spatial Analyst和Image Analyst模块的功能基本重叠。主要差异：

1. 功能范围：
   • Spatial Analyst：Spatial Analyst 模块提供了一系列用于空间分析的工具和函数。它包含了许多常用的栅格分析功能，如地形分析、栅格计算、栅格重分类、栅格统计等。Spatial Analyst 模块适用于各种栅格数据分析和建模任务。
   • Image Analysis：Image Analysis 模块是 ArcGIS 中专门用于遥感图像处理和分析的模块。它提供了一系列用于遥感图像处理的工具和函数，包括图像分类、变换、镶嵌、索引计算等。Image Analysis 模块适用于处理和分析遥感图像数据。
2. 数据类型支持：
   • Spatial Analyst：Spatial Analyst 模块支持处理各种栅格数据类型，包括数字高程模型（DEM）、栅格图像、栅格数据集等。它可以处理不同分辨率、不同像素类型的栅格数据。
   • Image Analysis：Image Analysis 模块专注于遥感图像数据的处理。它支持常见的遥感图像格式，如 TIFF、JPEG 等，并提供了用于读取、处理和分析这些图像的功能。
3. 工具和函数：
   • Spatial Analyst：Spatial Analyst 模块提供了一系列用于空间分析的工具和函数，如缓冲区分析、提取等值线、栅格重分类等。它还包含了一些专业的地形分析工具，如坡度计算、流域分析等。
   • Image Analysis：Image Analysis 模块提供了一系列用于遥感图像处理和分析的工具和函数，如图像分类、变换（如主成分分析、拉伸等）、索引计算（如 NDVI、NDWI 等）等。它还提供了一些用于图像增强和图像显示的工具。

总的来说，Spatial Analyst 和 Image Analysis 是 ArcPy 中用于处理栅格数据的两个模块。Spatial Analyst 模块适用于各种栅格数据的空间分析和建模任务，而 Image Analysis 模块专注于遥感图像的处理和分析。它们提供了不同的功能和工具，以满足不同的栅格数据处理需求。

以下举例如何用ArcPy进行坡度处理：

from arcpy.sa import *
out_aspectslope_raster = AspectSlope("elevation.tif", 3)
out_aspectslope_raster.save("C:/arcpyExamples/outputs/aspectslope.tif")

1.使用栅格对象

ArcPy通过引用栅格数据集的 Raster 类来作为栅格对象，有以下三种方式创建：

1. 通过引用磁盘上的现有栅格

   rasObject = Raster("C:/Data/elevation")

2. 通过运行地理处理工具

   rasObject = Slope("C:/Data/elevation")

3. 通过使用地图代数语句

Raster 对象将引用栅格数据集，而如果用于地图显示中，则还可能与内容列表中的栅格图层关联在一起。 大多数情况下，栅格数据集、Raster 对象和栅格图层之间的关系保持不变，但想要有效地使用 Spatial Analyst 地图代数，了解这些关系至关重要。

2.栅格对象的属性

栅格对象具有许多属性，包括 bandCount (波段数)、 compressionType(波段名称) 、 format (块大小)、 height 、 width 、 pixelType 、 spatialReference 等。这些属性大多是只读的。

以上属性也可以通过 Describe() 和 da.Describe() 函数获得。

import arcpy
from arcpy.sa import *
arcpy.env.workspace = "C:/Raster/Study.gdb"
ras = Raster("elevation")
print(ras.catalogPath) # 栅格数据集的完整路径和名称
print(ras.compressionType) # 压缩类型
print(ras.format) # 栅格格式
print(ras.pixelType) 
print(ras.spatialReference.name)
>>>返回值：
C:\Raster\Study.gdb\elevation
LZ77
FGDBR
S16
NAD_1983_Transverse_Mercator

一维的栅格对象只有一种方法（函数）：save()多维栅格对象具有获取多维数据。属性、名称和值的其他方法。

3.列出栅格

以上示例都是针对单个栅格数据集使用的。许多情况下我们需要处理多个栅格数据集。和之前列出空间数据一样，在处理栅格数据集时，将用ListRasters()，用法大同小异。用一个实例回忆一下：

import arcpy
arcpy.env.workspace = "C:/Raster"
rasterlist = arcpy.ListRasters("*", "IMG")
for raster in rasterlist:    
  print(raster)

4.描述栅格属性

可以使用 Describe() 和 da.Describe() 函数描述栅格数据集。

da.Describe() 后者返回字典。

import arcpy
raster = "A039697_20230128T024705"
desc = arcpy.Describe(raster)
print(desc.dataType)
>>> RasterDataset

栅格数据和前面栅格类的有相同的属性，也有不同的。有些复杂，略过。

下面的示例代码将打印单波段和多波段栅格的行数和列单元格数。对于多波段栅格，将打印每个波段的特性。我们以上海哨兵卫星图为例：

import arcpy
arcpy.env.workspace = r'./resource/data3/上海哨兵卫星.gdb'

rdesc = arcpy.da.Describe(raster)
bandcount = rdesc["bandCount"] # 波段数量

if bandcount == 1:
	print("Raster name: " + str(rdesc["baseName"]))    		
	print("No. rows: " + str(rdesc["height"]))
	print("No. columns: " + str(rdesc["width"]))
if bandcount > 1:
	counter = 1    
	while counter <= bandcount:
		band = "Band_" + str(counter)
		bdesc = arcpy.da.Describe(raster + "/" + band)
		print("Raster name: " + str(rdesc["baseName"]))
		print("Band name: " + str(bdesc["baseName"]))
		print("No. rows: " + str(bdesc["height"]))
		print("No. columns: " + str(bdesc["width"]))
		counter += 1

5.地图代数

地图代数是一种定义语法的语言，通过应用运算符创建新栅格图层来组合栅格图层。 arcpy.sa 和 arcpy.ia 模块将**地图代数(map algebra)**直接集成到Python环境中。这两个模块都提供对映射代数运算符的支持。这些运算符与栅格计算器等地理处理工具中使用的运算符相同，例如加法( + )、乘法( * )和布尔和( & )。

相当于栅格计算器。

6.其他工具

1）重分类

Reclassify

重新分类工具的语法为

Reclassify(in_raster, reclass_field, remap, {missing_values})

remap直接将重映射表的所有值作为参数键入会很复杂，因为这个表可能有许多不同的条目。相反，使用 remap 参数，以便重新映射表可以作为单独的对象创建。语法是RemapValue(remapTable)，此重映射表是一个python列表。

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

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值可以不提供

下面的代码演示了如何使用重映射对象对高程栅格进行重新分类：

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

第一个范围（1-1000）的结束值与第二个范围（1000-2000）的开始值相同，除“重分类”工具外，“加权叠加”工具还使用重映射表。

示例：下载哨兵2卫星地图并且计算NDVI

参考：哨兵Sentinel系列卫星介绍与下载教程

1.介绍

哨兵系列卫星是欧盟发射的一系列观测地球的卫星，我们此次使用的是哨兵2（Sentinel-2）由 2 颗卫星组成：2015 年发射的的 Sentinel 2A和 2017 年推出的 Sentinel 2b。

关于整个哨兵Sentinel系列卫星和下载介绍不做过多介绍，总结起来建议去欧空局官网，更新及时，可以可视化操作，只不过需要挂梯子。

本教材提供了上海地区的卫星图:S2A_MSIL1C_20230528T023531_N0509_R089_T51RUQ_20230528T045318.zip,请微信关注“renhai-lab”之后发送“哨兵卫星20230528”获取下载链接。

下载后直接解压到当前文件夹即可。

我们可以看看整个文件夹（.SAFE）的结构，下载好的文件在IMG_DATA文件夹中储存着分波段的遥感影像，为.jp2格式：

2.哨兵2卫星波段组合

我们重点来看哨兵2号卫星的波段组合：

1. Sentinel-2携带多光谱成像仪 (MSI)。该传感器提供13个光谱带，像素大小从10 米到60 米不等。

   • 蓝色 (B2)、绿色 (B3)、红色 (B4) 和近红外 (B8) 通道的分辨率为 10 米。

   • 红边（B5, red edge）、近红外NIR（B6、B7和B8A）和短波红外SWIR（B11和B12）的地面采样距离为20米。

   • 沿海气溶胶 (B1, coastal aerosol) 和卷云带 (B10, cirrus band) 的像素大小为 60 米。

在arcgis中我按照以下顺序来组合波段：

我们按如下规则对常用栅格波段进行栅格合成，可以全部合成也可以单独根据需求合成栅格

我们可以使用影像模块下的CompositeBand (raster)来合成栅格，只需要将栅格按指定顺序放入就能够生成各种专题图，对于复杂的合成，比如归一化差值植被指数 (NDVI)我们可以使用栅格计算器来操作，也可以用栅格函数来实现。

我们看看哪些常用的图的合成方法和效果：

自然色 Natural Color（B4、B3、B2）

自然色带组合使用红色 (B4)、绿色 (B3) 和蓝色 (B2) 通道。其目的是像我们的眼睛看世界一样显示图像。就像我们看到的那样，健康的植被是绿色的。其次，城市特征通常呈现白色和灰色。最后，水的深蓝色取决于水的清洁程度。

最小化的自然色栅格数据就是红绿蓝三个波段了，也是就可见波段，我们使用B04、B03、B02的波段顺序作为列表传入函数arcpy.ia.CompositeBand中：

# 定义工作空间
arcpy.env.workspace = r"resource\data4\S2A_MSIL1C_20230528T023531_N0509_R089_T51RUQ_20230528T045318.SAFE\GRANULE\L1C_T51RUQ_A041413_20230528T024830\IMG_DATA" 
# 使用波段合成函数
compband_raster = arcpy.ia.CompositeBand(['T51RUQ_20230528T023531_B04.jp2', 'T51RUQ_20230528T023531_B03.jp2', 'T51RUQ_20230528T023531_B02.jp2'])
# 可以使用save方法保存
compband_raster.save('resource/data4/上海哨兵2卫星.gdb/natural_raster')

彩色红外 Color Infrared（B8、B4、B3）

彩色红外波段组合用于强调健康和不健康的植被。通过使用近红外（B8）波段（用于反应叶绿素清空）。所以在彩色红外图像中，较茂密的植被呈红色。但城市地区是白色，水呈现充满活力的蓝色。

短波红外线 Short-Wave Infrared（B12、B8A、B4）

短波红外波段组合使用 SWIR (B12)、NIR (B8A) 和红光 (B4)。该合成图显示了各种绿色色调的植被。一般来说，较深的绿色表示植被较茂密。但棕色表示裸土和建筑区域。

农业 Agriculture（B11、B8、B2）

农业波段组合使用 SWIR-1 (B11)、近红外 (B8) 和蓝色 (B2)。由于它使用短波和近红外，因此主要用于监测农作物的健康状况。这两个波段特别擅长突出显示为深绿色的茂密植被。

地质学 Geology （B12、B11、B2）

地质波段组合是寻找地质特征的巧妙应用程序。这包括断层、岩性和地质构造。通过利用 SWIR-2 (B12)、SWIR-1 (B11) 和蓝色 (B2) 波段，地质学家倾向于使用这种前哨波段组合进行分析。

NDVI 植被指数

由于近红外光（植被强烈反射）和红光（植被吸收），植被指数有利于量化植被数量。归一化植被指数差异的公式为(B8-B4)/(B8+B4)。高值表示茂密的树冠，低值或负值表示城市和水景。

计算的黑白的NDVI（-1到1）：

我们对其应用名叫“NDVI3"的颜色映射，看起来更能分辨出哪里的植被更茂密：

NDVI3 - 用于可视化植被的色彩映射表。接近于零的值为蓝色。然后随着植被指数的由低到高，颜色也从红色逐渐变为橙色，然后变为绿色。

归一化差值水指数 NDWI （B03-B08）/（B03+B08）

NDWI使用绿色和近红外（NIR）波段。此指数的公式为： NDWI = (Green - NIR) / (Green + NIR)

水分指数 NDMI Moisture Index （B8A-B11）/（B8A+B11）

水分指数是发现植物水分胁迫的理想指标。它使用短波和近红外线来产生水分含量指数。一般来说，湿润的植被具有较高的值。但较低的水分指数值表明植物处于水分不足的压力下。

以上为栅格数据的基础内容，这里提供了一个实际案例：使用 Python 检测洪水影像的区域。

3.拓展阅读:使用numpy数组

Arcpy有两个常规的函数：NumPyArrayToRaster() 和 RasterToNumPyArray()。

可以编写一个脚本，将光栅转换为NumPy数组，然后从SciPy包调用专用函数。类似于：

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

除了支持专门的功能外，NumPy数组还用于更简单的任务以提高性能。NumPy数组被读入内存，处理速度非常快。另外，numpy还提供许多其他的属性，比如中位数、方差。例如，以下脚本在将栅格转换为NumPy数组后确定栅格的基本描述性统计：

import arcpy
import numpy
raster = "C:/Raster/elevation"
# 通过numpy读取栅格属性，通常会更快
array = arcpy.RasterToNumPyArray(raster)
print(array.min())
print(array.max())
print(array.mean())
print(array.std())
# 等同于栅格对象的属性：
dem = arcpy.Raster("C:/Raster/elevation")
print(dem.minimum)
print(dem.maximum)
print(dem.mean)
print(dem.standardDeviation)

numpy也能够直接读取栅格数据集和其他数据

在某些情况下，您可以直接将栅格数据集和其他数据集读入NumPy，而无需依赖ArcGIS转换函数。以ASCII格式的高程数据集为例。此文本文件遵循特定类型的格式。ASCII是在不同应用程序和平台之间传输数据的常用格式。这种类型的文件可以使用ArcGIS Pro中的ASCII转栅格工具转换为栅格数据集，但也可以使用 numpy.loadtext()函数直接读取到NumPy数组中。当在文本编辑器中打开时，表示光栅数据的典型ASCII文件看起来如图所示。

文章索引：

【ArcGIS Python系列】系列笔记为学习ArcGIS Pro和Arcpy过程中的总结，记下来方便回看，最新版本会优先发布在我的博客和GITHUB。

【ArcGIS Python系列】教程部分：

• 一、Arcpy介绍和安装【ArcGIS Python系列】
• 二、ArcGIS Pro和ArcMap的区别【ArcGIS Python系列】
• 三、Arcpy基础【ArcGIS Python系列】
• 四、探索空间数据【ArcGIS Python系列】
• 五、处理地理数据异常【ArcGIS Python系列】
• 六、处理几何数据【ArcGIS Python系列】
• 七、处理栅格数据【ArcGIS Python系列】
• 八、制图模块【ArcGIS Python系列】
• 九、自定义工具箱【ArcGIS Python系列】
• 十、ArcGIS_Pro常见问题【ArcGIS Python系列】
• 利用ArcGIS Python批量处理地理数据的坐标系
• 使用ArcGIS Python检测洪水影像的区域
• 利用ArcGIS_Python制作考虑路况的交通等时圈
• 利用ArcGIS Pro制作弧线OD图【ArcGIS Python系列】
• 使用ArcGIS Pro对卫星图进行建筑轮廓识别和车辆检测
• ArcGIS_Pro官方课程整理
• 持续更新…

【ArcGIS Python系列】jupyter notebook：

• 4.2.3-arcpy基础(代码练习).ipynb
• 4.2.4-探索空间数据(代码练习).ipynb
• 4.2.5-示例1：使用Arcpy进行GIS人口空间分布数据探索.ipynb
• 4.2.7-处理几何数据代码练习和示例2.ipynb
• 4.2.8-栅格数据(代码练习).ipynb
• 4.2.12-实操1-如何将ArcGIS Pro中的统一坐标系工作转换为Python自动化流程.ipynb
• 4.2.13-实操2-使用 Python 对图像中的洪水进行分类.ipynb
• 4.2.14-实操3-利用ArcGIS_Python制作考虑路况的交通等时圈.ipynb

如果你觉得本系列文章有用，欢迎关注博客，点赞和收藏，也欢迎在评论区讨论：

• 我的博客
• 我的GITHUB
• 我的GITEE
• 微信公众号: renhai-lab
• 我的知乎
• RSS订阅说明