基于ArcPy的RUSLE模型LS因子计算.docx
《基于ArcPy的RUSLE模型LS因子计算.docx》由会员分享，可在线阅读，更多相关《基于ArcPy的RUSLE模型LS因子计算.docx（17页珍藏版）》请在冰点文库上搜索。
基于ArcPy的RUSLE模型LS因子计算
基于ArcPy的RUSLE模型LS因子计算
〔FromAMLandarcgisscripting〕
半荒漠猪毛菜属
ArcPy自ArcGIS10.X后推出，是一个以成功的arcgisscripting模块为根底并继承了arcgisscripting功能进而构建而成的站点包，许多函数与类与前版本差异较大。
本文参考Hickey〔1994〕对RUSLE模型LS因子计算的AML代码的根底上，将飞天小组已移植至ArcGIS平台的python脚本，再次移植至ArcGIS10.X平台。
程序如下：
importarcpy
fromarcpyimport*
fromarcpy.saimport*
#即不需要arcgisscripting中导入，用arcpy包即可
arcpy.env.workspace="E:
\\temp"
dem_input="E:
\\temp\\DEM.tif"#输入栅格数据
wshed="E:
\\temp\\Boundary.shp"#输入流域边界数据
demunits="meters"
scf_lt5=
scf_ge5=
#定义信息提示函数
defsendmsg（msg）:
printmsg
arcpy.AddMessage（msg）
#定义一个函数，输入字符型坐标、cellsize、倍数，返回平移后的字符型坐标值，目的为保存原始小数位数不变
defStoS（s,cellsize,mult）:
stri=s.split（'.'）
inte=float（stri[0]）+mult*cellsize
returnstr（int（inte））+'.'+stri[1]
#可覆盖文件
arcpy.env.overwriteOutput=1
#判断输入DEM数据的水平和垂直方向的单位是否一致
ifdemunits==Noneordemunits.strip（）=="":
demunits="meters"
sendmsg（"使用默认单位：
meters"）
elifdemunits!
="meters"anddemunits!
="feet":
demunits="meters"
sendmsg（"DEM单位输入有误,使用默认单位meters"）
#设置结束/开始坡长累计的中断因子；为小于或大于等于5度的坡设置不同的参数
#输入坡度小于5度时，建议值为，大于等于5度时，建议值为
#scf_lt5,scf_ge5值均需小于，否那么赋予默认值
ifscf_lt5>=:
scf_lt5=
ifscf_ge5>=:
scf_ge5=
else:
ifscf_ge5>=:
scf_ge5=
sendmsg（str（scf_lt5）+","+str（scf_ge5））
#通过Describe方法获取输入DEM数据的范围和分辨率大小
dem_des=arcpy.Describe（dem_input）
cell_W=dem_des.MeanCellWidth
cell_H=dem_des.MeanCellHeight
cell_size=max（cell_W,cell_H）
cell_size=max（cell_W,cell_H）#如果格网高宽不一样，取最大值
#定义一个函数，输入字符型坐标、cellsize、倍数，返回平移后的字符型坐标值，目的为保存原始小数位数不变
extent=dem_des.extent
extent_buf=StoS（str（extent.XMin）,cell_size,-1）+""+StoS（str（extent.YMin）,cell_size,-1）+""+StoS（str（extent.XMax）,cell_size,1）+""+StoS（str（extent.YMax）
cell_size,1）
sendmsg（"做一个格网缓冲后的范围"+extent_buf）
sendmsg（"创立填充DEM——dem_fill"）
#检查Spatial工具权限，很重要的一步
arcpy.CheckOutExtension（"Spatial"）
#arcpy.Fill_sa（dem_input,"dem_fill"）
#使用Hickey对ArcGIS自带Fill功能的修改构建填充DEM;本算法使用一个格网的圆环用于单个洼地格网，用八邻域格网的最小值应用于洼地格网
arcpy.Extent="MAXOF"
arcpy.Extent=extent
arcpy.CellSize=cell_size
arcpy.CopyRaster_management（dem_input,"dem_fill2.tif"）
dem_flow=FocalFlow（"dem_fill2.tif"）
dem_flow.save（"dem_flow.tif"）
dem_fill=Con（"dem_flow"==255,FocalStatistics（"dem_fill2.tif",NbrAnnulus（1,1,"CELL"）,"MINORITY"）,"dem_fill2.tif"）
dem_fill.save（"dem_fill.tif"）#用八邻域格网的最小值应用于洼地格网
sendmsg（"根据八邻域格网值创立每个格网的流入或流出方向"）
flowdir_in=FocalFlow（"dem_fill.tif"）
flowdir_in.save（"flowdir_in.tif"）
flowdir_out=FlowDirection（"dem_fill.tif"）
flowdir_out.save（"flowdir_out.tif"）
#重新设置Environment的Extent为扩充一个格网大小的范围
arcpy.Extent="MAXOF"
arcpy.Extent=extent
sendmsg（"为dem_fill创立一个格网大小的缓冲"）
dem_fill_b=Con（IsNull（"dem_fill.tif"）,FocalStatistics（"dem_fill.tif",NbrAnnulus（1,1,"CELL"）,"MINORITY"）,"dem_fill.tif"）
dem_fill_b.save（"dem_fill_b.tif"）
#设置直角的和对角线的流向计算时的格网长度
cellorth=*cell_size
celldiag=cell_size*（2**）
#为每个格网计算坡降〔downslope〕角，修正了以前代码并重新设置平地格网〔默认度，即没有流出方向〕并<0.57（inv.tanof1%gradient）；
#建议值；新的假设是所有格网即使实际上是平的比方干湖，都有的坡度；这保证了所有格网都和流向网络有关，因而可以被赋坡度角和最终的LS因素值，
#然而它需要非常小。
sendmsg（"为每个格网计算坡降〔downslope〕角"）
deg=/math.pi
#使用shift计算不同流向的偏移量
Shift_management（"dem_fill_b.tif","dem_fill_b64.tif","0","-"+str（cell_size））
Shift_management（"dem_fill_b.tif","dem_fill_b128.tif","-"+str（cell_size）,"-"+str（cell_size））
Shift_management（"dem_fill_b.tif","dem_fill_b1.tif","-"+str（cell_size）,"0"）
Shift_management（"dem_fill_b.tif","dem_fill_b2.tif","-"+str（cell_size）,str（cell_size））
Shift_management（"dem_fill_b.tif","dem_fill_b4.tif","0",str（cell_size））
Shift_management（"dem_fill_b.tif","dem_fill_b8.tif",str（cell_size）,str（cell_size））
Shift_management（"dem_fill_b.tif","dem_fill_b16.tif",str（cell_size）,"0"）
Shift_management（"dem_fill_b.tif","dem_fill_b32.tif",str（cell_size）,"-"+str（cell_size））
#计算每个网格的坡降角
down_slp_ang2=Con（flowdir_out==64,deg*ATan（（Raster（"dem_fill_b.tif"）-Raster（"dem_fill_b64.tif"））/cellorth）,\
Con（flowdir_out==128,deg*ATan（（Raster（"dem_fill_b.tif"）-Raster（"dem_fill_b128.tif"））/cellorth）,\
Con（flowdir_out==1,deg*ATan（（Raster（"dem_fill_b.tif"）-Raster（"dem_fill_b1.tif"））/cellorth）,\
Con（flowdir_out==2,deg*ATan（（Raster（"dem_fill_b.tif"）-Raster（"dem_fill_b2.tif"））/cellorth）,\
Con（flowdir_out==4,deg*ATan（（Raster（"dem_fill_b.tif"）-Raster（"dem_fill_b4.tif"））/cellorth）,\
Con（flowdir_out==8,deg*ATan（（Raster（"dem_fill_b.tif"）-Raster（"dem_fill_b8.tif"））/cellorth）,\
Con（flowdir_out==16,deg*ATan（（Raster（"dem_fill_b.tif"）-Raster（"dem_fill_b16.tif"））/cellorth）,\
Con（flowdir_out==32,deg*ATan（（Raster（"dem_fill_b.tif"）-Raster（"dem_fill_b32.tif"））/cellorth）））））））））
down_slp_ang2.save（"down_slp_ang2.tif"）
#将等于的格网赋值为
down_slp_ang=Con（down_slp_ang2<=0,,down_slp_ang2）
down_slp_ang.save（"down_slp_ang.tif"）
#重新设置环境中Extent为原始大小，并裁减downslope格网，重命名为原始名称
arcpy.Extent="MAXOF"
arcpy.Extent=extent
sendmsg（"计算每个格网的非累计格网坡长slp_lgth_cell，考虑到直角或对角线流出方向〔暂没考虑局部高程点〕"）
slp_lgth_cell=Con（flowdir_out==2,celldiag,Con（flowdir_out==8,celldiag,Con（flowdir_out==32,celldiag,Con（flowdir_out==128,celldiag,cellorth））））
slp_lgth_cell.save（"slp_lgth_cell.tif"）
#再设置环境的Extent为缓冲范围，创立缓冲格网为0的流出方向格网
arcpy.Extent="MAXOF"
arcpy.Extent=extent
flowdir_out_b=Con（IsNull（flowdir_out）,0,flowdir_out）
flowdir_out_b.save（"flowdir_out_b.tif"）
#创立初始每个格网单元的非累计坡长〔NCSL〕，并对flowdir_in和flowdir_out做按位于运算，找到正常的流向格网，
#并设为Nodata，然后计算高点〔包括填充的洼地〕为1/2*slp_lgth_cell长度。
sendmsg（"创立初始累计坡长格网slp_lgth_cum"）
Shift_management（"flowdir_out_b.tif","flowdir_out_b64.tif","0","-"+str（cell_size））
Shift_management（"flowdir_out_b.tif","flowdir_out_b128.tif","-"+str（cell_size）,"-"+str（cell_size））
Shift_management（"flowdir_out_b.tif","flowdir_out_b1.tif","-"+str（cell_size）,"0"）
Shift_management（"flowdir_out_b.tif","flowdir_out_b2.tif","-"+str（cell_size）,str（cell_size））
Shift_management（"flowdir_out_b.tif","flowdir_out_b4.tif","0",str（cell_size））
Shift_management（"flowdir_out_b.tif","flowdir_out_b8.tif",str（cell_size）,str（cell_size））
Shift_management（"flowdir_out_b.tif","flowdir_out_b16.tif",str（cell_size）,"0"）
Shift_management（"flowdir_out_b.tif","flowdir_out_b32.tif",str（cell_size）,"-"+str（cell_size））
slp_lgth_cum=Con（BitwiseAnd（flowdir_in,64）&（Raster（"flowdir_out_b64.tif"）==4）,SetNull（BitwiseAnd（flowdir_in,64）&（Raster（"flowdir_out_b64.tif"）==4）,1）,
Con（BitwiseAnd（flowdir_in,128）&（Raster（"flowdir_out_b128.tif"）==8）,SetNull（BitwiseAnd（flowdir_in,128）&（Raster（"flowdir_out_b128.tif"）==8）,1）,
Con（BitwiseAnd（flowdir_in,1）&（Raster（"flowdir_out_b1.tif"）==16）,SetNull（BitwiseAnd（flowdir_in,1）&（Raster（"flowdir_out_b1.tif"）==16）,1）,
Con（BitwiseAnd（flowdir_in,2）&（Raster（"flowdir_out_b2.tif"）==32）,SetNull（BitwiseAnd（flowdir_in,2）&（Raster（"flowdir_out_b2.tif"）==32）,1）,
Con（BitwiseAnd（flowdir_in,4）&（Raster（"flowdir_out_b4.tif"）==64）,SetNull（BitwiseAnd（flowdir_in,4）&（Raster（"flowdir_out_b4.tif"）==64）,1）,
Con（BitwiseAnd（flowdir_in,8）&（Raster（"flowdir_out_b8.tif"）==128）,SetNull（BitwiseAnd（flowdir_in,8）&（Raster（"flowdir_out_b8.tif"）==128）,1）,
Con（BitwiseAnd（flowdir_in,16）&（Raster（"flowdir_out_b16.tif"）==1）,SetNull（BitwiseAnd（flowdir_in,16）&（Raster（"flowdir_out_b16.tif"）==1）,1）,
Con（BitwiseAnd（flowdir_in,32）&（Raster（"flowdir_out_b32.tif"）==2）,SetNull（BitwiseAnd（flowdir_in,32）&（Raster（"flowdir_out_b32.tif"）==2）,1）,*slp_lgth_cell））））））））
slp_lgth_cum.save（"slp_lgth_cum.tif"）
#设置起始坡长计算点〔高点和填充的洼地〕在所有其他格网坡长已经被决定进入每个迭代；
#起始点将有一个等于1/2它们坡长的值；起始点（局部高程点）就是周围没有其他格网单元流入，或有其他单元流入，
#但与入流单元之间坡角为零的格网单元，对应于DEM中的山顶、山脊线上的点及位于DEM边缘的点，这些点通过水流方向矩阵识别，
#识别的条件是格网单元周边各相邻点的水流方向均不知向该单元；修正了以前的代码，改变了“平地〞高点得到一个0~1/2格网坡长的值；
#新的假设是，最小累计坡长是1/2格网坡长，即使是填充洼地和“平地〞高点，从而确保每个格网LS因子的值
sendmsg（"设置起始坡长计算点slp_lgth_beg"）
slp_lgth_beg=Con（IsNull（slp_lgth_cum）,cell_size,slp_lgth_cum）
slp_lgth_beg.save（"slp_lgth_beg.tif"）
#指配坡度结束〔slope-end〕因素在累计坡长结束处；修正了以前的代码中利用RUSLE准那么建议的坡度临界弧度）来区分两个不同的侵蚀/沉积对特别小
#或特别陡的坡度；对<5%使用的参数比>=5%的大；这会使在浅滩处更容易结束侵蚀，开始沉积过程；比方，一个更高的临界值意味着需要更少的坡度降低就可以结束累计。
sendmsg（"创立结束坡长累计阈值的格网slp_lgth_fac"）
slp_end_fac=Con（down_slp_ang<,scf_lt5,scf_ge5）
slp_end_fac.save（"slp_end_fac.tif"）
#移除所有任何剩余的方向格网数据〔之前运行留下的〕
ifarcpy.Exists（"fromcell_n"）:
arcpy.Delete_management（"fromcell_n"）
ifarcpy.Exists（"fromcell_ne"）:
arcpy.Delete_management（"fromcell_ne"）
ifarcpy.Exists（"fromcell_e"）:
arcpy.Delete_management（"fromcell_e"）
ifarcpy.Exists（"fromcell_se"）:
arcpy.Delete_management（"fromcell_se"）
ifarcpy.Exists（"fromcell_s"）:
arcpy.Delete_management（"fromcell_s"）
ifarcpy.Exists（"fromcell_sw"）:
arcpy.Delete_management（"fromcell_sw"）
ifarcpy.Exists（"fromcell_w"）:
arcpy.Delete_management（"fromcell_w"）
ifarcpy.Exists（"fromcell_nw"）:
arcpy.Delete_management（"fromcell_nw"）
#修正以前版本代码中创立一系列测试nodata数据来跟踪运行过程；重新设置环境Extent为正常，用dem_fill格网作为掩膜检验缓冲格网
arcpy.Extent="MAXOF"
arcpy.Extent=extent
arcpy.Mask=dem_input
arcpy.CellSize=cell_size
ndcell=1
#修正了以前版本代码中设置迭代中nodata格网为0
arcpy.CopyRaster_management（"slp_end_fac.tif","slp_lgth_nd3.tif"）
slp_lgth_nd2=Con（Raster（"slp_lgth_nd3.tif"）>0,0）
slp_lgth_nd2.save（"slp_lgth_nd2.tif"）
warn=0
#开始为每个格网计算累计坡长的迭代循环：
依据格网单元流向，将流入当前格网单元的上游格网单元非累计坡长进行累加。
#如果当前格网单元的上游单元不知一个，那么取当前格网单元上游最大坡长值作为当前格网单元的上游累计坡长。
finished=0
n=1
slp_lgth_cum2=Raster（"slp_lgth_cum.tif"）
slp_lgth_cum2.save（"slp_lgth_cum2.tif"）
finished=0
n=1
whilenotfinished:
sendmsg（"现在开始每个格网坡长计算的第"+str（n）+"次循环！
"）
slp_lgth_prev=Raster（"slp_lgth_cum2.tif"）
slp_lgth_prev.save（"slp_lgth_prev.tif"）
count=range（1,9）
forcounterincount:
#为不同的条件设置不同的参数值
ifcounter==1:
dirfrom=4
dirpossto=64
cellcol=0
cellrow=-1
elifcounter==2:
fromcell_n=Raster（"fromcell_dir.tif"）
fromcell_n.save（"fromcell_n.tif"）
dirfrom=8
dirpossto=128
cellcol=1
cellrow=-1
elifcounter==3:
fromcell_ne=Raster（"fromcell_dir.tif"）
fromcell_ne.save（"fromcell_ne.tif"）
dirfrom=16
dirpossto=1
cellcol=1
cellrow=0
elifcounter==4:
fromcell_e=Raster（"fromcell_dir.tif"）
f
升级会员 