Python提取DEM数据的山脊线代码

#自动提取山脊线
import arcgisscripting
import os
import sys
#创建GP工具
gp=arcgisscripting.create(9.3)
#DEM参数设置   
inputdemlyr=sys.argv[1]
#原始DEM最大值
demmaxvalue=sys.argv[2]
#生成山脊线所在文件夹
folderpath=sys.argv[3]
gp.CheckOutExtension("Spatial")
#生成坡向图层A
OutRaster = folderpath+"\sub_demA"
try:
# Check out ArcGIS Spatial Analyst extension license
    # Process: Aspect
    print "正在生成临时坡向图层A....."
    gp.AddMessage("正在生成坡向图层A.....")
    gp.Aspect_sa(inputdemlyr, OutRaster)
    print "坡向图层A已生成"
    gp.AddMessage("坡向图层A已生成")
    #由A图层生成SOA1图层
    InMeasurementType = "DEGREE"
    soalyr1=folderpath+"\sub_demSOA1"
    gp.AddMessage("正在生成SOA1图层....")
    gp.Slope_sa(OutRaster,soalyr1,InMeasurementType)
    gp.AddMessage("SOA1图层已生成")
    #生成反地形DEM图层
    gp.AddMessage("正在根据原始DEM最大高程值生成反地形DEM....")
    outfandem=folderpath+"\sub_fandem"
    gp.Minus_sa(demmaxvalue,inputdemlyr,outfandem)
    gp.AddMessage("已生成反地形DEM")
    #求反地形坡向变率SOA2
    soalyr2=folderpath+"\sub_fandemA2"
    gp.AddMessage("正在生成反地形坡向图层A.....")
    gp.Aspect_sa(outfandem, soalyr2)
    gp.AddMessage("反地形坡向图层A已生成")
    outfandemSOA2=folderpath+"\sub_SOA2"
    gp.AddMessage("正在生成反地形坡向变率SOA2.....")
    gp.Slope_sa(soalyr2,outfandemSOA2,InMeasurementType)
    gp.AddMessage("已生成反地形坡向变率SOA2")
    #求取原始DEM中没有误差的坡向变率SOA
    #求取SOA1+SOA2
    gp.AddMessage("正在生成两图层相加....")
    soa1plussoa2=folderpath+"\soa1p2"
    gp.Plus_sa(soalyr1,outfandemSOA2,soa1plussoa2)
    gp.AddMessage("已完成栅格叠加！")
    #求取SOA1-SOA2的绝对值
    gp.AddMessage("求取SOA1-SOA2的绝对值....")
    soa1missoa2=folderpath+"\smis2"
    gp.Minus_sa(soalyr1,outfandemSOA2,soa1missoa2)
    gp.AddMessage("已完成SOA1-SOA2的绝对值")
    abssoa1minussoa2=folderpath+"\absmis2"
    gp.Abs_sa(soa1missoa2, abssoa1minussoa2)
    outminuslyr=folderpath+"\minuslyr"
    gp.Minus_sa(soa1plussoa2,abssoa1minussoa2,outminuslyr)
    lastcorrectlyr=folderpath+"\laSOA"
    #原始DEM中没有误差的坡向变率SOA
    gp.Divide_sa(outminuslyr, 2, lastcorrectlyr)
    #对原始DEM进行邻域平均值统计
    OutRaster1=folderpath+"\demstat"
    InNeighborhood = "Rectangle 11 11 Cell"
    InStatisticsType = "MEAN"
    # Check out Spatial Analyst extension license
    gp.CheckOutExtension("Spatial")
    # Process: BlockStatistics
    gp.BlockStatistics_sa(inputdemlyr,OutRaster1, InNeighborhood, InStatisticsType)
    #求取C图层
    clyer=folderpath+"\clyer"
    gp.Minus_sa(inputdemlyr,OutRaster1,clyer)
    #提取山脊线
    shanjilyer=folderpath+"\shanjixian"
    #中间图层1
    lyr1=folderpath+"\lyr1"
     #中间图层2
    lyr2=folderpath+"\lyr2"
    gp.GreaterThan_sa(clyer,0,lyr1)
    gp.GreaterThan_sa(lastcorrectlyr,70,lyr2)
    gp.AddMessage("正在生成山脊线....")
    gp.BooleanAnd_sa(lyr1,lyr2,shanjilyer)
    gp.AddMessage("正在删除全部中间图层......")
    #gp.delete_management(OutRaster)
    #gp.delete_management(soalyr1)
    #gp.delete_management(soalyr2)
    #gp.delete_management(outfandemSOA2)
    #gp.delete_management(soa1plussoa2)
    #gp.delete_management(soa1missoa2)

    #gp.delete_management(abssoa1minussoa2)
    #gp.delete_management(outminuslyr)
    #gp.delete_management(lastcorrectlyr)
    #gp.delete_management(OutRaster1)
    #gp.delete_management(clyer)
    #gp.delete_management(lyr1)
    #gp.delete_management(lyr2)
    gp.AddMessage("山脊线生成结束")
except:
    # If an error occurred while running a tool, then print the messages.
    print gp.GetMessages()