ArcGIS使用Python調(diào)用arcpy實現(xiàn)自動計算地形濕度指數(shù)(TWI)的方法

ArcGISPythonarcpy實現(xiàn)自動計算地形濕度指數(shù)（TopographicWetness的方法趙巖簡介TWITWI的計算同時使用地理信息系統(tǒng)(GIS)Python(Python)地形濕潤度指數(shù)(Schillaci。TWIetal.，1979):TWI=ln[CA/Slope]其中，CA是通過網(wǎng)格單元排水的局部上坡集水區(qū)，slope是每個網(wǎng)格單元最陡的向外坡度，以落差/距離衡量，即坡度角度的tan值(Tarboton,1997)。ArcGISPythonarcpyArcPy是一個Python站點包，可提供以實用高效的方式通過Python執(zhí)行地理數(shù)據(jù)分析、數(shù)據(jù)轉(zhuǎn)換、數(shù)據(jù)管理和地圖自動化。傳統(tǒng)方法步驟m（原文作者：0。①基于表面工具箱Surface計算Slope。②單位轉(zhuǎn)換：Scale_slope=Slope*pi/180。③填洼(Fill_dem)-Fill。④水流方向(FlowDir_Fill)-FlowDireation。⑤匯流累積量(FlowAcc_Flow)-FlowAccumulation。⑥單位面積的匯流量(SCA)-RasterCalculator。⑦地形濕度指數(shù)（TWI）-RasterCalculator：TWI=ln(SCA/Tan(Scale_slope))。Python實現(xiàn)自動化計算①購買咖啡，燒開水。②打開ArcGIS，打開Python窗口。"""

③輸入以下代碼：...TWI-filled-plus0.1constant.py...TopographicWetnessIndex...2013-11-12...JeffreyWolf(EEB);G.AndrewFricker(GEOG)...UCLA...ThisscriptwaswrittentobeusedasatoolinArcGIS.ThispythonscriptcanbeimportedtocreateaTWItool.InputsaretheworkspaceandtheinputDEM,outputistheTWIlayer...WechosetofillallsinksduetosomesmallsinksinalidarredivedDEM.Wealsoaddasmallconstanttothedenominatortoavoiddividngbyzero....TheoriginalscriptwasbasedoffthearcpyscriptwrittenbyPrasadPathak..../details.asp?dbid=16750...Thisrevisedscriptconvertstheterrainslopeindegreestoradians...Thisscriptalsousesthedefaultsettingsfortheflowaccumulationraster...howeverdifferentmethodstocalculateflowaccumulationcandramticallychangetheresultsoftheTWI..."""importarcpy,mathif name ==' main ':... arcpy.CheckOutExtension("Spatial")...... #Defineworkspaceandsetinputandoutputfiles... arcpy.env.workspace="D:/test.gdb"... inDEM"E:/DATA/Goulinping/GLPDEM剪切.tif"...... #Intermediates... print("FillingDEM.")... DEM_filled=arcpy.sa.Fill(inDEM)...... print("Creatingflowdirection.")... outFlowDirection=arcpy.sa.FlowDirection(DEM_filled,"FORCE")...... print("Creatingflowaccumulation.")... #outFlowAccumulation=arcpy.sa.FlowAccumulation(outFlowDirection,"","FLOAT")+... outFlowAccumulation=arcpy.sa.FlowAccumulation(outFlowDirection,"","INTEGER")+...... print("Creatingslope.")... slope=arcpy.sa.Slope(DEM_filled)...... print("Convertingslopeindegreestoslopeinradians")... #2Piradians=360degrees... #Piradians=180degrees... #conversion:Piradians/180degress... slope_radians=slope*math.pi/180.0...... #Output... print("CreatingTWI")... TWI=arcpy.sa.Ln(outFlowAccumulation(arcpy.sa.Tan(slope_radians)+.01))... TWI.save("E:/DATA/Goulinping/GLP_TWI.tif")... print("SavedTWI.Done.")注意修改輸入的DEM和輸出的文件路徑和名稱。④點擊回車運行代碼。⑤喝咖啡，O(∩_∩)O哈哈~。⑥出現(xiàn)提示：FillingDEM.Creatingflowdirection.Creatingflowaccumulation.Creatingslope.ConvertingslopeindegreestoslopeinradiansCreatingTWISavedTWI.Done.表示計算完成。References:Beven,K.J.,M.J.Kirkby,andJ.Seibert.1979.Aphysicallybased,variablecontributingareamodelofbasinhydrology.HydrologicalScienceBulletin24:43-69.Schillaci,C.,A.Braun,andJ.Kropacek.2015.Terrainanalysisandlandformrecognition;Chapter2.4.2,inGeomorphologicalTechniques;BritishSocietyforGeomorphology.18pp.Tarbo