栅格数据是一个独特的空间数据类型。ArcPy中有一个名为arcpy.sa的空间分析模块，该模块将地图代数全部整合到Python环境中，从而提高了脚本运行效率
1.列出栅格要素
ListRaster函数是以Python列表的形式返回工作空间中的栅格要素，该函数语法如下：
ListRasters({wild_crad}, {raster_type})
可选参数wild_crad通过名称限制返回的结果。参数raster_type通过栅格数据的类型显示返回的结果。
例如：
import arcpy
from arcpy import env
env.workspace = "C:/rasters"
rasterlist = arcpy.ListRasters()
for raster in rasterList:
print raster
ListRasters函数可以根据参数过滤一部分数据。下面的代码就只输出了ERDAS IMAGINE格式的文件名称。
import arcpy
from arcpy import env
env.workspace = "C:/raster"
rasterlist = arcpy.ListRasters("*", "IMG")
for raster in rasterlist:
print raster
一旦获取到栅格数据，就可以将相关函数应用于这些栅格数据。
2.描述栅格数据
Describe函数将返回指定数据元素的各种属性。针对不同的栅格元素，还会有一些具有针对性的属性。
三种不同类型的栅格元素如下：
（1）栅格数据集–一种栅格数据模型，可以存储在磁盘或者地理数据库中。又多种存储格式，包括TIFF、JPEG、IMAGINE、Ersi GRID和MrSID。可以是单波段，也可是多波段组成。
（2）栅格波段–栅格数据集中的一个图层，代表电磁光普某个范围内或波段内的值。
（3）栅格目录–以表格形式定义的栅格数据集的集合，目录中的每条记录表示一个栅格数据集。通常用于显示相邻、完全重叠的栅格数据集。
2.1栅格数据集
上述3中栅格元素具有不同的属性。例如，数据格式（TIFF，JPRG等）是栅格数据集的属性，栅格单元的大小是栅格波段的属性。可以通过dataType属性来确定栅格元素的类型。
import arcpy
from arcpy import env
env.workspace = "C:/raster"
raster = "landcover.tif"
desc = arcpy.Describe(raster)
print desc.dataType
上面的例子使用了TIFF格式的栅格数据。dataType属性的返回值为RasterDataset。栅格数据集的属性包含了一下内容：
2.2栅格波段
大多数与栅格相关属性只能通过栅格波段获取。例如分辨率。需要属性都是栅格波段所特有的，包括以下内容：
对于单波段的栅格数据集，并不需要指定波段（毕竟它只有一个波段）。因此，在描述栅格数据集时可以直接访问其属性。
import arcpy
from arcpy import env
env.workspace = "C:/Data"
rasterband = "landcover.tif"
desc = arcpy.Describe(raster)
print desc.meanCellHeight
print desc.meanCellWidth
print desc.pixelType
度量单位的类型并不包含在这些属性中，它是Spatial Reference的属性。例如下面的代码获取了一个空间参考的名称及度量单位：
spatialref = desc.spatialReference
print spatialref.name
print spatialref.linearUnitName
>>> WGS_1984_UTM_zone_50N
>>> Meter
对于多波段栅格数据来说，需要指定特定的波段，可以类似于Bnad_1、Band_2，之后才能访问栅格分辨率等属性。
import arcpy
from arcpy import env
env.workspace = "C:/raster"
rasterband = "img.tif/Band_1"
desc = arcpy.Describe(rasterband)
print desc.meanCellHeight
print desc.meanCellWidth
print desc.PixelType
注释：多波段栅格数据的某一波段有时用Layer_1表示，而不用Band_1。
3.处理栅格对象
3.1.创建栅格数据
在ArcPy中，有两种方式创建栅格对象：
（1）直接引用磁盘上已有的栅格数据
import arcpy
myraster = arcpy.Raster("C:/raster/elevation")
（2）使用地图代数语句创建栅格对象。
import arcpy
outraster = arcpy.sa.Slope(“C:/raster/evelation”)
通过这两种方式得到的栅格对象既可用于Python语句中，也可以用于地图代数的表达式中。
3.2save方法
栅格对象只有一个方法：save方法。由地图代数语句返回的栅格对象默认为临时数据，这意味着他们将在使用后被删除，save方法可以使用使栅格对象永久的保存在磁盘中，其语法如下：
save({name})
在先前的例子中，栅格对象是临时的，但是可以使用下面的代码使之变成永久的数据：
import arcpy
outrater = arcpy.sa.Slope("C:/raster/elevation")
outraster.save("C:/raster/slope")
4.Spatial Analyst模块
ArcPy中的Spatial Analyst模块arcpy.sa可以执行地图代数和其他相关操作。Spatial Analyst模块为所有的栅格处理工具提供了访问接口。通过这种方式比使用Spatial Analyst工具箱更高效。下面的代码调用Slope工具：
import arcpy
elevraster = arcpy.Raster("C:/raster/elevation")
outraster = arcpy.sa.Slope(elevraster)
5.地图代数
arcpy.sa模块处理可以访问Spatial Analyst工具接口，还提供了一系列用于执行地图代数的运算符。下面的代码使用Times工具将高程值从英尺转换成米：
import arcpy
from arcpy.sa import *
elevraster = arcpy.Raster("C:/raster/elevation")
outraster = Times(elevraster, "0.3048")
outraster.save("C:/raster/elev_m")
使用地图代数运算符可以代替Times工具。将倒数第二行代码改为：
outraster = elevraster * 0.3048"
当某些地图代数操作符不能直接用在Python表达式中，可以使用Raster Calculator工具：
RasterCalculator(expression, output_raster)
arcpy.sa模块中可用的地图代数操作符，可以分为四类：算数、按位、布尔、关系.
6.ApplyEnvironment函数
Spatial Analyst工具处理地理处理外，还有一个ApplyEnvironment函数。语法如下：
ApplyEnvironment(in_raster)
这个函数可以改变输入数据的范围、分辨率，或者为数据添加分析掩膜。下面的代码使用ApplyEnvironment函数将输出数据分辨率改为30，并将一个已有的shapefile设置为分析掩膜。
import arcpy
from arcpy import env
from arcpy.sa import *
elevfeet = arcpy.Raster("C:/raster/elevation")
elevmeter = elevfeet * 0.3048
env.cellsize = 30
env.mask = "C:/raster/studyarea.shp"
elevrasterclip = ApplyEnvironment(elevmeter)
elevrasterclip.save("C:/raster/dem30_m")
并非所有的环境设置参数都适用于ApplyEnvironment函数，而是是仅限于Cell Size、Current Workspace、Extent、Mask、Output Coordinate System、scratch Workspace和snap Raster。这些都是处理栅格数据最常用的环境参数。
7.arcpy.as模块中的类
arcpy.sa模块中也包含很多用于定义栅格工具参数的类。这些类以简洁的方式描述参数，避免了一些复杂的字符串。
7.1 remap类
Reclassification可以根据重分类表更改栅格中的值。
Reclassify(in_raster, reclass_field, remap, {missing_values})
remap是一个remap对象。根据重分类对象的性质，有两种不同的Remap类：
1.RemapValue–以单个输入值作为重分类项。
2.RemapRange–以输入值范围作为重分类项。
7.1.1RemapValue
语法：
remapValue(remapTable)
一个remapTable对象定义了一个Python列表。每个列表中包含了栅格数据的旧值和新值。在RemapValue对象中定义一个重分类的语法如下：
[[oldValue1, newValue1], [oldValue2, newValue2],…]
下面的代码以土地利用数据为例，说明了RemapValue对象的用法：
import arcpy
from arcpy import env
from arcpy.sa import *
env.workspace = "C:/raster"
myremap = remapValue(["Barren", 1], ["Mixed Forest", 4], ["Coniferous Fores", 0], ["Cropland", 2], ["Grassland", 3], ["water", 0])
outreclass = Reclassify("landuse", "S-VALUE", myremap)
outreclass.save("C:/raster/lu_recl")
7.1.2 RemapRange
RemapRange对象使用方法类似于RemapValue的使用方法。但是RemapRange处理的对象是数值范围，在RemapRange对象中定义一个重分类表的语法如下：
[[startValue, endValue, newValue], ...]
下面的代码以高程数据为例，说明了RemapRange对象的用法：
import arcpy
from arcpy import env
from arcpy.sa import *
env.workspace = ("C:/raster")
myremap = RemapRange([1, 1000, 0], [1000, 2000, 1], [2000, 3000,2], [3000, 4000, 4])
outreclass = reclassify("elevation", "TYPE", myremap)
outreclass.save("C:/raster/elev_recl")
除了Reclassify工具，在Weighted Overlay工具中也会用到重分类表。在arcpy.sa模块中还有许多其他的类，根据逻辑功能，可以将它们分布以下的几类：
2.Neighborhood类
除了Ramap类以外，Neighborhood类也被广泛使用。Neighborhood类定义了不同的形状和区域。比如Focal Statistics工具以及Neighborhood工具箱中其他的工具，都需要指定一个具体的领域。
由于存在不同类型的邻域，所以在邻域函数中需要定义一个邻域对象来设置邻域参数。例如Focal Statisticas工具的语法如下：
FocalStatisticas(in_raster, {neighbordhood}, {statistics_type}, {ignore_nodata})
有六种类型的邻域对象，每种类型对应一种参数。
例如NbrRectangle对象的语法如下：
NbrRectangle({width}, {height}, {units})
下面的代码定义了一个邻域对象，并将它作为FocalStatisticas函数的参数：
import arcpy
from arcpy import env
from arcpy.sa import *
env.workspace = "C:/raster"
mynbr = NbrRectangle(5, 5, "CELL")
outraster = FocalStatistics("landuse", mynbr, "VARIETY")
outraster.save("C:/raster/lu_var")
运行代码，将使用5*5邻域对土地覆盖类型进行统计。
8.Numpy数组
NumPyArrayToRaster和RasterToNumPyArray是常规的ArcPy函数，并不是arcpy.sa模块里的函数。这两个函数可以在栅格数据和NumPy数组之间实现转换。
数据转化的代码如下：
import arcpy, scipy
from arcpy.sa import *
inRaster = arcpy.Raster("C:/raster/myraster")
my_array = RasterToNumPyArray(inRaster)
outarray = scipy.<module>.<function>(inRaster)
outRaster = NumPyArrayToRaster(outarray)
outraster.save("C:/raster/result")
扫一扫
401
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
