基于Arcpy自动提取清涧河流域水系特征

何灿灿,徐翠玲,宋敦江

(1.长安大学地球科学与资源学院,陕西西安710054；2.中国科学院科技政策与管理研究所,北京100190)

基于Arcpy自动提取清涧河流域水系特征

何灿灿1,徐翠玲1,宋敦江2

(1.长安大学地球科学与资源学院,陕西西安710054；2.中国科学院科技政策与管理研究所,北京100190)

水系特征的提取是进行流域水文分析的基础,快速、有效地获取流域水系信息是水文分析的内在需求。依据Arcpy可以用Python语言调用ArcGIS地理工具的特点,提出了一种基于Arcpy自动提取水系特征的方法。首先,确定水系特征提取工具和模块；其次,根据ArcHydro提取流程编写脚本；最后,运行程序成功提取出清涧河流域水系特征。实验表明,基于Arcpy自动提取水系的方法是切实可行的。提取不同阈值下的清涧河水系与 1∶25万数字化水系对比结果也显示,当阈值为4.05 km2时提取的河网与实际水系较吻合。

水系特征；自动提取；Arcpy；ArcHydroTools；清涧河流域

0 引 言

流域作为我国水资源管理的基本单元之一,其水文信息的获取十分重要。随着地理信息系统以及水文模型的发展,ArcHydro成为提取流域水系特征的有效工具。相较于传统手工提取,目前常用的ArcHydro提取方法是极大的进步,但这只是提取流域水系特征半自动化的手段。随着ArcGIS应用程序的不断更新和升级,使用Arcpy在独立的脚本环境中执行水文分析已经成为可能,也将大大节省人力和时间进行流域水系特征的提取,提高地理的自动化水平。

1 水系提取工具与模块

1.1 Arcpy

Python因其免费、开源性和跨平台等特点,成为ArcGIS地理处理中的脚本语言。Arcpy站点包将ArcGIS和Python紧密地结合在一起,为Python增加地理处理功能的函数库,即为Python提供了所有地理处理工具和用于GIS数据的多种有用函数[1]。

在ArcGIS中使用脚本处理数据可以划分为两种基本方式：一种是将脚本添加到模型或ArcToolBox中运行；另一种是在独立的脚本环境中执行,即在操作系统提示符下执行脚本。后者在独立环境下执行的方式,运用更加灵活。使用以Python编写的ArcGIS应用程序和脚本,可以访问并使用由来自多种不同领域的GIS专业人员和程序员开发的大量Python模块[2]。

1.2 ArcHydroTools

美国环境系统研究所(ESRI)联合得克萨斯大学的水资源研究中心(CRWR),根据GeoDatabase数据格式开发了ArcHydro模型[3]。基于ArcHydro模型,CRWR提供了ArcGIS环境下运行的ArcHydroTools。其主要功能是就输入的栅格DEM数据获取描述流域整体的汇流特征数据集。

2 研究区概况

清涧河是黄河右岸的一级支流,发源于陕西省安塞县,自西向东流经安塞、子长、清涧和延川等4县,至延川县苏亚河村附近汇入黄河。清涧河流域[4]为树枝状水系,左右基本对称分布,河流全长167.8 km,流域面积为4 080 km2,其中集水面积超过100 km2的支流有8条[5],其余支流集水面积不超过10 km2。该流域气候特征为大陆性暖温带季风半干旱气候,具有明显的大陆性季风气候特点,冬季寒冷漫长,夏季高温暴雨；年平均气温为9.5 ℃,多年平均降水量约500 mm；其为典型的黄土丘陵沟壑区,流域内天然植被覆盖差,土质疏松,降水较少且暴雨集中,所以水土流失现象严重。

3 水系特征提取原理

水系特征提取依据ArcHydro水文分析流程,主要包括洼地填充、水流方向确定、汇流累积量计算、河流定义与分段、流域划分和河网矢量化等6个方面。

3.1 洼地填充

洼地是指DEM中的凹陷区域,其形成原因主要有两方面：一是DEM内插或栅格精度的人为失误形成的“伪地形”；二是真实存在的洼地地形。洼地填充就是要填充数据误差所导致的洼地,而真实的洼地不需要填充。在ArcGIS水文分析模块中表现为执行Fill Sinks指令,输入参数为原始DEM数据,输出参数为填平后的栅格图层Fil。

