沟谷网络的提取及沟壑密度的计算（ArcPy实现）

一、背景

沟壑密度是描述地面被水道切割破碎程度的一个指标。沟壑密度是气候、地形、岩性、植被等因素综合影响的反映。沟壑密度越大，地面越破碎，平均坡度增大，地表物质稳定性降低，且易形成地表径流，土壤侵蚀加剧。因此，沟壑密度的测定，对于了解区域地形发育特征，水土流失监测、水土保持规划有着重要的意义。

沟壑密度也称沟谷密度或沟道密度，指单位面积内沟壑的总长度。以公里/平方公里为单位。数学表达为

Ds为沟壑密度﹔>L为研究区域内的沟壑总长度（单位:公里);A为特定研究区域的面积（单位:平方公里)。

二、目的

掌握利用水文分析工具提取沟谷网络的原理及方法。

三、要求

(1)利用水文分析工具提取研究区域的沟谷网络。

(2)计算出该研究区域的沟壑密度。

四、数据

25m分辨率的DEM数据，区域面积约有59km²(\ChP11 \Ex3)。

五、模型构建器

六、ArcPy实现

# -*- coding: utf-8 -*-

# ---------------------------------------------------------------------------

# 11-3 沟谷网络的提取及沟壑密度的计算.py

# Created on: 2021-10-12 09:38:51.00000

#   (generated by ArcGIS/ModelBuilder)

# Description:

# ---------------------------------------------------------------------------

# Import arcpy module

import arcpy

from arcpy.sa import Test

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"

Output_descent_rate_raster_data1 = "Descent_data1"

FlowDir = "FlowDir"

Sink = "Sink"

Fill_dem = "Fill_dem"

Output_descent_rate_raster_data2 = "Descent_data2"

COLUMNCOUNT = "696"

ROWCOUNT = "622"

CELLSIZEX = "25"

CELLSIZEY = "25"

FlowDir_Fill = "FlowDir_Fill"

FlowAcc = "FlowAcc"

rastercalc = "rastercalc"

Reclass_rast1 = "Reclass_rast1"

StreamT1 = "StreamT1.shp"

StreamT1_Statistics = "StreamT1_Statistics"

# 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.FlowDirection_sa(dem, FlowDir, "NORMAL", Output_descent_rate_raster_data1, "D8")

# Process: 汇

print "Process: 汇"

arcpy.gp.Sink_sa(FlowDir, Sink)

# Process: 填洼

print "Process: 填洼"

arcpy.gp.Fill_sa(dem, Fill_dem, "")

# Process: 流向 (2)

print "Process: 流向 (2)"

arcpy.gp.FlowDirection_sa(Fill_dem, FlowDir_Fill, "NORMAL", Output_descent_rate_raster_data2, "D8")

# Process: 流量

print "Process: 流量"

arcpy.gp.FlowAccumulation_sa(FlowDir_Fill, FlowAcc, "", "FLOAT", "D8")

# Process: 栅格计算器

# print "Process: 栅格计算器"

# arcpy.gp.RasterCalculator_sa("\"%FlowAcc%\" > 100", rastercalc)

# Process: 条件测试

print "Process: 条件测试"

arcpy.gp.Test_sa(FlowAcc, "value>100", rastercalc)

# Process: 重分类

print "Process: 重分类"

arcpy.gp.Reclassify_sa(rastercalc, "Value", "0 NODATA;1 1", Reclass_rast1, "DATA")

# Process: 栅格河网矢量化

print "Process: 栅格河网矢量化"

arcpy.gp.StreamToFeature_sa(Reclass_rast1, FlowDir_Fill, StreamT1, "SIMPLIFY")

# Process: 添加几何属性

print "Process: 添加几何属性"

arcpy.AddGeometryAttributes_management(StreamT1, "LENGTH", "KILOMETERS", "", "")

# Process: 汇总统计数据

print "Process: 汇总统计数据"

arcpy.Statistics_analysis(StreamT1, StreamT1_Statistics, "Length SUM", "")

# Process: 获取栅格属性

print "Process: 获取栅格属性COLUMNCOUNT"

arcpy.GetRasterProperties_management(dem, "COLUMNCOUNT", "")

# Process: 获取栅格属性 (2)

print "Process: 获取栅格属性ROWCOUNT"

arcpy.GetRasterProperties_management(dem, "ROWCOUNT", "")

# Process: 获取栅格属性 (3)

print "Process: 获取栅格属性CELLSIZEX"

arcpy.GetRasterProperties_management(dem, "CELLSIZEX", "")

# Process: 获取栅格属性 (4)

print "Process: 获取栅格属性CELLSIZEY"

arcpy.GetRasterProperties_management(dem, "CELLSIZEY", "")

# Process: 添加字段

print "Process: 添加字段"

arcpy.AddField_management(StreamT1_Statistics, "area", "DOUBLE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

# Process: 计算字段

print "Process: 计算字段，计算面积,除以1000000,是转换单位,将m²转为km²"

# arcgisscripting.ExecuteError: ERROR 999999: 执行函数时出错。无效字符

# arcpy.CalculateField_management(StreamT1_Statistics, "area", "%CELLSIZEY%*%CELLSIZEX%*%ROWCOUNT%*%COLUMNCOUNT% /1000000", "VB", "")

arcpy.CalculateField_management(StreamT1_Statistics, "area", "{}*{}*{}*{} /1000000".format(CELLSIZEY,CELLSIZEX,ROWCOUNT,COLUMNCOUNT), "VB", "")

# Process: 添加字段 (2)

print "Process: 添加字段 (2)"

arcpy.AddField_management(StreamT1_Statistics, "ds", "DOUBLE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

# Process: 计算字段 (2)

print "Process: 计算字段 单位为km/km²"

arcpy.CalculateField_management(StreamT1_Statistics, "ds", "[SUM_Length] / [area]", "VB", "")

save = ["streamt1", "streamt1_statistics"]

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)