[ArcGIS] 模型构造器中使用重分类 ArcPy与Numpy 相等间隔重分类

问题

ArcGIS模型构建器中对二次生成的数据无法进行重分类

如，图中dem数据生成坡度数据，这里坡度数据为二次生成的数据，造成在之后【重分类】工具中无法对它进行重分类

原因：二次生成的数据还没有准确的value，无法设置相应间隔

解决

分割工具

可以使用【重分类】中【分割工具】：
按照相等间隔区域、相等面积或自然间断点分级法分割或重分类输入像元值的范围。

分割方法	介绍
EQUAL_INTERVAL	确定输入值的范围然后将该范围分割为指定数量的输出区域。分割后输出栅格中每个区域的输入像元值都可能与极值范围相同。
EQUAL_AREA	指定输入值将被划分为指定数量的输出区域，且每个区域的像元数相同。每个区域所代表的面积大小相等。
NATURAL_BREAKS	指定各类将基于数据中固有的自然分组。中断点将通过选择分类间隔识别，这些分类间隔可对相似值进行最恰当地分组并使各类之间的差异最大化。像元值将被划分到各个类，如果数据值中出现相对较大的跳跃性，可为这些类设置界限。

使用栅格计算器实现相等间距重分类

我们可以使用一个公式，然后通过栅格计算器来计算出相等间距重分类

变量：栅格数据raster，分成groupCnt类，编号从startingNum开始

测试：groupCnt=10,startingNum=1

1. 求间距d=(max-min)/10
2. 分类结果”reclass” = Int( “raster”/d ) + 1
3. 如果按照第二步的计算结果，会出现第11类（即图像最大值的地方会归到第11类中）。但是在【重分类】工具中【最大值】是归在最后一类（第10类中）
4. 我们只需加个判断（最大值时赋予10，其他情况还是第二步的分类公式）即可，最后的公式：”reclass” = Con(“raster”==max, 10, Int( “raster”/d ) + 1)

公式：Con(“raster”==max, startingNum+groupCnt-1, Int(“raster”/ ((max-min)/groupCnt))+1 )

获取图像最大值最小值的工具：

总的模型图：groupCnt=10，startingNum=1的情况

其中，栅格计算器公式：

结果：成功分类

使用ArcPy实现相等间隔重分类

分割工具中没有对栅格数据进行相等间隔分类的类别，这里使用ArcPy实现，然后创建自定义工具，就能解决在模型中的使用。

相关资源下载：http://download.csdn.net/download/summer_dew/10131760

功能分析

将栅格等间距分成groupCnt个组数，组数从startingNum开始
例如：将栅格分成10个组数，组数编码从1开始，也就是1-10

结果：

ArcPy代码

import arcpy,sys
from arcpy.sa import *

arcpy.CheckOutExtension("sptial") # check license
# 得到在工具箱中输入的参数
rasterPath = sys.argv[1] # the path of raster data
outPath = sys.argv[2] # the path of the resulting raster data
startingNum = int(sys.argv[3]) # the starting num of reclassing
groupCnt = int(sys.argv[4]) # the number of the group
# 读取栅格
r = Raster(rasterPath) # open the raster
# 计算
max = r.maximum # get the maximum of this raster #栅格中的最大值
min = r.minimum # get the minimum of this raster #栅格中的最小值
d = (max - min)/groupCnt # 等差数值

array = arcpy.RasterToNumPyArray(r) #栅格转成numpy，方便我们对每个像元进行处理
rowNum,colNum = array.shape #获得栅格的行数，列数
for i in range(0,rowNum): #对每个像元处理
    for j in range(0,colNum):
        value = array[i][j]
        if value==max: #最大值的像元分成最后一类
            array[i][j] = groupCnt + startingNum - 1
            continue
        cnt = groupCnt - 1
        while cnt!=-1:
            if value >= min+d*cnt:
                array[i][j] = startingNum + cnt
                break
            cnt=cnt-1

#将numpy保存为栅格
lowerLeft = arcpy.Point(r.extent.XMin, r.extent.YMin)
cellWidth = r.meanCellWidth #栅格像元宽度
cellHeight = r.meanCellHeight #栅格像元高度
newRaster = arcpy.NumPyArrayToRaster(array,lowerLeft,cellWidth,cellHeight )
newRaster.save(outPath)

创建自定义工具

1. 在工具箱中添加脚本
   

2. 脚本信息填写
   

3. 脚本路径
   

4. 4个参数设置
   输入栅格：
   
   输出栅格：
   
   起始数值：
   
   组数：
   

工具使用

将坡度数据分成10份，类别数从1开始，即1-10

工具测试

将我们【自定义工具】的结果【reclass_test】与【重分类】工具的结果【reclass_right】进行相减

结果：
获得全0的图像，说明两幅图相同，工具正确。