3.2 水流方向确定

水流方向在网格中是指水流离开每个网格时的指向。计算水流方向最常用的是单流向法的D8算法[6]。假设每个网格中的水流方向只有8种可能的流向,并依照最陡坡降法,只从一个方向流出网格,分别用数字1、2、4、8、16、32、64、128来表示东、东南、南、西南、西、西北、北、东北8个方向[7]。当计算出所有单元格的方向,即可以得到流域的流向数据。在ArcHydro中,执行Flow Direction工具,输入Fil图层,即可获得流向栅格Fdr。

3.3 汇流累积量计算

在流域流向数据的基础上,计算流经每个单元格的上游所有单元格的总数可以得到汇流累积量。流域内一个网格的汇流量反映了其汇聚水流能力的强弱,汇流累积量的数值越大,代表着该区域越容易形成地表径流[7]。Flow Accumulation工具满足汇流累积量的计算,输入流向栅格Fdr,可自动生成流量栅格Fac。

3.4 河段定义与分段

当网格的汇流累积量大于等于给定的集水面积阈值时,这些网格则可定义为是河道[8]；在径流处理Stream Definition指令中,输入流量栅格Fac和给定的阈值,即可提取出河流栅格Stream。集水面积阈值是决定河网提取精度的关键参数,集水面积阈值越小,提取的河流越密集；反之,提取的河网越稀疏。

河网分级是对一个线性的河流网络进行分级的数字标识。利用地表径流模拟[9]的思想,不同级别的河网所代表的汇流累积量不同,级别越高的河网,其汇流累积量也越大。在水文研究中,级别较高的河流往往被认为是主流,级别较低的河网则是支流[10]。ArcGIS采用Strahler 分级法[11],执行Stream Segmentation命令,输入河流栅格Stream和流向栅格Fdr,即可得到河段栅格Strlnk。

3.5 流域划分

流域又称集水区域,是指流经其中的水流和其他物质从一个公共的出水口排出从而形成的一个集中的排水区域[12]。ArcGIS依据流域内每一条子流域中只有一条流经且所有栅格都指向流经的思想,先确定一个出水口点,然后结合水流方向数据,分析搜索出流经该出水口的上游所有栅格,直至所有的栅格确定相应的集水区位置。流域划分的工具指令为Catchment grid delineation,需要输入流向栅格Fdr和河段栅格Strlnk,才能生成集水区栅格Cat。

3.6 河网矢量化

Catchment polygons processing的功能是将集水区栅格Cat转化为的集水区面要素Catchment。Drainage Line Processing的功能是河段栅格Strlnk转化为的DrainageLine要素。Drainage Point Processing的功能是生成和集水区域相关的DrainageLine要素。Catchment,DrainageLine,DrainagePoint是水文分析中水文网络建立的基础。

4 水系特征自动提取工作流与脚本编写

4.1 水系特征自动提取工作流

根据Python脚本编写规则以及水系特征提取流程,可以将工作流分成6个步骤(见图1)。

图1 水文提取代码编写流程资料准备

对清涧河流域进行水系提取,首先需要清涧河流域的DEM影像资料,并将其数据命名为“qj90m.tif”,此为后续工作开展的关键；其次需要ArcGIS 10 and ArcHydro 2.0 及以上版本的应用程序；此外,文中的实验皆为在独立环境下执行脚本,所选用的编辑器为ArcGIS程序中默认安装的IDLE,将在脚本窗口编写代码,在交互式解释器中输出结果。

4.2 导入函数、工具和模块

在使用Arcpy之前,需要先导入Arcpy站点包,才能在Python中调用ArcGIS标准工具箱中所有地理处理工具。同理,在水系特征提取中所有需调用的函数、工具和模块都必须先行导入。以Arcpy站点包导入为例,代码编写如下：

import arcpy

4.2.1 工作区间设置

在Arcpy中根据系统和程序的识别,分为两种类型的路径：一系统路径,即Windows操作系统可以识别的路径,例如E: estLayers为系统路径,既能被Windows识别,又同时能被Python识别；二目录路径,仅ArcGIS可以识别的路径,如E: estLayersQJ90.GDB则是目录路径,只能在ArcGIS中识别和使用[13]。

