何燦燦,徐翠玲,宋敦江

(1.長安大學地球科學與資源學院,陜西西安710054；2.中國科學院科技政策與管理研究所,北京100190)

基于Arcpy自動提取清澗河流域水系特征

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

(1.長安大學地球科學與資源學院,陜西西安710054；2.中國科學院科技政策與管理研究所,北京100190)

水系特征的提取是進行流域水文分析的基礎,快速、有效地獲取流域水系信息是水文分析的內(nèi)在需求。依據(jù)Arcpy可以用Python語言調(diào)用ArcGIS地理工具的特點,提出了一種基于Arcpy自動提取水系特征的方法。首先,確定水系特征提取工具和模塊；其次,根據(jù)ArcHydro提取流程編寫腳本；最后,運行程序成功提取出清澗河流域水系特征。實驗表明,基于Arcpy自動提取水系的方法是切實可行的。提取不同閾值下的清澗河水系與 1∶25萬數(shù)字化水系對比結(jié)果也顯示,當閾值為4.05 km2時提取的河網(wǎng)與實際水系較吻合。

水系特征；自動提??；Arcpy；ArcHydroTools；清澗河流域

0 引 言

流域作為我國水資源管理的基本單元之一,其水文信息的獲取十分重要。隨著地理信息系統(tǒng)以及水文模型的發(fā)展,ArcHydro成為提取流域水系特征的有效工具。相較于傳統(tǒng)手工提取,目前常用的ArcHydro提取方法是極大的進步,但這只是提取流域水系特征半自動化的手段。隨著ArcGIS應用程序的不斷更新和升級,使用Arcpy在獨立的腳本環(huán)境中執(zhí)行水文分析已經(jīng)成為可能,也將大大節(jié)省人力和時間進行流域水系特征的提取,提高地理的自動化水平。

1 水系提取工具與模塊

1.1 Arcpy

Python因其免費、開源性和跨平臺等特點,成為ArcGIS地理處理中的腳本語言。Arcpy站點包將ArcGIS和Python緊密地結(jié)合在一起,為Python增加地理處理功能的函數(shù)庫,即為Python提供了所有地理處理工具和用于GIS數(shù)據(jù)的多種有用函數(shù)[1]。

在ArcGIS中使用腳本處理數(shù)據(jù)可以劃分為兩種基本方式：一種是將腳本添加到模型或ArcToolBox中運行；另一種是在獨立的腳本環(huán)境中執(zhí)行,即在操作系統(tǒng)提示符下執(zhí)行腳本。后者在獨立環(huán)境下執(zhí)行的方式,運用更加靈活。使用以Python編寫的ArcGIS應用程序和腳本,可以訪問并使用由來自多種不同領域的GIS專業(yè)人員和程序員開發(fā)的大量Python模塊[2]。

1.2 ArcHydroTools

美國環(huán)境系統(tǒng)研究所(ESRI)聯(lián)合得克薩斯大學的水資源研究中心(CRWR),根據(jù)GeoDatabase數(shù)據(jù)格式開發(fā)了ArcHydro模型[3]?；贏rcHydro模型,CRWR提供了ArcGIS環(huán)境下運行的ArcHydroTools。其主要功能是就輸入的柵格DEM數(shù)據(jù)獲取描述流域整體的匯流特征數(shù)據(jù)集。

2 研究區(qū)概況

清澗河是黃河右岸的一級支流,發(fā)源于陜西省安塞縣,自西向東流經(jīng)安塞、子長、清澗和延川等4縣,至延川縣蘇亞河村附近匯入黃河。清澗河流域[4]為樹枝狀水系,左右基本對稱分布,河流全長167.8 km,流域面積為4 080 km2,其中集水面積超過100 km2的支流有8條[5],其余支流集水面積不超過10 km2。該流域氣候特征為大陸性暖溫帶季風半干旱氣候,具有明顯的大陸性季風氣候特點,冬季寒冷漫長,夏季高溫暴雨；年平均氣溫為9.5 ℃,多年平均降水量約500 mm；其為典型的黃土丘陵溝壑區(qū),流域內(nèi)天然植被覆蓋差,土質(zhì)疏松,降水較少且暴雨集中,所以水土流失現(xiàn)象嚴重。

3 水系特征提取原理

水系特征提取依據(jù)ArcHydro水文分析流程,主要包括洼地填充、水流方向確定、匯流累積量計算、河流定義與分段、流域劃分和河網(wǎng)矢量化等6個方面。

3.1 洼地填充

洼地是指DEM中的凹陷區(qū)域,其形成原因主要有兩方面：一是DEM內(nèi)插或柵格精度的人為失誤形成的“偽地形”；二是真實存在的洼地地形。洼地填充就是要填充數(shù)據(jù)誤差所導致的洼地,而真實的洼地不需要填充。在ArcGIS水文分析模塊中表現(xiàn)為執(zhí)行Fill Sinks指令,輸入?yún)?shù)為原始DEM數(shù)據(jù),輸出參數(shù)為填平后的柵格圖層Fil。

