python提取斜坡结构

python提取斜坡结构

介绍

开始是在帮师妹处理某个试验流程中发现需要进行斜坡结构的提取，后面百度找了教程一步一步的做，发现挺麻烦的，所有写了一段代码，所以里面文件夹名字可能emmmm，不重要，这些步骤主要还是python二开，比较简单，如果有什么写的不好的请大家多多包涵。这是我的第一篇博客，希望能有个好的开始把。

代码

import arcpy
from arcpy import env
from arcpy.sa import *
import math
def Xiepojiegou():
    env.workspace = "E:\\gisforpyexcel\\mandata\\arcpyresult"
    # 这是输入的坡向和倾向（倾角和倾向都是用产状点插值生成的）
    inRaster1 = arcpy.Raster("aspect0516.tif")
    inRaster2 = arcpy.Raster("倾向裁剪.tif")
    # Check out the ArcGIS Spatial Analyst extension license
    arcpy.CheckOutExtension("Spatial")

    # 进行栅格计算
    outPlus1 = inRaster1 - inRaster2
    outPlus2= Abs(outPlus1)
    # Save the output
    outPlus2.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\pojianqinxiang.tif")
    # 切换工作空间
    env.workspace = "E:\\gisforpyexcel\\mandata\\斜坡结构"
    # 这个是重分类的设置，首先不考虑倾角大小和坡向-倾向的绝对值大小，得出顺向坡、顺斜坡、横向坡、逆斜坡、逆向坡
    a=[[0,30,0],[30,60,4],[60,120,5],[120,150,6],[150,210,7],[210,240,6],[240,300,5],[300,330,4],[330,360,0]]
    outReclass1 = Reclassify("pojianqinxiang.tif", "Value", RemapRange(a))
    outReclass1.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\reclass04567.tif")
    env.workspace = "E:\\gisforpyexcel\\mandata\\arcpyresult"
    # 对坡度和倾角大小进行判断，得出顺向飘倾斜坡、顺向层面坡、顺向伏倾坡
    inRaster3 = arcpy.Raster("slope0516.tif")
    inRaster4 = arcpy.Raster("倾角裁剪1.tif")

    outPlus2 = inRaster3 - inRaster4

    # Save the output
    outPlus2.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\podujianqinjiao.tif")
    # 对倾角以5°为间隔重分类，后面要用
    outReclassqinjiao = Reclassify("倾角裁剪1.tif", "Value",
                             RemapRange([[0, 5, 0], [5, 90, 1]]))
    outReclassqinjiao.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\Remapqinjiao05.tif")
    print 3
    env.workspace = "E:\\gisforpyexcel\\mandata\\斜坡结构"
    # 读取输出的最大值和最小值，来判别关系
    mmax = arcpy.GetRasterProperties_management("E:\\gisforpyexcel\\mandata\\斜坡结构\\podujianqinjiao.tif", "MAXIMUM")
    mmin= arcpy.GetRasterProperties_management("E:\\gisforpyexcel\\mandata\\斜坡结构\\podujianqinjiao.tif","MINIMUM")
    # 进行格式转换
    max = int(math.ceil(float(mmax[0])))
    min = int(math.floor(float(mmin[0])))
    print min
    print max
    print 4
    # 这个重分类我感觉我写的有问题，但是我不知道怎么同时运用RemapRange和RemapValue，实际应该是[min,0,1]、[0,0,2]、[0,max,3],但是会报错，因为range是[0,0]，错误的范围
    outReclass2 = Reclassify("podujianqinjiao.tif", "Value", RemapRange([[min,-0.000001,1],[-0.000001,0.000001,2],[0.000001,max,3]]))
    outReclass2.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\reclass123.tif")

    print 5
    inRaster5 = arcpy.Raster("reclass04567.tif")
    inRaster6 = arcpy.Raster("reclass123.tif")

    outPlus3 = inRaster5 + inRaster6

    # Save the output
    outPlus3.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\shunxiangpo12345678910.tif")

    # 进行重分类的目的是为了之后的栅格合并不会出现值的交叉错误，利于区分
    a = [[0, 1, 1], [1, 2, 2], [2, 3, 3],[3,10,0]]
    outReclass3 = Reclassify("shunxiangpo12345678910.tif", "Value", RemapRange(a))
    outReclass3.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\shunxiangpo.tif")

    inRaster7 = arcpy.Raster("shunxiangpo.tif")
    inRaster8 = arcpy.Raster("reclass04567.tif")

    outPlus4 = inRaster7 + inRaster8

     # Save the output，这里得到的就是不考虑倾角关系下的其余7种，但是包含近水平状坡，后面排除
    outPlus4.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\all1234567.tif")

    inRaster9 = arcpy.Raster("Remapqinjiao05.tif")
    inRaster10 = arcpy.Raster("all1234567.tif")

    outPlus4 = inRaster9 * inRaster10

    # Save the output
    outPlus4.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\zzall01234567.tif")

    a = [[-0.5, 0.5, 0], [0.5, 1.5, 1], [1.5, 2.5, 2],[2.5,3.5,3],[3.5,4.5,4],[4.5,5.5,5],[5.5,6.5,6],[6.5,7.5,7]]
    outReclass3 = Reclassify("zzall01234567.tif", "Value", RemapRange(a))
    outReclass3.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\all01234567.tif")

    arcpy.AddField_management("all01234567.tif", 'Class', "TEXT")
    # 为了避免忘记，属性表加入字段名字
    try:
        cur = arcpy.UpdateCursor("all01234567.tif")
        for row in cur:
            # print row
            if(row.Value==0):
                row.Class = "近水平层状坡"
            if(row.Value==1):
                row.Class = "顺向飘倾坡"
            if (row.Value == 2):
                row.Class = "顺向层面坡"
            if (row.Value == 3):
                row.Class = "顺向伏倾坡"
            if (row.Value == 4):
                row.Class = "顺斜坡"
            if (row.Value == 5):
                row.Class = "横向坡"
            if (row.Value == 6):
                row.Class = "逆斜坡"
            if (row.Value == 7):
                row.Class = "逆向坡"
            cur.updateRow(row)
        del cur, row
    except:
        print arcpy.GetMessages()

注意事项

在进行处理之前要确保各要素的坐标一致，否则就会出现栅格点缺失的情况。
在各位进行代码复制的时候，修改文件夹的名字就可以了。我的基础数据是放在第一个文件夹里面。