在脚本编写时,文件路径的设置要考虑：一,如何从初始路径下调用原始影像；二,ArcHydro程序运行时所生成数据的输出路径如何创造。此外,新生成的数据有栅格数据和矢量数据,其中栅格数据存储位置应存储在系统路径,矢量数据的存储位置必须设置为目录路径。

4.2.2 水文工具调用

水系提取时,主要需调用ArcHydro工具集中[Fill Sinks],[Flow Direction],[Flow Accumulation]等多项命令,调用水文分析工具的通用语法如下：

arcpy.AddMessage(“ …”)

ArcHydroTools.(in_features,out_features)

4.2.3 执行脚本

脚本编写完成并保存以后,即可在通过点击Run>Run Module F5运行。

4.2.4 结束

如果未在脚本中添加程序结束标识,并不影响脚本的正常执行,却无法得知程序何时已运行结束。为了提示程序已经运行完成,通常情况下,会在代码编写的最后一行添加print语句例如：print ‘success’。则当脚本执行完毕后,Python Shell窗口会显示success,即为提示运行成功。

4.3 完整代码编写

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

#导入函数、工具和模块

import ArcHydroTools

import os, arcpy,sys

#工作区间设置

#调用原始影像路径

root=‘e:\test\’

dem_file=os.path.join(root, ‘qj90.tif’)

arcpy.CheckOutExtension("Spatial")

#栅格和矢量数据输出路径

DEM=‘e:\test\qj90.tif’

workspace=os.path.split(DEM)[0]+‘\’

arcpy.CreateFolder_management(workspace, "Layers")

arcpy.AddMessage(" Setting script variables...")

interws=workspace +"Layers"+os.sep

sshed_gdb_name="qj90"+".gdb"

sshed_gdb=interws+sshed_gdb_name

fdr=interws+"fdr"

fac=interws+"fac"

streams=interws+"streams"

strlnk=interws+"strlnk"

cat=interws+"cat"

catchment=sshed_gdb+os.sep+"catchment"

drainageline=sshed_gdb+os.sep+"drainageline"

drainagepoint=sshed_gdb+os.sep+"drainagepoint"

arcpy.workspace=interws

toolboxList=arcpy.ListToolboxes("archydro")

ArcHydroTools.SetTargetLocations("HydroConfig", "Layers", interws, sshed_gdb)

# 调用ArcHydro工具

fil_dem=os.path.join(workspace,‘fil.img’)

arcpy.AddMessage(" FillSinks...")

ArcHydroTools.FillSinks(DEM, fil_dem,"","")

arcpy.AddMessage(" Flow direction...")

ArcHydroTools.FlowDirection(fil_dem, fdr)

arcpy.AddMessage(" Flow accumulation...")

ArcHydroTools.FlowAccumulation(fdr, fac)

arcpy.AddMessage(" Stream definition...")

#阈值输入(此处阈值为栅格的个数)

ArcHydroTools.StreamDefinition(fac, 500, streams, "")

arcpy.AddMessage(" Stream segmentation...")

ArcHydroTools.StreamSegmentation(streams, fdr, strlnk, "", "")

arcpy.AddMessage(" Catchment grid delineation...")

ArcHydroTools.CatchmentGridDelineation(fdr, strlnk, cat)

arcpy.AddMessage(" Catchment polygons...")

ArcHydroTools.CatchmentPolyProcessing(cat, catchment)

arcpy.AddMessage(" DrainageLineProcessing...")

ArcHydroTools.DrainageLineProcessing(strlnk, fdr, drainageline)

arcpy.AddMessage(" DrainagePointProcessing...")

ArcHydroTools.DrainagePointProcessing(fac, cat, catchment, drainagepoint)

#运行结束标识

print ‘success’

5 对比分析

集水面积阈值的设定是决定河网提取精度的关键参数,它不仅决定了河流的源头并直接决定提取的河网形态[14]。本文分别取流域集水面积阈值为3.24、3.65、4.05、4.46、4.86、5.27 km2时提取河网,并得到不同阈值下清涧河的河网形态(图略),与清涧河1∶25万地形图数字化水系(图略)进行对比的结果是,集水面积阈值在4.05 km2时,所提取出的河网与实际河网比较符合。

6 结 语

