土壤稳定性评估（ArcPy实现）

土壤稳定性评估（ArcPy实现）

在进行区域土地开发时，需要对整个区域的土壤稳定性评估。应用GIS空间分析方法，能够快速有效地对影响土壤稳定性的因子进行制图并评估打分，通过构建评价体系，利用叠加分析，形成土壤稳定性专题图，以为土地开发保护提供决策支持。

一、背景

在进行区域土地开发时，需要对整个区域的土壤稳定性评估。应用GIS空间分析方法，能够快速有效地对影响土壤稳定性的因子进行制图并评估打分，通过构建评价体系，利用叠加分析，形成土壤稳定性专题图，以为土地开发保护提供决策支持。

二、数据

某区域的数字高程模型和土地利用图，数字高程模型为GRID格式数据，土地利用数据为landuse. shp（\Chp13\Ex3\）。

三、要求

经专家研究，土壤稳定性评估原则如下:

(1）坡度越陡，稳定性越低。坡度分级临界值分别为:3、6°、11°、20°和30°；

(2）阴坡比阳坡稳定；

(3）土地利用类型的稳定性级别由高到低分别为:森林、水域、草原、居住用地和农耕地。

各个因子的量化分值随地理位置、重要性程度、所占比例等的不同而需分别制定，本例中使用的分值及权重见下文。最后要求做出土壤稳定性级别图，级别的制定也需人为确定，这里不做具体规定，重在了解问题解决思路和过程。

四、工作流程

基于DEM提取坡度数据，按照分级临界值进行重分类，并对每个坡度区间设定权重值;基于DEM提取坡向数据，重分类划分阴阳坡，并对两个坡向设定权重值。

将土地利用的矢量数据按照土地利用类型转换为栅格数据，再重分类设定每种土地利用类型的权重值。

综合坡度、阴阳坡和土地利用类型进行空间叠加分析加权求和，得到该区域土壤稳定性数据，最终划分等级制作土壤稳定性专题图。

工作流程如下图所示：

五、模型构建器

六、ArcPy实现

# -*- coding: utf-8 -*-# ---------------------------------------------------------------------------# 13-3 土壤稳定性评估.py# Created on: 2021-10-12 15:28:28.00000# (generated by ArcGIS/ModelBuilder)# Description: # ---------------------------------------------------------------------------# Import arcpy moduleimport arcpyimport osimport shutilimport timeprint time.asctime()path = raw_input("请输入数据所在文件夹的绝对路径：").decode("utf-8")# 开始计时time_start = time.time()paths = path + "\\result"if not os.path.exists(paths): os.mkdir(paths)else: shutil.rmtree(paths) os.mkdir(paths)# Local variables:dem = path + "\\dem"landuse = path + "\\landuse.shp" # 注意，直接读取shapefile文件，需要加后缀格式，不然后报错Slope_dem = "Slope_dem"Rec_Slope = "Rec_Slope"Aspect_dem = "Aspect_dem"Rec_Aspect = "Rec_Aspect"land_ToRaster = "land_ToRaster"Rec_landuse = "Rec_landuse"stability = "stability"# Set Geoprocessing environmentsprint "Set Geoprocessing environments"arcpy.env.scratchWorkspace = paths # 临时工作空间arcpy.env.workspace = paths # 工作空间arcpy.env.extent = dem # 处理范围arcpy.env.cellSize = dem # 像元大小arcpy.env.mask = dem # 掩膜# Process: 坡度print "Process: 坡度"arcpy.gp.Slope_sa(dem, Slope_dem, "DEGREE", "1", "PLANAR", "METER")# Process: 重分类print "Process: 重分类"arcpy.gp.Reclassify_sa(Slope_dem, "Value", "0 3 10;3 6 8;6 11 7;11 20 5;20 30 3;30 60 1", Rec_Slope, "DATA")# Process: 坡向print "Process: 坡向"arcpy.gp.Aspect_sa(dem, Aspect_dem, "PLANAR", "METER")# Process: 重分类 (2)print "Process: 重分类 (2)"arcpy.gp.Reclassify_sa(Aspect_dem, "Value", "-1 5;0 90 10;90 270 1;270 360 10", Rec_Aspect, "DATA")# Process: 面转栅格print "Process: 面转栅格"arcpy.PolygonToRaster_conversion(landuse, "LANDUSE", land_ToRaster, "CELL_CENTER", "NONE", dem)# Process: 重分类 (3)print "Process: 重分类 (3)"arcpy.gp.Reclassify_sa(land_ToRaster, "LANDUSE", "agriculture 2;river 8;grass 6;forest 10;urban 4", Rec_landuse, "DATA")# Process: 加权总和print "Process: 加权总和"arcpy.gp.WeightedSum_sa("Rec_Slope Value 0.3;Rec_Aspect Value 0.3;Rec_landuse VALUE 0.4", stability)save = ["stability"]rasters = arcpy.ListRasters()for raster in rasters: if raster.lower() not in save: print u"正在删除{}图层".format(raster) arcpy.Delete_management(raster)# 结束计时time_end = time.time()# 计算所用时间time_all = time_end - time_startprint time.asctime()print "执行完毕！>>><<< 共耗时{:.0f}分{:.2f}秒".format(time_all // 60, time_all % 60)

七、结果

再简单丰富一下色彩：

实验结束 byebye~~~

箴言：因为这些东西是非常简单的。不要抱怨自己学不会，那是因为你没有足够用心。

版权声明：本文内容由网络用户投稿，版权归原作者所有，本站不拥有其著作权，亦不承担相应法律责任。如果您发现本站中有涉嫌抄袭或描述失实的内容，请联系我们jiasou666@gmail.com 处理，核实后本网站将在24小时内删除侵权内容。