[ArcPy] 1 ArcPy与栅格(Raster)

原文链接：https://www.yuque.com/higis/gisinpython/zhqgbm
GIS In Python系列文章：https://www.yuque.com/higis/gisinpython

Raster栅格对象

常用方法

打开栅格数据

Raster(inRaster) #数据类型：Raster

# 例子
r = Raster("c:/data/dem") # 绝对路径
r = Raster("19960909.img") #相对路径，当不是ArcGIS的栅格数据时，要加上后缀

保存栅格数据

rasterObj.save(path) #使用：Raster对象.save(路径字符串)

# 例子
r.save("c:/data/dem_1") # 绝对路径保存

列出工作目录下的所有栅格

arcpy.ListRasters({wild_card},{raster_type})

# 例子：列出工作空间中的Grid栅格名称
import arcpy
arcpy.env.workspace = "c:/data/DEMS"
rasters = arcpy.ListRasters("*","GRID")
for raster in rasters:
	print raster

参数	说明	数据类型
wild_card	通配符可限制返回的结果，例如匹配前面有A的文件名(“A*”)	String
raster_type	栅格格式	String

栅格转换为NumPy数组

使用例子：https://blog.csdn.net/summer_dew/article/details/78867410

转换成NumPy便于我们对像元进行操作，详细请点我

RasterToNumPyArray (in_raster, {lower_left_corner}, {ncols}, {nrows}, {nodata_to_value})

RasterToNumPyArray支持将多波段栅格直接转换成多维数组(ndarray)

1. 如果输入Raster实例基于多波段栅格，则会返回 ndarry，其中第一维的长度表示波段数。ndarray 将具有维度（波段、行、列）
2. 如果输入Raster实例基于单个栅格或多波段栅格中的特定波段，则会返回含维度（行、列）的二维数组。

转换时隐含规则：

1. 如果数组的定义（左下角以及行数和列数）超出 in_raster 的范围，则数组值将分配为 NoData
2. 如果 lower_left_corner 与像元角不重合，则会自动捕捉到最近像元角的左下角，该捕捉方法采用的规则与“捕捉栅格”环境设置中的规则相同。RasterToNumPy 函数内的捕捉操作不会与“捕捉栅格”环境设置相混淆；该函数只使用相同的交互

属性

属性	说明	数据类型	属性	说明	数据类型
bandCount	波段数量	Integer	pixelType	像素类型(U32:Unsigned 32 bit integers)	String
name	数据名称	String	spatialReference	空间参考	SpatialReference
path	完整路径和名称	String	catalogPath	全路径和名称的字符串	String
compressionType	压缩类型	String	format	数据格式	String
extent	栅格数据的范围	Extent	hasRAT	存在关联的属性表	Boolean
height	行数	Integer	width	列数	Integer
isInteger	数据具有整数类型	Boolean	isTemporary	数据是临时的	Boolean
maximum	最大值	Double	minimum	最小值	Double
mean	平均值	Double	standardDeviation	标准差	Double
uncompressedSize	磁盘大小	Double	noDataValue	在数据中NoData的值	Double

操作

【官方文档】

• 栅格计算器
• 像元统计

out_raster=in_raster1*in_raster2 #直接运算
out_raster = Slope("dem") #指定函数：参数为字符串，全路径或相对路径
out_raster = Con(in_raster >= 2, 1, 0) #栅格计算器中的函数
outStats = CellStatistics(["inraster1", "inraster2", "inraster3"], "MEAN", "DATA") #平均值

案例

名称	说明
多个栅格文件相加	创建一个相同范围，像元值都为0的栅格文件aoi_value_0，递归相加，保存

镶嵌（求多幅影像的平均值）

def mosaic(in_dir, out_fp):
    # input
    in_datas = []
    fns = os.listdir(in_dir)
    for fn in fns:
        if fn.endswith("tif"):
            fp = os.path.join(in_dir, fn)
            in_datas.append(fp)
            print "\t[data] {}".format(fp)

    if len(in_datas) == 0:
        print "\t[No data] {}".format(in_dir)
        return None
    else:
        print "\t[find data] {}".format(len(fns) )

    # copy out_data
    print "\t[new] create the out file"
    first_data = in_datas[0]
    out_raster = Raster(first_data)
    out_raster.save(out_fp)

    # Mosaic
    print "\t[doing]"
    arcpy.Mosaic_management(in_datas[1:], out_fp, "MEAN")
    print "\t[done]"

    return out_fp

求多幅影像的平均值（去除最大值、最小值）