3.2 水流方向確定

水流方向在網(wǎng)格中是指水流離開每個網(wǎng)格時的指向。計算水流方向最常用的是單流向法的D8算法[6]。假設每個網(wǎng)格中的水流方向只有8種可能的流向,并依照最陡坡降法,只從一個方向流出網(wǎng)格,分別用數(shù)字1、2、4、8、16、32、64、128來表示東、東南、南、西南、西、西北、北、東北8個方向[7]。當計算出所有單元格的方向,即可以得到流域的流向數(shù)據(jù)。在ArcHydro中,執(zhí)行Flow Direction工具,輸入Fil圖層,即可獲得流向柵格Fdr。

3.3 匯流累積量計算

在流域流向數(shù)據(jù)的基礎上,計算流經(jīng)每個單元格的上游所有單元格的總數(shù)可以得到匯流累積量。流域內(nèi)一個網(wǎng)格的匯流量反映了其匯聚水流能力的強弱,匯流累積量的數(shù)值越大,代表著該區(qū)域越容易形成地表徑流[7]。Flow Accumulation工具滿足匯流累積量的計算,輸入流向柵格Fdr,可自動生成流量柵格Fac。

3.4 河段定義與分段

當網(wǎng)格的匯流累積量大于等于給定的集水面積閾值時,這些網(wǎng)格則可定義為是河道[8]；在徑流處理Stream Definition指令中,輸入流量柵格Fac和給定的閾值,即可提取出河流柵格Stream。集水面積閾值是決定河網(wǎng)提取精度的關鍵參數(shù),集水面積閾值越小,提取的河流越密集；反之,提取的河網(wǎng)越稀疏。

河網(wǎng)分級是對一個線性的河流網(wǎng)絡進行分級的數(shù)字標識。利用地表徑流模擬[9]的思想,不同級別的河網(wǎng)所代表的匯流累積量不同,級別越高的河網(wǎng),其匯流累積量也越大。在水文研究中,級別較高的河流往往被認為是主流,級別較低的河網(wǎng)則是支流[10]。ArcGIS采用Strahler 分級法[11],執(zhí)行Stream Segmentation命令,輸入河流柵格Stream和流向柵格Fdr,即可得到河段柵格Strlnk。

3.5 流域劃分

流域又稱集水區(qū)域,是指流經(jīng)其中的水流和其他物質(zhì)從一個公共的出水口排出從而形成的一個集中的排水區(qū)域[12]。ArcGIS依據(jù)流域內(nèi)每一條子流域中只有一條流經(jīng)且所有柵格都指向流經(jīng)的思想,先確定一個出水口點,然后結(jié)合水流方向數(shù)據(jù),分析搜索出流經(jīng)該出水口的上游所有柵格,直至所有的柵格確定相應的集水區(qū)位置。流域劃分的工具指令為Catchment grid delineation,需要輸入流向柵格Fdr和河段柵格Strlnk,才能生成集水區(qū)柵格Cat。

3.6 河網(wǎng)矢量化

Catchment polygons processing的功能是將集水區(qū)柵格Cat轉(zhuǎn)化為的集水區(qū)面要素Catchment。Drainage Line Processing的功能是河段柵格Strlnk轉(zhuǎn)化為的DrainageLine要素。Drainage Point Processing的功能是生成和集水區(qū)域相關的DrainageLine要素。Catchment,DrainageLine,DrainagePoint是水文分析中水文網(wǎng)絡建立的基礎。

4 水系特征自動提取工作流與腳本編寫

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

根據(jù)Python腳本編寫規(guī)則以及水系特征提取流程,可以將工作流分成6個步驟(見圖1)。

圖1 水文提取代碼編寫流程資料準備

對清澗河流域進行水系提取,首先需要清澗河流域的DEM影像資料,并將其數(shù)據(jù)命名為“qj90m.tif”,此為后續(xù)工作開展的關鍵；其次需要ArcGIS 10 and ArcHydro 2.0 及以上版本的應用程序；此外,文中的實驗皆為在獨立環(huán)境下執(zhí)行腳本,所選用的編輯器為ArcGIS程序中默認安裝的IDLE,將在腳本窗口編寫代碼,在交互式解釋器中輸出結(jié)果。

4.2 導入函數(shù)、工具和模塊

在使用Arcpy之前,需要先導入Arcpy站點包,才能在Python中調(diào)用ArcGIS標準工具箱中所有地理處理工具。同理,在水系特征提取中所有需調(diào)用的函數(shù)、工具和模塊都必須先行導入。以Arcpy站點包導入為例,代碼編寫如下：

import arcpy

4.2.1 工作區(qū)間設置

在Arcpy中根據(jù)系統(tǒng)和程序的識別,分為兩種類型的路徑：一系統(tǒng)路徑,即Windows操作系統(tǒng)可以識別的路徑,例如E: estLayers為系統(tǒng)路徑,既能被Windows識別,又同時能被Python識別；二目錄路徑,僅ArcGIS可以識別的路徑,如E: estLayersQJ90.GDB則是目錄路徑,只能在ArcGIS中識別和使用[13]。

