基于Python脚本语言RUSLE_LS因子计算.pdf

“RUSLE_LS.py”翻译移植过程(From AML)基于RUSLE准则和DEM数据计算LS因子。编码准备：Rick 0.Van Remortel,洛克希德马丁环境服务部，美国内华达 州拉斯维加斯，最新草稿时间为2003年。其他主要贡献人Robert J.Hickey，中央华盛顿大学，华盛顿州艾伦斯堡；Mathew E.Hamilton和Robert W.Maichle,洛克希德马丁环境服务部，美国 内华达州拉斯维加斯。RU5LE第四版修正了以前版本中S因子各成分的计算顺序，从而计算出用 于RU5LE模型更精确的L5因子。RU5LE第三版(2002年5月，2003年8月修订了最终LS栅格数据的四 舍五人问题)通过倒置坡长重新初始化代码提升了速度；调整了坡度角代码，从 而在分配最小坡度上得到更加流畅的结果；调整了单个格网长度代码，使结果更 具普遍性，解决了 ESRI Arclnfo7的“in”功能的错误，从而可移植到Arcinfo 8 上。RUSLE第二版添加了更多关于基于输入DEM数据的集水区(Watershed Catchment)配置防止误解的警告，修改了对第一版中产生格网的Nodata检查 数。原始的基于USLE的AML代码由Robert Hickey编写，U5LE第一版文档 由Hickey及其他人在1994年完成，USLE第二版由Hickey在2000年完成。1/15USLE 第二版代码由 Rick Van Remortel 和 Matt Hamilton 修改,改变了一些 关于填充洼地（fill sinks）和平地面积（flat areas）的假设，解决了处理集水 区边界上Nodata的问题，允许为不同的坡度系列分配不同的坡度临界值，以及 禾U用LS计算算法依照众多RUSLE改进文档（McCool等人1997）昨为RUSLE Handbook的第四章（Renard等人（1997）卜一篇描述“基于RUSLE的AML 程序”的期刊论文已经发表:Van Remortelz R.D.,M.E.Hamilton,and R.J.Hickey.20 0 1.Estimating the LS factor for RUSLE through iterative slope length processing of digital elevation data within Arcinfo Grid.Cartography Vol.30,No.1,Pg.27-35.测试平台：Arcinfo Workstation 8.2,WinXP移植代码运行平台：ArcG工5 9.3M9.3,WinXPAWin7使用说明：越陡越长的坡比地面上产生更高的流速，但是土壤侵蚀对5（坡度）比对L（坡长）有更强的敏感度。不规则和分割的坡型对RUSLE的影响没有体现在本 程序中。LS计算的算法时基于McCool等人在1997年的研究成果校正了坡长的水 平投影；这对GI5中通过格网或地图（x,y）而不是字段（x,y,z）计算坡长很有 用。程序中坡长的计算从高点（high points,如平缓的脊顶ridgetops）到低点（low points，比如,流域倾点 watershed pour point 或其他出口点 outlet b 行政上定义的流域可能并不合适，除非其也是水文上定义的集水区域。计算L5 因子的理想输入数据是已经裁剪过或者已经扩大到感兴趣的区域及任何相关的 2/15集水区的DEM数据。为了避免任何由于尺度导致的边界影响，最大范围应该稍 微比目标区域大一点。确保DEM的高程单位和水平距离单位是一致的（默认是 米卜L和5因子计算的输出数据需要仔细检查以确保和DEM相符以及没有重 大格式问题。如果在处理中因为浮点型格式出现问题，截短或取整或许可以解决。DEM数据中水平或垂直的点、行或边界匹配异常均能使程序不稳定或产生 不连续的坡长数据。这里的平滑算法可能会解决一些DEM不规则问题，但是同 时也会使不需要这种改正的其他DEM高程格网产生不希望的平滑或泛化。如果 使用了，DEM增强算法可能有足够的证据并应用来避免坡长的显著过分扩展。定义坡度角（theta）单位为度（坡度的反正切卜坡度临界值（）相对 坡度角变化，导致累积坡长的结束和从下一个下坡（downslope）格网开始；高 值比低值更容易导致累积坡长的结束，比如，0.7比0.5需要更小的坡度差别使 累积结束。这是一个非常值得考虑的初始设置，所以小心使用。实际地形中，当 坡度降低到一定程度时会发生由侵蚀到沉积的转变，在这里用中断因子（cutoff factor）来表现。中断因子定义为格网单元对应的坡度值与其上游流向方向的格 网单元坡度比值的临界值，这一比值反映了当前格网坡度相对于上游格网坡度的 减小程度，若当前格网与上游格网坡度比值小于中断因子（临界值），则水流到 此产生沉积，坡面结束。可见中断因子越小，对坡度减小的程度要求越高。中断 因子通常由熟悉研究区特征的专家确定，当无法得到中断因子研究成果时，也可 以使用Hickey提供的默认值，坡度大于5%时设为0.5,对于坡度小于5%的情 况，为了表现缓坡水流搬运能力较弱，坡度的微笑降低就可能导致沉积出现，将 中断因子设为0.7,即格网单元坡度小于其上游单元坡度的70%时就中止坡长累 加。3/15按照惯例，将DEM做1个格网的缓冲以避免Nodata出现在边界；这样做 足够避免很多精度裁剪后的边界错误；然而，推荐为流域DEM添加10个格网 的缓冲，从而确保任何 trapped pools 或条带状 Nodata在流域边界出现后可 以通过裁剪操作去掉。下面开始程序：import arcgisscripting,sys,os,math,string#注意：这里不耍定义gp二arcgisscripting.create(9.3),否则下面一些语句(如gp.getrasterproperties。)获取的值是。bject,而非数值型gp=arcgisscripting.create()#接收脚本工具传来的参数ws=gp.getparameterastext(O)dem_input=gp.getparameterastext(l)#输入 DEM 高程格网数据名称wshed=gp.getparameterastext(2)#输入流域边界格网数据，最后裁剪LS格网之用 demunits=gp.getparameterastext(3)#输入DEM单位，meters或feet,默认为meters scf_lt5=gp.getparameterastext(4)#scf Jt5=0.7 scf_ge5=gp.getparameterastext(5)#scf_ge5=0.5gp.workspace=ws#定义返回错误函数def sendmsg(msg):print msggp.addmessage(msg)gp.overwriteoutput=l#判断指定的目录是否存在 if not gp.Exists(ws):4/15sendmsg(指定目录不存”)else:sendmsg(”指定目录存在，可以继续。”)#判断输入DEM数据的水平和垂直方向的单位是否一致if demunits=None or demunits.strip()=：demunits=meterssendmsg(使用默认单位：meters)elif demunits!=meters and demunits!=feet：demunits=meterssendmsg(DEM单位输入有误，使用默认单位meters”)#设置结束/开始坡长累计的中断因子；为小于或大于等于5度的坡设置不同的参数#输入坡度小于5度时.，建议值为0.7,大于等于5度时，建议值为0.5#S415a1:_905值均需小于1.1,否则赋予默认值if scf_lt5=l.l:scf_lt5=0.7if scf_ge5=l.l:scf_ge5=0.5else:if scf_ge5=l.l:scf_ge5=0.5sendmsg(str(scf_lt5)+,+str(scf_ge5)#通过Describe方法获取输入DEM数据的范围和分辨率大小dem_des=gp.describe(dem_input)cell_W=dem_des.MeanCell Widthcell_H=dem_des.MeanCellHeightcell_size=max(cell_W,cell_H)#如果格网高宽不一样，取最大值#定义一个函数，输入字符型坐标、cellsize,倍数，返回平移后的字符型坐标值，目的为保留原始小数位 数不变def StoS(s,cellsize,mult):stri=s.split(.)inte=float(stri0)+mult*cellsizereturn str(int(inte)+.+strilextent=dem_des.Extent.split()extent_nor=dem_des.Extentextent_buf=StoS(extent0,cell_size,-l)+StoS(extent 1,cell_size,-l)+StoS(extent2,cell_size,l)+StoS(extent3,cell_size,l)sendmsg(原始 dem 范围为+dem_des.Extent)sendmsg(”做一个格网缓冲后的范围”+extent_buf)sendmsg(创建填充 DEM-dem_fill)#检查Spatial工具权限，很重要的一步 gp.CheckOutExtension(Spatial)5/15#gp.Fill_sa(dem_input,dem_fiir,)#使用Hickey对ArcG工S自带Fill功能的修改构建填充DEM;本算法使用一个格网的圆环用于单个洼地格网，用八邻域格网的最小值应用于洼地格网gp.Extent=maxofgp.Extent=extent_nor gp.CellSize=cell_size if gp.ExistS(dem_fill)：gp.delete_management(dem_fiii)if gp.ExistS(dem_f ill2)：gp.Delete_management(ndem_fiii2)gp.MultiOutputMapAlgebra_sa(dem_fili=+dem_input)finished=Owhile not finished=l:finished=lgp.rename_managernent(dem_fiii,dem_fiii2)gp.MultiOutputMapAlgebra_sa(dem_fili=con(focalflow(dem_fili2)=255,focalmin(dem_fill2,annulus,1,1),dem_fill2)gp.MultiOutputMapAlgebra_sa(test_grid=con(focalflow(dem_fill2)=255,0,1)gp.Delete_management(dem_fiii2)try:finished=gp.GetRasterProperties(ntest_gr id,minimum)except:sendmsg(gp.GetMessages(2)gp.Delete_management(ntest_grid)sendmsg(”根据八邻域格网值创建每个格网的流入或流出方向”)if gp.ExistS(flowdir_in)：gp.Delete_management(fiowdir_in)gp.FocalFlow_sa(dem_f lowdir_in)if gp.ExistS(flowdir_out)：gp.Delete_management(fiowdir_out)gp.FlowDirection_sa(ndem_f ill,f iowdir_out)#重新设置Environment的Extenl为扩充一个格网大小的范围 gp.Extent=MAXOFgp.Extent=extent_bufsendmsg(为创建一个格网大小的缓冲”)if gp.ExistS(dem_fill_b)：gp.Delete_management(ndem_fiii_b)gp.MultiOutputMapAlgebra_sa(dem_fiii_b=con(isnull(dem_fill),focalmin(dem_fill),dem_fill)gp.Delete_management(dem_fiii)#设置直角的和对角线的流向计算时的格网长度6/15cellorth=1.00*cell_sizecelldiag=cell_size*(2*0.5)#为每个格网计算坡降(downslope)角，修正了以前代码并重新设置平地格网(默认0.0度，即没有流出 方向)0.00并0.57(inv.tan of 1%gradient)；建议值0.1；新的假设是所有格网即使实际上是平的比如干 湖，都有0.00的坡度；这保证了所有格网都和流向网络有关，因而可以被赋坡度角和最终的L5因素值，然而它需要非常小。sendmsg(为每个格网计算坡降(downslope)角”)deg=180.0/math.piif gp.ExiStS(down_slp_ang)：gp.Delete_management(down_sip_ang)dem_fill_des=gp.describe(dem_fiii_b)#方便下面使用SHIFT函数时需要的左下角坐标 fill_extent=dem_fill_des.extent.split()#F而利用C。n()函数代替if结构进行do wn_s I p_ang的计算gp.MultiOutputMapAlgebra_sa(down_slp_ang=con(flowdir_out=64,+str(deg)+*atan(dem_fill_b-shift(+dem_fiii_b+,+fill_extentO+,+StoS(fill_extentl,cell_size,-1)+)/+str(cellorth)+),con(fiowdir_out=128,+str(deg)+*atan(dem_fill_b-shift(+dem_fiii_b+,+StoS(fill_extent0zcell_size,-l)+z+StoS(fill_ex tentl,cell_size,-l)+)/+str(celldiag)+),con(fiowdir_out=1,+str(deg)+*atan(dem_fill_b-shift(+dem_fiii_b+,+StoS(fill_extent0,cell_size,-l)+,+fill_extent1 +)/+str(cellorth)+),con(fiowdir_out=2,+str(deg)+*atan(dem_fill_b-shift(+dem_fiii_b+,+StoS(fill_extent0zcell_size,-l)+,+StoS(fill_ex tentl,cell_size,l)+)/+str(celldiag)+),con(fiowdir_out=4,+str(deg)+*atan(dem_fill_b-shift(+dem_fiii_b+,+fill_extentO+,+StoS(fill_extentl,cell_size,l)+)/+str(cellorth)+n),con(fiowdir_out=8,+str(deg)+*atan(dem_fill_b-shift(+dem_fiii_b+,+StoS(fill_extent0,cell_size,l)+,+StoS(fill_ext entl,cell_size,l)+)/+str(celldiag)+),con(fiowdir_out=16,+str(deg)+*atan(dem_fill_b-shift(+dem_fiii_b+,+StoS(fill_extent0,cell_size,l)+,+fill_extentl+)/+str(cellorth)+),con(fiowdir_out=32,+str(deg)+*atan(dem_fill_b-shift(+dem_fiii_b+,+StoS(fill_extent0,cell_size,l)+,+StoS(fill_ext entl,cell_size,-l)+)/+str(celldiag)+),o.1)#将等于0.0的格网赋值为0.1if gp.ExistS(down_slp_ang2)：gp.Delete_management(down_sip_ang2)gp.MultiOutputMapAlgebra_sa(down_slp_ang2=con(down_slp_ang=0,0.1,down_slp_ang)gp.Delete_management(down_sip_ang)7/15gp.rename_management(down_sip_ang2,down_sip_ang)#重新设置环境中Extent为原始大小，并裁减downslope格网，重命名为原始名称gp.Extent=MAXOFgp.Extent=extent_norif gp.ExistS(down_slp_ang2)：gp.Delete_management(down_sip_ang2)gp.MultiOutputMapAlgebra_sa(down_slp_ang2=down_s lp_ang)gp.Delete_management(down_sip_ang)gp.rename_management(down_sip_ang2,down_sip_ang)sendmsg(”计算每个格网的非累计格网坡长sp_igth_ceii,考虑到直角或对角线流出方向(暂没考虑局部高程点)if gp.ExistS(slp_lgth_cell)：gp.Delete_management(nsip_igth_ceii)gp.MultiOutputMapAlgebra_sa(slp_lgth_cell=con(flowdir_out=2 or flowdir_out=8 or flowdir_out=32 or flowdir_out=128,+str(celldiag)+,+str(cellorth)+)#再设置环境的Extent为缓冲范围，创建缓冲格网为。的流出方向格网gp.Extent=maxofgp.Extent=extent_bufif gp.ExistS(flowdir_out_b)：gp.Delete_management(fiowdir_out_b)gp.MultiOutputMapAlgebra_sa(f lowdir_out_b=con(isnull(flowdir_out),0,flowdir_out)gp.Delete_management(fiowdir_out)#创建初始每个格网单元的非累计坡长(NCSL),并对flowdir_in和flowdijout做按位于运算，找到正 常的流向格网，并设为Nodata,然后计算高点(包括填充的洼地)为l/2*slpgth_cell长度。sendmsg(”创建初始累计坡长格网sip_igth_cum”)if gp.ExistS(slp_lgth_cum)：gp.Delete_management(nsip_igth_cum)gp.MultiOutputMapAlgebra_sa(slp_lgth_cum=con(flowdir_in&64)and(shift(f iowdir_out_b,+fill_extentO+,+str(float(fill_extent 1)-cell_size)+)=4)or(flowdir_in&128)and(shift(fiowdir_out_bz+str(float(fill_extent0)-cell_size)+,+str(float(fill_ extentl)-cell_Size)+)=8)or (flowdir_in&1)and(shift(fiowdir out b,+str(float(fill_extent0)-cell_size)+z+fill_extentl+)=16)or(flowdir_in&2)and(shift(f iowdir_out_b,+str(float(fill_extent0)-cell_size)+,+str(float(fill_ extentl)+cell_size)+)=32)or(flowdir_in&4)and(shift(f iowdir_out_b,+fill_extentO+,+str(float(fill_extent 1)+cell_size)+)=64)or(flowdir_in&8)and(shift(fiowdir_out_b,+str(float(fill_extent0)+cell_size)+,+str(float(fill _extentl)+cell_Size)+)=128)or (flowdir_in&16)and8/15(shift(f iowdir_out_b,+str(float(fill_extent0)+cell_size)+,+fill_extentl+)=1)or(flowdir_in&32)and(shift(fiowdir_out_b,+str(float(fill_extent0)+cell_size)+,+str(float(fill _extent 1 )-cel l_size)+)=2),setnull(1=1),0.5*sip_igth_ceii)#设置起始坡长计算点(高点和填充的洼地)在所有其他格网坡长已经被决定进入每个迭代；起始点将有 一个等于1/2它们坡长的值；起始点(局部高程点)就是周围没有其他格网单元流入，或有其他单元流入，但 与入流单元之间坡角为零的格网单元，对应于DEM中的山顶、山脊线上的点及位于DEM边缘的点，这些点 通过水流方向矩阵识别，识别的条件是格网单元周边各相邻点的水流方向均不知向该单元；修正了以前的 代码，改变了“平地”高点得到一个01/2格网坡长的值；新的假设是，最小累计坡长是1/2格网坡长，即 使是填充洼地和“平地”高点，从而确保每个格网LS因子的值0.00sendmsg(”设置起始坡长计算点slp_1gth_beg)if gp.ExistS(slp_lgth_beg)：gp.Delete_management(sip_igth_beg)gp.MultiOutputMapAlgebra_sa(sip_igth_beg=con(isnull(slp_lgth_cum),4-Str(cell_Size)+,slp_lgth_cum)#指配坡度结束(slope-end)因素在累计坡长结束处；修正了以前的代码中利用RUSLE准则建议的坡度临 界5%(2.8624弧度)来区分两个不同的侵蚀/沉积对特别小或特别陡的坡度;对5%使用的参数比=5%的大;这会使在浅滩处更容易结束侵蚀，开始沉积过程；比如，一个更高的临界值意味着需要更少的坡度降低就 可以结束累计。sendmsg(”创建结束坡长累计阈值的格网sip_igth_fac”)if gp.ExistS(slp_end_fac)：gp.Delete_management(sip_end_fac)gp.MultiOutputMapAlgebra_sa(slp_end_fac=con(down_slp_ang 2.8624,+str(scf_lt5)+,+str(scf_ge5)+)#移除所有任何剩余的方向格网数据(之前运行留下的)if gp.ExistS(fromcell_n)：gp.Delete_management(fromceii_n)if gp.ExistS(fromcell_ne)：gp.Delete_management(fromceii_ne)if gp.ExistS(fromcell_e)：gp.Delete_management(fromceii_e)if gp.ExistS(fromcell_se)：gp.Delete_management(fromceii_se)if gp.ExistS(fromcell_s)：gp.Delete_management(fromceii_s)if gp.ExistS(fromcell_sw)：gp.Delete_management(fromceii_sw)if gp.ExistS(fromcell_w)：gp.Delete_management(fromceii_w)if gp.ExistS(fromcell_nw)：gp.Delete_management(fromceii_nw)#修正以前版本代码中创建一系列测试nodata数据来跟踪运行过程；重新设置环境Extent为正常，用 格网作为掩膜检验缓冲格网gp.Extent=MAXOF9/15gp.Extent=extent_norgp.Mask=dem_inputgp.CellSize=cell_size ndcell=l#修正了以前版本代码中设置迭代中nodata格网为0if gp.Existsf slp_lgth_nd2)：gp.Delete_management(sip_igth_nd2)gp.MultiOutputMapAlgebra_sa(sip_igth_nd2=0)wam=0#开始为每个格网计算累计坡长的迭代循环：依据格网单元流向，将流入当前格网单元的上游格网单元非 累计坡长进行累加。如果当前格网单元的上游单元不知一个，则取当前格网单元上游最大坡长值作为当前 格网单元的上游累计坡长。finished=On=lwhile not finished:sendmsg(现在开始每个格网坡长计算的第”+s“(n)+”次循环！”)if gp.ExistS(slp_lgth_prev)：gp.Delete_management(sip_igth_prev)gp.MultiOutputMapAlgebra_sa(sip_igth_prev=sip_igth_cum)count=range(l,9)for counter in count:#为不同的条件设置不同的参数值if counter=l:dirfrom=4dirpossto=64cellcol=0 cellrow=-lelif counter=2:gp.rename_management(fromceii_dir,fromceii_n)dirfrom=8dirpossto=128 cellcol=l cellrow=-lelif counter=3:gp.rename_management(fromceii_dir,fromceii_ne)dirfrom=16 dirpossto=l cellcol=l cellrow=0elif counter=4:gp.rename_management(frornceii_ciir,fromceii_e)dirfrom=32 dirpossto=2 cellcol=l10/15cellrow=lelif counter=5:gp.rename_management(fromceii_dir,fromceii_se)dirfrom=64dirpossto=4cellcol=0cellrow=lelif counter=6:gp.rename_management(frornceii_ciir,fromceii_s)dirfrom=128 dirpossto=8 cellcol=-l cellrow=lelif counter=7:gp.rename_management(frornceii_ciir,fromceii_sw)dirfrom=ldirpossto=16cellcol=lcellrow=0else:gp.rename_management(frornceii_ciir,fromceii_w)dirfrom=2dirpossto=32cellcol=lcellrow=-lgp.MultiOutputMapAlgebra_sa(fromcell_dir=con(A(flowdir_in&+str(dirpossto)+)or(shift(fiowdir out b,+StoS(extent0,cell_size,-l*cellcol)+z n+StoS(exte ntl,cell_size,cellrow)+)A=+str(dirfrom)+)or(down_sip_ang=(shift(down slp ang,+StoS(extent0,cell_size,-l*cellcol)+z+StoS(extent l,cell_size,cellrow)+)*sip_end_fac),shift(sip igth prev,+StoS(extent0,cell_size,-l*cellcol)+z+StoS(extent l,cell_size,cellrow)+)+shift(sip_igth_ceii,+StoS(extent0,cell_size,-l*cellcol)+,+StoS(extent l,cell_size,cellrow)+),con(isnuii(shift(sip_igth_prev,+StoS(extent0,cell_sizez-l*cellcol)+,+StoS(extentl,cell_size,cellrow)+),setnull(1=1),0)if counter=8:gp.rename_management(fromceii_dir,fromceii_nw)#在忏。0|各方向中选择最大的累计坡长if gp.ExistS(slp_lgth_cum)：gp.Delete_management(sip_igth_cum)11/15gp.MultiOutputMapAlgebra_sa(slp_lgth_cum=max(fromcell_n,fromcell_ne,fromcell_e,fromcell_se,fromcell_s,fromcell_sw,fromcell_w,fromcell_nw,slp_lgth_beg)#检查最后一次循环所有格网都有nodata=ndcellif nodata=0:fanished=l#检查其他残余nodata格网if gp.exists(slp_lgth_nd)：gp.Delete_management(sip_igth_nd)gp.MultiOutputMapAlgebra_sa(slp_lgth_nd=con(isnull(slp_lgth_cum)and(A isnull(flowdir_out_b),1,0)ndcell=0#获取slp_lgth_nd的最大值给ndcelltry:ndcel I=int(gp.getrasterproperties_management(sip_igth_nd,maximum)sendmsg(str(ndcell)except:sendmsg(gp.GetMessages(2)#修正了以前版本中监视每个循环中nodata格网是否有增加，如果两个循环后没有增加，则结束循环 开始创建LS格网。在这种情况下，有可能一个或多个很小的nodata条带出现在边界，可能在输入DEM数据 10个格网缓冲区内，而不在实际研究区域内。if gp.ExiStS(nd_chg2)：gp.Delete_management(nd_chg2)gp.MultiOutputMapAlgebra_sa(nd_chg2=con(sip_igth_nd=slp_lgth_nd2),0,1)ndchg2=0try:ndchg2=int(gp.getrasterproperties_management(nd_chg2,maximum)sendmsg(str(ndchg2)except:sendmsg(gp.GetMessages(2)nd2=ndchg2if nd2=0:finished=lwarn=l#删除上一个循环的临时方向格网。gp.Delete_management(fromceii_n)gp.Delete_management(fromceii_ne)gp.Delete_management(nfromceii_e)gp.Delete_management(fromceii_se)12/15gp.Delete_management(fromceii_s)gp.Delete_management(fromceii_sw)gp.Delete_management(fromceii_w)gp.Delete_management(fromceii_nw)if gp.Existsf slp_lgth_nd2)：gp.Delete_management(sip_igth_nd2)gp.MultiOutputMapAlgebra_sa(sip_igth_nd2=slp_lgth_nd)gp.Delete_management(sip_igth_nd)n=n+l#将最后一次循环产生的累计格网重命名为max,裁剪后再重命名回去gp.rename_management(sip_igth_cum,sip_igth_max)gp.Extent=maxofgp.Extent=extent_norif gp.Existsf slp_lgth_max2)：gp.Delete_management(nsip_igth_max2)gp.rename_management(sip_igth_max,sip_igth_rnax2)gp.MultiOutputMapAlgebra_sa(sip_igth_max=s Ip_igth_max2)gp.Delete_management(sip_igth_max2)#如果有必要的话将坡长单位从meters转换为feetif gp.Exists(slp_lgth_ft):gp.Delete_management(sip_igth_ft)if demunits=meters：gp.MultiOutputMapAlgebra_sa(sip_igth_ft=sip_igth_max/0.3048)else:gp.MultiOutputMapAlgebra_sa(sip_igth_ft=slp_igth_maxn)#修正了以前版本中根据细沟/细沟侵蚀率为RUSLE中坡长分配指数;耍明确的是草地/森林有低的敏感度;根据McCool等人准则中表格4-5(1997)0sendmsg(”计算坡度基指数m_sipexp”)if gp.Exists(m_sipexp):gp.Delete_management(m_sipexp)gp.MultiOutputMapAlgebra_sa(m_slpexp=con(down_slp_ang 0.1)and(down_slp_ang=0.2)and(down_slp_ang=0.4)and(down_slp_ang 0.85),0.08,con(down_slp_ang=0.85)and(down_slp_ang=1.4)con(down_slp_ang=2.0)con(down_slp_ang=2.6)and(down_slp_ang 2.0)and(down_slp_ang 2.6)and(down_slp_ang=3.1con(down_slp_ang=3.7and(down_slp_ang 3.7),0.28,and(down_slp_ang 5.2),0.32,con(down_slp_ang=