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

2023-10-12 23:31:58

一、背景

在进行区域土地开发时，需要对整个区域的土壤稳定性评估。应用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 module
import arcpy
import os
import shutil
import time

print 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 environments
print "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_start
print time.asctime()
print "执行完毕！>>><<< 共耗时{:.0f}分{:.2f}秒".format(time_all // 60, time_all % 60)