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 ，