在腳本編寫時,文件路徑的設置要考慮：一,如何從初始路徑下調(diào)用原始影像；二,ArcHydro程序運行時所生成數(shù)據(jù)的輸出路徑如何創(chuàng)造。此外,新生成的數(shù)據(jù)有柵格數(shù)據(jù)和矢量數(shù)據(jù),其中柵格數(shù)據(jù)存儲位置應存儲在系統(tǒng)路徑,矢量數(shù)據(jù)的存儲位置必須設置為目錄路徑。

4.2.2 水文工具調(diào)用

水系提取時,主要需調(diào)用ArcHydro工具集中[Fill Sinks],[Flow Direction],[Flow Accumulation]等多項命令,調(diào)用水文分析工具的通用語法如下：

arcpy.AddMessage(“ …”)

ArcHydroTools.(in_features,out_features)

4.2.3 執(zhí)行腳本

腳本編寫完成并保存以后,即可在通過點擊Run>Run Module F5運行。

4.2.4 結(jié)束

如果未在腳本中添加程序結(jié)束標識,并不影響腳本的正常執(zhí)行,卻無法得知程序何時已運行結(jié)束。為了提示程序已經(jīng)運行完成,通常情況下,會在代碼編寫的最后一行添加print語句例如：print ‘success’。則當腳本執(zhí)行完畢后,Python Shell窗口會顯示success,即為提示運行成功。

4.3 完整代碼編寫

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

#導入函數(shù)、工具和模塊

import ArcHydroTools

import os, arcpy,sys

#工作區(qū)間設置

#調(diào)用原始影像路徑

root=‘e:\test\’

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

arcpy.CheckOutExtension("Spatial")

#柵格和矢量數(shù)據(jù)輸出路徑

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)

# 調(diào)用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...")

#閾值輸入(此處閾值為柵格的個數(shù))

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)

#運行結(jié)束標識

print ‘success’

5 對比分析

集水面積閾值的設定是決定河網(wǎng)提取精度的關鍵參數(shù),它不僅決定了河流的源頭并直接決定提取的河網(wǎng)形態(tài)[14]。本文分別取流域集水面積閾值為3.24、3.65、4.05、4.46、4.86、5.27 km2時提取河網(wǎng),并得到不同閾值下清澗河的河網(wǎng)形態(tài)(圖略),與清澗河1∶25萬地形圖數(shù)字化水系(圖略)進行對比的結(jié)果是,集水面積閾值在4.05 km2時,所提取出的河網(wǎng)與實際河網(wǎng)比較符合。

6 結(jié) 語

使用Arcpy自動提取流域水系具明顯的優(yōu)勢：首先,Arcpy不需要打開ArcGIS軟件,即可自動提取流域水系,不但可減少手動操作步驟和時間,還能降低系統(tǒng)運行空間,尤其在大流域和高分辨率影像圖的處理中效果更顯著。其次,Arcpy不僅可以應用于流域水系提取方面,可以進行流域分析、構(gòu)建水文網(wǎng)絡模型等。第三,Python語言相較于其他程序語言,易學易讀,即使是不了解Python語言的研究人員,只需將文中代碼簡單修改,即可用于提取其他流域的水系信息。此外,本文的不足之處在于集水面積閾值的合理選取需要多次對比分析獲得；而本研究側(cè)重于自動提取,故未對閾值選取的具體步驟進行詳細的論述。

[1]邵保華, 田學志. 談Python在Arcgis地理處理中的應用[J]. 林業(yè)勘查設計, 2012(2)： 99- 100．

[2]余詠勝, 彭艷麗, 尹言軍, 等. 基于Arcpy的影像地圖自動處理技術研究[J]. 測繪通報, 2015(3)： 82- 85．

[3]朱思蓉, 吳華意. Arc Hydro水文數(shù)據(jù)模型[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環(huán)境下基于DEM的流域特征提取[J]. 湖北民族學院學報： 自然科學版, 2010, 28(4)： 439- 441．

[8]趙健, 賈忠華, 羅紈. ARCGIS環(huán)境下基于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地理信息系統(tǒng)空間分析實驗教程[M]. 北京： 科學出版社, 2007: 429- 445．

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

[12]李剛, 鄂文峰, 張紅紅, 等. ArcGIS環(huán)境下基于DEM的信息提取及應用[J]. 吉林地質(zhì), 2010, 29(4)： 163- 166．

[13]PUAL A Z. 面向ArcGIS的Python腳本編程[M]. 李明巨, 等譯. 1版. 北京： 人民郵電出版社, 2014: 112- 113．

[14]宋向陽, 吳發(fā)啟, 趙龍山, 等. 基于DEM的延河流域水文特征提取與分析[J]. 干旱地區(qū)農(nóng)業(yè)研究, 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

中央高?；鹳Y助項目(310827161014)

何燦燦(1992—),女,湖北隨州人,碩士研究生,研究方向為水文分析；徐翠玲(通訊作者)．

P331(241)

A

0559- 9342(2017)03- 0030- 04