arcpy处理DEM栅格数据

#提取山谷线
import arcpy
from arcpy import env
from arcpy import sa
workpath=''
env.workpath=workpath
env.scratchWorkspace=workpath#临时工作空间
env.overwriteOutput=True
env.extent=''#左下坐标，右上坐标
dem=arcpy.Raster(workpath+'\\DEMMerge.tif')
A=sa.Aspect(dem)#计算坡向
SOA1=sa.Slope(A,'DEGREE')#计算坡向变率
H=dem.maximum
RDEM=H-dem
B=sa.Aspect(RDEM)
SOA2=sa.Slope(B,'DEGREE')
SOA=((SOA1+SOA2)-abs(SOA1-SOA2))/2
NbrRValley=sa.NbrCircle(30)
BS=sa.BlockStatistics(dem,NbrRValley,'MEAN')
C=dem-BS
ValleyR=(C<0)&(SOA>60)
ValleyR.save(workpath+'\\ValleyR')

#根据影响因子计算最佳路径
import arcpy
from arcpy import env
from arcpy import sa
workpath=''
env.workpath=workpath
env.scratchWorkspace=workpath
env.overwriteOutput=True
env.extent=''
dem=arcpy.Raster(workpath+'\\DEMMerge.tif')
density=arcpy.Raster(workpath+'\\densityrB')
ValleyR=arcpy.Raster(workpath+'\\ValleyR')
STP='WanFoTang'
EDP='BeiChanFang'
#村庄密度重分类
def groupR(s,e,count):
    s=float(s)
    e=float(e)
    step=(e-s)/count
    lst=[]
    while True:
        lst.append(s)
        s+=step
        if s>e:
            lst[-1]=e
            break
    tlst=[]
    for i in range(len(lst)-1):
        tlst.append([lst[i],lst[i+1],i+1])
    return tlst
denmax=density.maximum
denmin=density.minimum
denre=groupR(denmin,denmax,10)
Denremape=sa.RemapRange(denre)#参数二维列表[[,,],[,,]]
DenReclassi=sa.Reclassify(density,'Value',Denremape)

#坡度数据重分类
def group(s,e,count):
    s=float(s)
    e=float(e)
    step=(e-s)/count
    lst=[]
    while True:
        lst.append(s)
        s+=step
        if s>e:
            lst[-1]=e
            break
    tlst=[]
    for i in range(len(lst)-1):
        tlst.append([lst[i],lst[i+1],len(lst)-1-i])#区间值越小，重分类越大
    tlst.reverse()
    return tlst
Slopemin=outSlope.minimum
Slopemax=outSlope.maximum
degc=group(Slopemin,Slopemax,10)
sloperemap=sa.RemapRange(degc)
slopeclassi=sa.Reclassify(outSlope,'Value',sloperemap,'NODATA')

#起伏度数据重分类
#改进group
def group(s,e,count):
    s=float(s)
    e=float(e)
    step=(e-s)/count
    print(s,e,step)
    lst=[]
    while True:
        lst.append(s)
        s+=step
        print(s)
        if '%.3f'%s=='%.3f'%e:
            lst.append(e)
            break
        elif s>e:
            lst[-1]=e
            break
    print(lst)
    tlst=[]
    for i in range(len(lst)-1):
        tlst.append([lst[i],lst[i+1],len(lst)-1-i])
    tlst.reverse()
    return tlst
NbrR=sa.NbrRectangle(10,10)
outBS=sa.BlockStatistics(dem,NbrR,'RANGE','NODATA')
sBS=outBS.minimum
eBS=outBS.maximum
BSG=group(sBS,eBS,count)
BSremap=sa.RemapRange(BSG)
BSclassi=sa.Reclassify(outBS,'Value',BSremap)

#对DEM重分类，裁剪，不具有最大值最小值属性，使用游标搜索
demclip=workpath+'\\demclip'
extent=''
demclip=arcpy.Clip_management('DEMMerge.tif',extent,demclip)
outEValue=[]
cursor=arcpy.da.SearchCursor(demclip,['VALUE'])
for row in cursor:
    outEValue.append(row.getvalue('VALUE'))
del cursor
ElevetionL=min(outEValue)
ElevetionH=max(outEValue)
EleveRemaplst=group(ElevetionL,ElevetionH,count)
DEMremap=sa.RemapRange(EleveRemaplst)
EleveReclassi=sa.Reclassify(demclip,'Value',DEMremap)


#权重分配与成本计算
cost=(slopeclassi*0.6+BSclassi*0.4)+EleveReclassi*1.5+ValleyR*5
#成本数据值越高成本越低，转化成本数据值越低成本越低
costR=cost.maximum-cost
outBKLinkRaster=workpath+'\\oblr'#成本回溯连接栅格的保存
outCostDistance=sa.CostDistance(STP,costR,'',outBKLinkRaster)
path=sa.CostPath(EDP,outCostDistance,outBKLinkRaster,'','ID')
path.save()
#村落密度栅格存在NoData值，可以通过条件函数工具处理


























在这里插入代码片