1、基于ArcPy的RUSLE模型LS因子计算基于ArcPy的RUSLE模型LS因子计算From AML and arcgisscripting半荒漠猪毛菜属 ArcPy 自ArcGIS 10.X后推出，是一个以成功的 arcgisscripting 模块为根底并继承了 arcgisscripting 功能进而构建而成的站点包，许多函数与类与前版本差异较大。本文参考Hickey1994对RUSLE模型LS因子计算的AML代码的根底上，将飞天小组已移植至ArcGIS平台的python脚本，再次移植至ArcGIS 10.X平台。程序如下：import arcpyfrom arcpy import *f
2、rom arcpy.sa import *#即不需要arcgisscripting中导入，用arcpy包即可arcpy.env.workspace=E:tempdem_input=E:tempDEM.tif #输入栅格数据wshed=E:tempBoundary.shp #输入流域边界数据demunits=metersscf_lt5=scf_ge5=#定义信息提示函数def sendmsg(msg): print msg arcpy.AddMessage(msg)#定义一个函数，输入字符型坐标、cellsize、倍数，返回平移后的字符型坐标值，目的为保存原始小数位数不变def StoS(s,c
3、ellsize,mult): stri=s.split(.) inte=float(stri0)+mult*cellsize return str(int(inte)+.+stri1#可覆盖文件arcpy.env.overwriteOutput = 1#判断输入DEM数据的水平和垂直方向的单位是否一致if demunits=None or demunits.strip()=: demunits=meters sendmsg(使用默认单位：meters)elif demunits!=meters and demunits!=feet: demunits=meters sendmsg( DEM单位
4、输入有误,使用默认单位meters)#设置结束/开始坡长累计的中断因子；为小于或大于等于5度的坡设置不同的参数#输入坡度小于5度时，建议值为，大于等于5度时，建议值为#scf_lt5,scf_ge5值均需小于，否那么赋予默认值if scf_lt5=: scf_lt5= if scf_ge5=: scf_ge5=else: if scf_ge5=: scf_ge5=sendmsg(str(scf_lt5)+,+str(scf_ge5)#通过Describe方法获取输入DEM数据的范围和分辨率大小dem_des=arcpy.Describe(dem_input)cell_W=dem_des.Mea
5、nCellWidthcell_H=dem_des.MeanCellHeightcell_size=max(cell_W,cell_H)cell_size=max(cell_W,cell_H) #如果格网高宽不一样，取最大值#定义一个函数，输入字符型坐标、cellsize、倍数，返回平移后的字符型坐标值，目的为保存原始小数位数不变extent=dem_des.extentextent_buf=StoS(str(extent.XMin),cell_size,-1)+ +StoS(str(extent.YMin),cell_size,-1)+ +StoS(str(extent.XMax),cell_
6、size,1)+ +StoS(str(extent.YMax),cell_size,1)sendmsg(做一个格网缓冲后的范围+extent_buf)sendmsg(创立填充DEMdem_fill)#检查Spatial工具权限，很重要的一步arcpy.CheckOutExtension(Spatial)#arcpy.Fill_sa(dem_input,dem_fill)#使用Hickey对ArcGIS自带Fill功能的修改构建填充DEM;本算法使用一个格网的圆环用于单个洼地格网，用八邻域格网的最小值应用于洼地格网arcpy.Extent=MAXOFarcpy.Extent=extentarcp
7、y.CellSize=cell_sizearcpy.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(根据
8、八邻域格网值创立每个格网的流入或流出方向)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=MAXOFarcpy.Extent=extentsendmsg(为dem_fill创立一个格网大小的缓冲) dem_fill_b = Con(IsNull(dem_fill.tif
9、),FocalStatistics(dem_fill.tif , NbrAnnulus(1, 1, CELL), MINORITY),dem_fill.tif)dem_fill_b.save(dem_fill_b.tif)#设置直角的和对角线的流向计算时的格网长度cellorth=*cell_sizecelldiag=cell_size*(2*)#为每个格网计算坡降downslope角，修正了以前代码并重新设置平地格网默认度，即没有流出方向并0.57(inv. tan of 1% gradient)；#建议值；新的假设是所有格网即使实际上是平的比方干湖，都有的坡度；这保证了所有格网都和流向网络
10、有关，因而可以被赋坡度角和最终的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_
11、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_
12、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
13、= 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
14、)-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 ,
15、 , down_slp_ang2)down_slp_ang.save(down_slp_ang.tif)#重新设置环境中Extent为原始大小，并裁减downslope格网，重命名为原始名称arcpy.Extent=MAXOFarcpy.Extent=extentsendmsg(计算每个格网的非累计格网坡长slp_lgth_cell，考虑到直角或对角线流出方向暂没考虑局部高程点)slp_lgth_cell = Con(flowdir_out = 2,celldiag,Con(flowdir_out = 8,celldiag,Con(flowdir_out = 32,celldiag,Con(f
16、lowdir_out = 128,celldiag, cellorth)slp_lgth_cell.save(slp_lgth_cell.tif)#再设置环境的Extent为缓冲范围，创立缓冲格网为0的流出方向格网arcpy.Extent=MAXOFarcpy.Extent=extentflowdir_out_b = Con(IsNull(flowdir_out) , 0 ,flowdir_out)flowdir_out_b.save(flowdir_out_b.tif)#创立初始每个格网单元的非累计坡长NCSL，并对flowdir_in和flowdir_out做按位于运算，找到正常的流向格网
17、，#并设为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(c
18、ell_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
19、.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) & ( R
20、aster(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_
21、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)
22、, 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(Bitw
23、iseAnd(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边缘的点，这些点通过水流方向矩阵识别，#识别的条件是格网单元周边各相邻点的水流方向均不知向该单元；修
24、正了以前的代码，改变了“平地高点得到一个01/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%的大；
25、这会使在浅滩处更容易结束侵蚀，开始沉积过程；比方，一个更高的临界值意味着需要更少的坡度降低就可以结束累计。sendmsg(创立结束坡长累计阈值的格网slp_lgth_fac)slp_end_fac = Con(down_slp_ang 0,0)slp_lgth_nd2.save(slp_lgth_nd2.tif)warn=0#开始为每个格网计算累计坡长的迭代循环：依据格网单元流向，将流入当前格网单元的上游格网单元非累计坡长进行累加。#如果当前格网单元的上游单元不知一个，那么取当前格网单元上游最大坡长值作为当前格网单元的上游累计坡长。finished=0n=1slp_lgth_cum2=Rast
26、er(slp_lgth_cum.tif)slp_lgth_cum2.save(slp_lgth_cum2.tif)finished=0n=1while not finished: sendmsg(现在开始每个格网坡长计算的第+str(n)+次循环！) slp_lgth_prev=Raster(slp_lgth_cum2.tif) slp_lgth_prev.save(slp_lgth_prev.tif) count=range(1,9) for counter in count: #为不同的条件设置不同的参数值 if counter=1: dirfrom=4 dirpossto=64 cell
27、col=0 cellrow=-1 elif counter=2: fromcell_n=Raster(fromcell_dir.tif) fromcell_n.save(fromcell_n.tif) dirfrom=8 dirpossto=128 cellcol=1 cellrow=-1 elif counter=3: fromcell_ne=Raster(fromcell_dir.tif) fromcell_ne.save(fromcell_ne.tif) dirfrom=16 dirpossto=1 cellcol=1 cellrow=0 elif counter=4: fromcell_e=Raster(fromcell_dir.tif) f