使用Arcpy自动提取流域水系具明显的优势：首先,Arcpy不需要打开ArcGIS软件,即可自动提取流域水系,不但可减少手动操作步骤和时间,还能降低系统运行空间,尤其在大流域和高分辨率影像图的处理中效果更显著。其次,Arcpy不仅可以应用于流域水系提取方面,可以进行流域分析、构建水文网络模型等。第三,Python语言相较于其他程序语言,易学易读,即使是不了解Python语言的研究人员,只需将文中代码简单修改,即可用于提取其他流域的水系信息。此外,本文的不足之处在于集水面积阈值的合理选取需要多次对比分析获得；而本研究侧重于自动提取,故未对阈值选取的具体步骤进行详细的论述。

[1]邵保华, 田学志. 谈Python在Arcgis地理处理中的应用[J]. 林业勘查设计, 2012(2)： 99- 100．

[2]余咏胜, 彭艳丽, 尹言军, 等. 基于Arcpy的影像地图自动处理技术研究[J]. 测绘通报, 2015(3)： 82- 85．

[3]朱思蓉, 吴华意. Arc Hydro水文数据模型[J]. 测绘与空间地理信息, 2006, 29(5)： 87- 90．

[4]刘慧荣, 周维博, 李云排, 等. 清涧河流域近50年降水变化特征分析[J]. 水资源与水工程学报, 2013, 24(5)： 124- 127．

[5]王国庆, 庞慧, 荆新爱, 等. 清涧河流域的水文情势变化阶段及其特征[J]. 中国水土保持科学, 2005, 3(2)： 23- 27．

[6]O’CALLAGHAN J F, MARK D M. The Extraction of Drainage Networks from Digital Elevation Data[J]. Computer Vision, Graphics and Image Processing, 1984, 28(3)： 323- 344．

[7]王玉富, 王翰钊. ArcGIS环境下基于DEM的流域特征提取[J]. 湖北民族学院学报： 自然科学版, 2010, 28(4)： 439- 441．

[8]赵健, 贾忠华, 罗纨. ARCGIS环境下基于DEM的流域特征提取[J]. 水资源与水工程学报, 2006, 17(1)： 74- 76．

[9]TRIBE A S. Automated Recognition of Valley Lines and Drainage Networks from Grid Digital Elevation Models: A Review and A New Method[J]. Journal of Hydrology, 1992, 139(1- 4)： 263- 293．

[10]汤国安, 杨昕. ArcGIS地理信息系统空间分析实验教程[M]. 北京： 科学出版社, 2007: 429- 445．

[11]STRAHLER A N. Computer Vision[J]. Transaction, Amercian. Geophysical. Union, 1957, 38(6)： 913- 920.

[12]李刚, 鄂文峰, 张红红, 等. ArcGIS环境下基于DEM的信息提取及应用[J]. 吉林地质, 2010, 29(4)： 163- 166．

[13]PUAL A Z. 面向ArcGIS的Python脚本编程[M]. 李明巨, 等译. 1版. 北京： 人民邮电出版社, 2014: 112- 113．

[14]宋向阳, 吴发启, 赵龙山, 等. 基于DEM的延河流域水文特征提取与分析[J]. 干旱地区农业研究, 2012, 30(4)： 200- 206．

(责任编辑 陈 萍)

Automatic Extraction of Drainage Characteristics of Qingjian River Basin Based on Arcpy

HE Cancan1, XU Cuiling1, SONG Dunjiang2

(1. School of Earth Science and Resources, Chang’an University, Xi’an 710054, Shaanxi, China2. Institute of Policy and Management, Chinese Academy of Sciences, Beijing 100190, China)

The extraction of drainage characteristics is the basis of basin hydrologic analysis and the quick and effective basin information acquisition is the inherent requirements of hydrologic analysis. As Arcpy can use ArcGIS processing tools with Python, an Arcpy-based method that can automatically extract drainage characteristics is proposed. Firstly, the method determines drainage tools and modules. Then, the script is written according to ArcHydro. Finally, the drainage characteristics of Qingjian River are got by running the program. The experiment proves that the method is feasible. Compared the drainage system which is under six thresholds with 1∶250 000 digital drainage, it is concluded that the most suitable threshold is 4.05 km2．

drainage characteristics; automatic extraction; Arcpy; ArcHydroTools; Qingjian River Basin

2016- 05- 18

中央高校基金资助项目(310827161014)

何灿灿(1992—),女,湖北随州人,硕士研究生,研究方向为水文分析；徐翠玲(通讯作者)．

P331(241)

A

0559- 9342(2017)03- 0030- 04