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)/countlst=[]while True:lst.append(s)s+=stepif s>e:lst[-1]=ebreaktlst=[]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)/countlst=[]while True:lst.append(s)s+=stepif s>e:lst[-1]=ebreaktlst=[]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)/countprint(s,e,step)lst=[]while True:lst.append(s)s+=stepprint(s)if '%.3f'%s=='%.3f'%e:lst.append(e)breakelif s>e:lst[-1]=ebreakprint(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值，可以通过条件函数工具处理在这里插入代码片

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！