基于ArcPy的RUSLE模型LS因子计算

1/12基于ArcPy的RUSLE模型LS因子计算（FromAMLandarcgisscripting）半荒漠猪毛菜属ArcPy自ArcGIS10.X后推出，是一个以成功的arcgisscripting模块为基础并继承了arcgisscripting功能进而构建而成的站点包，许多函数与类与前版本差异较大。本文参考Hickey（1994）对RUSLE模型LS因子计算的AML代码的基础上，将飞天小组已移植至ArcGIS9.3平台的python脚本，再次移植至ArcGIS10.X平台。程序如下：importarcpyfromarcpyimport*fromarcpy.saimport*#即不需要arcgisscripting中导入，用arcpy包即可arcpy.env.workspace=E:\\tempdem_input=E:\\temp\\DEM.tif#输入栅格数据wshed=E:\\temp\\Boundary.shp#输入流域边界数据demunits=metersscf_lt5=0.7scf_ge5=0.5#定义信息提示函数defsendmsg(msg):printmsgarcpy.AddMessage(msg)#定义一个函数，输入字符型坐标、cellsize、倍数，返回平移后的字符型坐标值，目的为保留原始小数位数不变defStoS(s,cellsize,mult):stri=s.split('.')inte=float(stri[0])+mult*cellsize2/12returnstr(int(inte))+'.'+stri[1]#可覆盖文件arcpy.env.overwriteOutput=1#判断输入DEM数据的水平和垂直方向的单位是否一致ifdemunits==Noneordemunits.strip()==:demunits=meterssendmsg(使用默认单位：meters)elifdemunits!=metersanddemunits!=feet:demunits=meterssendmsg(DEM单位输入有误,使用默认单位meters)#设置结束/开始坡长累计的中断因子；为小于或大于等于5度的坡设置不同的参数#输入坡度小于5度时，建议值为0.7，大于等于5度时，建议值为0.5#scf_lt5,scf_ge5值均需小于1.1，否则赋予默认值ifscf_lt5=1.1:scf_lt5=0.7ifscf_ge5=1.1:scf_ge5=0.5else:ifscf_ge5=1.1:scf_ge5=0.5sendmsg(str(scf_lt5)+,+str(scf_ge5))#通过Describe方法获取输入DEM数据的范围和分辨率大小dem_des=arcpy.Describe(dem_input)cell_W=dem_des.MeanCellWidthcell_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_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=MAXOF3/12arcpy.Extent=extentarcpy.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(根据八邻域格网值创建每个格网的流入或流出方向)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),FocalStatistics(dem_fill.tif,NbrAnnulus(1,1,CELL),MINORITY),dem_fill.tif)dem_fill_b.save(dem_fill_b.tif)#设置直角的和对角线的流向计算时的格网长度cellorth=1.00*cell_sizecelldiag=cell_size*(2**0.5)#为每个格网计算坡降（downslope）角，修正了以前代码并重新设置平地格网（默认0.0度，即没有流出方向）0.00并0.57(inv.tanof1%gradient)；#建议值0.1；新的假设是所有格网即使实际上是平的比如干湖，都有0.00的坡度；这保证了所有格网都和流向网络有关，因而可以被赋坡度角和最终的LS因素值，#然而它需要非常小。sendmsg(为每个格网计算坡降（downslope）角)deg=180.0/math.pi#使用shift计算不同流向的偏移量Shift_management(dem_fill_b.tif,dem_fill_b64.tif,0,-+str(cell4/12_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)#将等于0.0的格网赋值为0.1down_slp_ang=Con(down_slp_ang2=0,0.1,down_slp_ang2)5/12down_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(flowdir_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做按位于运算，找到正常的流向格网，#并设为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(flo