土壤穩定性評估（ArcPy實現）

一、背景

在進行區域土地開發時，需要對整個區域的土壤穩定性評估。應用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)

七、結果

再簡單豐富一下色彩：

制作土壤穩定性專題地圖
點擊視圖區左下角的【布局視圖】按鈕點或點擊【菜單欄】|【視圖】|【布局視圖】，將視圖切換為布局視圖。僅顯示 stability.圖層;將stability圖層的標簽名稱改為“土壤穩三定性等級”,選中圖層數據框后逐個添加:“指北針”、“比例尺”“圖例”及“圖名”，制作過程略。

實驗結束 byebye~~~