def mosaic_nomaxmin(in_dir, out_fp):
    # input
    env.workspace = in_dir
    env.extent = "MAXOF"

    in_datas = []
    fns = os.listdir(in_dir)
    for fn in fns:
        if fn.endswith("tif"):
            in_datas.append(fn)
            print "\t[data] {}".format(fn)

    if len(in_datas) == 0:
        print "\t[No data] {}".format(in_dir)
        return None
    else:
        print "\t[find data] {}".format(len(fns))

    # max and min raster
    sum_raster = CellStatistics(in_datas, "SUM", "DATA")
    max_raster = CellStatistics(in_datas, "MAXIMUM", "DATA")
    min_raster = CellStatistics(in_datas, "MINIMUM", "DATA")

    # the number of Data
    cnt_raster = Con(IsNull(sum_raster), 0, 0)
    for data in in_datas:
        cnt_raster = Con( IsNull(data), cnt_raster, cnt_raster+1 )
    # cnt_raster.save(r"F:\out\cnt.tif")

    # compute
    print "\t[doing]"
    raster = (sum_raster - max_raster - min_raster) / ( cnt_raster - 2)
    raster.save(out_fp)
    print "\t[done]"

    return out_fp

求多幅影像的中值

def mosaic_median(in_dir, out_fp):
    # input
    env.workspace = in_dir
    env.extent = "MAXOF"

    in_datas = []
    fns = os.listdir(in_dir)
    for fn in fns:
        if fn.endswith("tif"):
            in_datas.append(fn)
            print "\t[data] {}".format(fn)

    if len(in_datas) == 0:
        print "\t[No data] {}".format(in_dir)
        return None
    else:
        print "\t[find data] {}".format(len(fn))

    # compute
    print "\t[doing]"
    raster = CellStatistics(in_datas, "MEDIAN", "DATA")
    raster.save(out_fp)
    print "\t[done]"

    return out_fp

错误集合

运行py文件出错

运行py文件出错，没有为脚本提供授权

错误：

Traceback (most recent call last):
  File "G:/workspace/python/arcpy/arcgis_running/file_add.py", line 6, in <module>
    out = Raster('19990101.img') + Raster('19990111.img')
  File "D:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\sa\Functions.py", line 4143, in Plus
    in_raster_or_constant2)
  File "D:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\sa\Utils.py", line 47, in swapper
    result = wrapper(*args, **kwargs)
  File "D:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\sa\Functions.py", line 4140, in Wrapper
    ["Plus", in_raster_or_constant1, in_raster_or_constant2])
RuntimeError

解决：
添加授权：arcpy.CheckOutExtension("spatial")

import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("spatial")
env.workspace = "E:/user/Desktop/date"
out = Raster('19990101.img') + Raster('19990111.img')

统计多个栅格的平均值 报错

【代码】

def mosaic_median(in_dir, out_fp):
    # input
    env.workspace = in_dir
    # env.extent = "MAXOF"

    in_datas = []
    fns = os.listdir(in_dir)
    for fn in fns:
        if fn.endswith("tif"):
            in_datas.append(fn)
            print "\t[data] {}".format(fn)

    if len(in_datas) == 0:
        print "\t[No data] {}".format(in_dir)
        return None
    else:
        print "\t[find data] {}".format(len(fn))

    # compute
    print "\t[doing]"
    raster = CellStatistics(in_datas, "MEDIAN", "DATA")
    raster.save(out_fp)
    print "\t[done]"

    return out_fp

【报错】

Traceback (most recent call last):
  File "D:/mycode/GISandPython/2Arcpy/ArcGISTools/2MOD04-L2/process.py", line 223, in <module>
    mosaic_all(ret, case)
  File "D:/mycode/GISandPython/2Arcpy/ArcGISTools/2MOD04-L2/process.py", line 216, in mosaic_all
    mosaic_median(month_dir, out_fp)
  File "D:/mycode/GISandPython/2Arcpy/ArcGISTools/2MOD04-L2/process.py", line 189, in mosaic_median
    raster = CellStatistics(in_datas, "MEDIAN", "DATA")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Functions.py", line 3149, in CellStatistics
    ignore_nodata)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Utils.py", line 53, in swapper
    result = wrapper(*args, **kwargs)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Functions.py", line 3145, in Wrapper
    [function] + Utils.flattenLists(in_rasters_or_constants))
RuntimeError: <exception str() failed>

【错误分析】多个栅格的范围不一致导致的错误
【解决】多加一句环境配置
env.extent = "MAXOF"