據(jù)處理與流域面積提取流程詳解)
簡(jiǎn)介一份黃河流域綜合地理信息數(shù)據(jù)集面向GIS、水利、環(huán)境及城鄉(xiāng)規(guī)劃領(lǐng)域的科研人員和從業(yè)者可用于地形分析、水文建模、流域管理與制圖表達(dá)。資源共178個(gè)文件以shp矢量數(shù)據(jù)、dbf屬性表、tif高程?hào)鸥?、gdb地理數(shù)據(jù)庫(kù)及xls表格為主壓縮包約36.07MB結(jié)構(gòu)清晰便于分類(lèi)使用。目前已有1068人學(xué)習(xí)下載。數(shù)據(jù)集內(nèi)含黃河流域數(shù)字高程模型、干流與支流走向、湖泊分布、城市點(diǎn)位、流域界線及黃土高原地區(qū)界等多類(lèi)地理要素依托這些數(shù)據(jù)可計(jì)算坡度坡向、提取水系網(wǎng)絡(luò)、劃定匯水范圍也可結(jié)合GIS軟件進(jìn)行污染源分布、水資源配置和生態(tài)敏感性分析。對(duì)于需要開(kāi)展黃河流域空間規(guī)劃、科研制圖或本底調(diào)查的用戶這份資料能節(jié)省大量數(shù)據(jù)預(yù)處理時(shí)間直接支撐多場(chǎng)景研究工作。1. 先聊聊我看到這套黃河流域數(shù)據(jù)的第一反應(yīng)網(wǎng)上有人分享了一套黃河流域的數(shù)據(jù)集評(píng)論區(qū)一片點(diǎn)贊。我點(diǎn)進(jìn)去翻了翻文件列表第一反應(yīng)是這種把DEM和流域資料打包整理好的活兒終于有人干了。以前自己想湊齊黃河流域的基礎(chǔ)數(shù)據(jù)得跑好幾個(gè)平臺(tái)下載幾十個(gè)分塊文件再花一整個(gè)下午拼接對(duì)齊還不一定保證坐標(biāo)系統(tǒng)一。為什么黃河流域數(shù)據(jù)歷來(lái)難湊齊這里先解釋一下。黃河從青藏高原出發(fā)穿過(guò)黃土高原、河套平原、華北平原橫跨好幾個(gè)大的地貌單元。不同地形區(qū)域的數(shù)據(jù)來(lái)源不同精度差異大發(fā)布機(jī)構(gòu)也不同。比如上游需要高精度的高山峽谷數(shù)據(jù)中游黃土高原需要能看清沖溝和塬面邊界的精細(xì)DEM下游平原又得注意建筑和灌渠的干擾。過(guò)去做一次全流域的GIS分析光前期數(shù)據(jù)整理就能把人磨到?jīng)]脾氣。這套數(shù)據(jù)集至少幫后來(lái)者省下了最繁瑣的一步說(shuō)一句“資料很詳細(xì)”真不是客套話。但拿到數(shù)據(jù)只是開(kāi)始。黃河流域的數(shù)據(jù)集里到底包含什么、怎么判斷精度夠不夠、DEM提取流域面積的流程該怎么走、黃土高原這種特殊地形有什么隱藏的坑這些問(wèn)題才是更值得展開(kāi)聊的。這篇文章我就從這套數(shù)據(jù)集出發(fā)把實(shí)際使用的經(jīng)驗(yàn)和方法一次說(shuō)清楚。2. 數(shù)據(jù)集內(nèi)容盤(pán)點(diǎn)不只有DEM2.1 數(shù)據(jù)分層能對(duì)應(yīng)哪些實(shí)際需求從標(biāo)題描述看這套數(shù)據(jù)至少包含以下圖層黃河流域城市點(diǎn)、干流線、湖泊面、流域界線、重點(diǎn)支流線、黃土高原地區(qū)界線以及DEM柵格數(shù)據(jù)。逐個(gè)分析一下它們的實(shí)際用途。城市點(diǎn)是流域內(nèi)主要城鎮(zhèn)的位置看起來(lái)簡(jiǎn)單卻是很多分析的基準(zhǔn)層。比如判斷一條支流周邊有哪些縣城、評(píng)估洪水影響范圍涉及多少居民點(diǎn)、分析污染源與城市排污口的空間關(guān)系都離不開(kāi)這個(gè)點(diǎn)層。沒(méi)有準(zhǔn)確的城市點(diǎn)空間關(guān)聯(lián)分析就缺了骨架。干流和重點(diǎn)支流線作用非常大。我拿到DEM后第一件事就是做河網(wǎng)提取但提取結(jié)果是否可靠不能自己說(shuō)了算需要參照權(quán)威的矢量水系來(lái)做對(duì)比。官方干支流線就是最直接的“標(biāo)尺”。如果提取的河網(wǎng)與矢量水系大體吻合說(shuō)明填洼、流向這些環(huán)節(jié)處理得當(dāng)如果明顯偏出那就要回頭查參數(shù)設(shè)置。所以這套數(shù)據(jù)里的干支流線既可以直接用于成圖也可以當(dāng)作水力分析的質(zhì)量驗(yàn)收參照。湖泊面層適合做蓄滯洪區(qū)和濕地分布分析。黃河流域的湖泊雖然不如南方密集但像烏梁素海、紅堿淖這些水體在區(qū)域生態(tài)研究中是繞不開(kāi)的。流域界線則決定了分析范圍匯水統(tǒng)計(jì)、面積計(jì)算都要以它為準(zhǔn)繩。讓我比較意外的是黃土高原地區(qū)界也單獨(dú)成層這種邊界通常藏在學(xué)術(shù)報(bào)告附件里能直接拿到矢量形式做水土流失分析時(shí)可以少畫(huà)不少線。2.2 DEM的底細(xì)與坐標(biāo)系問(wèn)題DEM是整個(gè)數(shù)據(jù)集的核心。按現(xiàn)在主流的公開(kāi)數(shù)據(jù)源來(lái)看黃河流域常用的DEM有三種SRTM 90米、ASTER 30米和ALOS 12.5米。標(biāo)題沒(méi)說(shuō)具體分辨率但從“資料很詳細(xì)”這個(gè)評(píng)價(jià)推斷大概率是12.5米或30米的產(chǎn)品。12.5米分辨率的數(shù)據(jù)一般由ALOS PALSAR傳感器生成可以從ASF和OpenTopography平臺(tái)下載。這個(gè)分辨率在黃土高原地區(qū)表現(xiàn)相當(dāng)出色沖溝、塬面邊界、破碎地形都能看出層次。不過(guò)要提醒一句ALOS 12.5米數(shù)據(jù)的坐標(biāo)系通常是WGS84有些國(guó)內(nèi)數(shù)據(jù)源的成果卻用的是CGCS2000或者Xian 1980。拿到手先別急著分析打開(kāi)屬性面板確認(rèn)坐標(biāo)系再統(tǒng)一投影到適合黃河流域范圍的中部Albers等積圓錐投影或UTM 49/50帶。很多人處理DEM時(shí)出現(xiàn)面積結(jié)果偏差根子往往不在算法而在坐標(biāo)系沒(méi)統(tǒng)一。2.3 12.5米分辨率是一把雙刃劍12.5米聽(tīng)上去肯定比30米更精細(xì)但實(shí)際使用是有代價(jià)的。第一是文件體量全流域范圍的12.5米DEM鑲嵌和渲染都極其吃?xún)?nèi)存筆記本電腦經(jīng)常轉(zhuǎn)圈轉(zhuǎn)到懷疑人生。第二是在河套灌區(qū)這類(lèi)人工地貌明顯的區(qū)域田壟、灌渠、排水溝在12.5米DEM上都會(huì)被識(shí)別成線性凹槽提取水系時(shí)容易混入大量人工溝渠造成偽河網(wǎng)。我的建議是差異化的用法。黃土高原溝壑區(qū)、峽谷段做重點(diǎn)分析時(shí)優(yōu)先使用12.5米數(shù)據(jù)做全流域性的概覽和制圖時(shí)30米數(shù)據(jù)更穩(wěn)計(jì)算效率也高。兩層DEM互補(bǔ)使用比一味追求高分辨率更理性。另外經(jīng)常有人搜索“DSM生成DEM”順帶說(shuō)一句DSM是包含建筑物和樹(shù)冠的地表模型要生成真正可用于水文分析的DEM需要先做濾波把非地面點(diǎn)剔除再插值成連續(xù)地面。這套數(shù)據(jù)集里如果直接提供的是DEM說(shuō)明已經(jīng)有人幫你處理好了這一層省心不少。3. DEM在流域分析中的用途從下載到出圖3.1 下載DEM的靠譜渠道與預(yù)處理心態(tài)如果手上沒(méi)有現(xiàn)成DEM需要自己下載我常用的渠道是OpenTopography和ASF。OpenTopography支持按框選范圍下載SRTM、ASTER等多種產(chǎn)品頁(yè)面操作直觀適合需要批量下載的人。ASF則是ALOS 12.5米數(shù)據(jù)的主要分發(fā)渠道注冊(cè)后可以按軌道號(hào)或地理范圍檢索下載速度也比較穩(wěn)定。下載時(shí)有幾個(gè)細(xì)節(jié)值得注意。一是要先用流域界線的外接矩形去圈定范圍不要下載整景不留余量否則數(shù)據(jù)量太大后續(xù)鑲嵌和裁剪都變慢。二是下載前確認(rèn)輸出格式和坐標(biāo)系盡量統(tǒng)一下載成GeoTIFF避免后續(xù)格式轉(zhuǎn)換引入誤差。三是記錄好每塊數(shù)據(jù)的來(lái)源和時(shí)相不同時(shí)期的DEM在河谷地區(qū)可能因水庫(kù)水位變化存在高程差異這會(huì)影響后續(xù)鑲嵌質(zhì)量。3.2 鑲嵌與裁剪先把數(shù)據(jù)拼齊從平臺(tái)下載的DEM通常是一塊一塊的拼接是繞不開(kāi)的步驟。我習(xí)慣用ArcMap里的Mosaic to New Raster做鑲嵌但有兩個(gè)經(jīng)驗(yàn)補(bǔ)充分享一下。第一鑲嵌前必須把所有分塊的像素類(lèi)型、波段數(shù)和坐標(biāo)系統(tǒng)一。比如有的塊是16位整型有的是32位浮點(diǎn)型直接拼接會(huì)在接縫處形成錯(cuò)誤的像素值過(guò)渡看起來(lái)像一個(gè)亮一個(gè)暗。第二重疊區(qū)域要檢查。如果兩塊數(shù)據(jù)的重疊部分高程值不一致說(shuō)明可能來(lái)自不同軌道或不同時(shí)相這類(lèi)區(qū)域拼接后容易出現(xiàn)網(wǎng)格狀接縫。簡(jiǎn)單的做法是鑲嵌完成后在出現(xiàn)接縫的地方用局部平滑或重采樣處理。裁剪時(shí)也有一點(diǎn)心得先用流域界線的外接矩形切一遍減少數(shù)據(jù)范圍再用精確邊界掩膜提取這樣速度更快而且邊界區(qū)域的可控性更好。3.3 用DEM計(jì)算河流流域面積的操作邏輯很多人搜索“ArcGIS用DEM計(jì)算河流流域面積”我來(lái)把這個(gè)過(guò)程拆開(kāi)講清楚。它本質(zhì)上是ArcToolbox中Hydrology工具集的流水線組合步驟并不復(fù)雜但每一步都有思維陷阱。第一步是填洼Fill。DEM里存在很多由高程噪聲造成的虛假凹陷不填掉的話水流會(huì)在里面打轉(zhuǎn)流不到下游。第二步是流向Flow Direction常用的D8算法會(huì)把每個(gè)像元的水流方向指定到最陡的下坡鄰域。第三步是流量累積Flow Accumulation統(tǒng)計(jì)每個(gè)像元上游到底匯入了多少個(gè)像元的水量這個(gè)值越大越可能是河道。第四步是設(shè)定閾值把流量累積超過(guò)某個(gè)值的柵格提出來(lái)轉(zhuǎn)換成矢量河網(wǎng)。最后一步是用Watershed工具指定出水口位置回溯生成子流域邊界。流域面積其實(shí)不是直接由一個(gè)公式算出來(lái)的而是Watershed生成的子流域矢量之后用幾何計(jì)算工具求的。關(guān)鍵在于出水口要定準(zhǔn)出水口偏移幾個(gè)像元上游集水區(qū)范圍就會(huì)明顯變形。整個(gè)流程里填洼的質(zhì)量決定了流向的正確性流向正確性又決定了后續(xù)所有結(jié)果所以預(yù)處理階段多花時(shí)間是值得的。4. 我實(shí)際做過(guò)的DEM處理流程含Python4.1 填洼、流向、累積量的一連串步驟我在處理黃河流域資料時(shí)用的是一套混合流程先用ArcMap做可視化預(yù)覽和人工檢查再用Python做批量計(jì)算遇到特殊地貌再手動(dòng)調(diào)參。整套流程按順序執(zhí)行下來(lái)基本可以覆蓋大部分分析需求。填洼不能無(wú)腦用默認(rèn)參數(shù)。默認(rèn)填洼會(huì)把所有低于周邊一個(gè)像元的坑全部填平這在較平坦的地形上沒(méi)問(wèn)題但黃土高原上有很多天然的壩地和侵蝕洼地一旦被全部填掉后續(xù)水系提取就會(huì)丟掉真實(shí)地貌信息。一個(gè)折中的辦法是設(shè)置限制條件只填充深度小于特定閾值的洼地或者采用不完全填洼的算法保留那些確實(shí)有水文意義的洼地。流向計(jì)算完成后流量累積量是一個(gè)大柵格直接看會(huì)非常花需要用累積量值的對(duì)數(shù)拉伸或分級(jí)渲染來(lái)觀察。4.2 用Python批量疊加DXF等矢量數(shù)據(jù)“DEM和DXF疊加”也是高頻搜索詞這套需求常見(jiàn)于工程項(xiàng)目手上有大量CAD格式的邊界線、房屋點(diǎn)、鉆孔位置需要疊加到地形表面上。核心問(wèn)題就是坐標(biāo)基準(zhǔn)不統(tǒng)一。我之前處理過(guò)一個(gè)項(xiàng)目甲方提供了幾百個(gè)DXF點(diǎn)文件需要與12.5米DEM疊加顯示手動(dòng)在GIS里一個(gè)個(gè)導(dǎo)入再配準(zhǔn)工作量很大。后來(lái)改成用Python的rasterio庫(kù)讀取DEM取得它的仿射變換參數(shù)和投影信息再用osgeo.ogr讀取DXF矢量先檢查投影是否一致不一致就用osr.CoordinateTransformation做投影轉(zhuǎn)換最后把矢量坐標(biāo)換算成DEM像素坐標(biāo)。這樣不僅能在圖像上疊加顯示還能直接提取出每個(gè)點(diǎn)所在位置的高程值批量輸出成表格。這套流程的核心是投影轉(zhuǎn)換必須明確不能靠肉眼對(duì)齊。4.3 ArcMap和QGIS中鑲嵌DEM的不同手感ArcMap里鑲嵌DEM遇到“顏色很怪”的問(wèn)題是常態(tài)。我調(diào)試過(guò)很多次發(fā)現(xiàn)大多數(shù)情況下不是數(shù)據(jù)壞了而是默認(rèn)渲染方式坑人。DEM灰度拉伸的范圍如果不對(duì)不同分塊的亮度反差會(huì)被放大看起來(lái)像拼接出了巨大的高程斷層。解決辦法很簡(jiǎn)單在圖層符號(hào)系統(tǒng)里選擇Percent Clip拉伸或者切換到一個(gè)統(tǒng)一的色帶馬上就會(huì)正常。QGIS處理鑲嵌更輕量一些我常用Raster菜單里的Build Virtual Raster搭建一個(gè)虛擬柵格它不復(fù)制數(shù)據(jù)只在后臺(tái)建立索引打開(kāi)大范圍DEM時(shí)幾乎不占內(nèi)存。虛擬柵格可以直接參與后續(xù)分析也可以等到最終結(jié)果需要固化時(shí)再實(shí)際合并。對(duì)于動(dòng)輒幾個(gè)GB的黃河流域DEM來(lái)說(shuō)這個(gè)特性實(shí)在太有用了。5. 黃土高原特殊地貌下的數(shù)據(jù)處理經(jīng)驗(yàn)5.1 大起伏地形中填洼參數(shù)的調(diào)整黃土高原的地貌特點(diǎn)是塬、墚、峁和沖溝交錯(cuò)高差大、細(xì)溝發(fā)育。最典型的一個(gè)坑是默認(rèn)填洼完成后提取出來(lái)的河網(wǎng)變得一段一段的像被切斷一樣。排查原因后發(fā)現(xiàn)深切的沖溝中存在大量“偽洼地”這些區(qū)域在DEM里高程比周邊低但并非真正的集水洼地而是溝道里的噪聲點(diǎn)。默認(rèn)填洼把它們?nèi)钇礁淖兞司植繀R水方向河網(wǎng)自然就亂了。解決思路有幾個(gè)方向。一是把填洼深度限制在一個(gè)較小范圍只填掉低于周?chē)欢ㄖ档南裨A魷系纼?nèi)的自然起伏。二是先用低通濾波對(duì)DEM做平滑去除高頻噪聲再填洼。我實(shí)際用下來(lái)先輕微平滑再限制深度填洼提取的河網(wǎng)連續(xù)性好很多。還有一個(gè)偏技巧的做法給DEM加上一個(gè)微小的梯度增量讓水流更傾向于朝某一個(gè)方向流動(dòng)在平緩地區(qū)尤其有效。5.2 干流與支流提取的閾值選擇提取河網(wǎng)時(shí)閾值的選擇講究一個(gè)“度”。閾值太小生成的河網(wǎng)密集得嚇人很多坡面漫流都被當(dāng)成河道閾值太大又只剩干流支流細(xì)節(jié)全丟。我不建議拍腦袋定一個(gè)常數(shù)更靠譜的方法是用參照矢量去反推。具體操作是這樣的先算出流量累積柵格的最大值記為Max然后用0.5%、1%、1.5%、2%四個(gè)閾值分別生成河網(wǎng)再和數(shù)據(jù)集里的干流及重點(diǎn)支流線做空間相交率分析。相交率最高的那個(gè)閾值就是當(dāng)前范圍內(nèi)最貼近真實(shí)河網(wǎng)的參數(shù)。這一招寫(xiě)進(jìn)自動(dòng)化腳本里可以批量處理多個(gè)子流域。在河套平原這類(lèi)人工地貌強(qiáng)烈的區(qū)域還要加一個(gè)動(dòng)作就是把灌溉渠系所在區(qū)域做成掩膜提前排除掉否則提取出的河網(wǎng)會(huì)呈現(xiàn)網(wǎng)格狀和真實(shí)自然水系完全兩回事。5.3 邊界對(duì)齊流域界線與DEM的邊緣效應(yīng)數(shù)據(jù)集里的流域界線、黃土高原界線和DEM柵格邊界往往不是嚴(yán)絲合縫的。用界線去裁剪DEM時(shí)邊界外可能會(huì)有少量DEM像元?dú)埩粼谘谀?nèi)雖然單個(gè)像元很小但周長(zhǎng)一長(zhǎng)統(tǒng)計(jì)面積和長(zhǎng)度時(shí)誤差就能累積到幾平方公里。我的辦法是在用Extract by Mask提取之后把結(jié)果和流域界線做一次柵格轉(zhuǎn)矢量后的邊界平滑剔除面積明顯小于一個(gè)像元的碎塊。如果分析項(xiàng)目對(duì)面積統(tǒng)計(jì)要求極高比如要出正式的徑流系數(shù)報(bào)告我還建議在邊界上人工檢查幾個(gè)關(guān)鍵節(jié)點(diǎn)確保這些位置沒(méi)有出現(xiàn)鋸齒偏差。這個(gè)步驟細(xì)碎但對(duì)成果質(zhì)量的影響非常直接。6. 拿到這套數(shù)據(jù)后可以怎么玩出花來(lái)6.1 干流、湖泊與城市點(diǎn)的空間關(guān)聯(lián)數(shù)據(jù)集里同時(shí)具備干流線、湖泊面和城市點(diǎn)這三層疊加起來(lái)能做的事情很多。最常規(guī)的操作是把城市點(diǎn)和干流線做一次Near分析計(jì)算每個(gè)城市到最近河道的距離這個(gè)指標(biāo)在流域水環(huán)境評(píng)價(jià)里很常用。再把湖泊面和流域界線疊加識(shí)別出位于流域邊界附近的湖泊這些湖泊對(duì)洪水?dāng)r蓄的作用和流域內(nèi)部湖泊完全不同值得單獨(dú)標(biāo)記出來(lái)做水文分析。如果想把分析做深一層還可以把DEM轉(zhuǎn)換出坡度、坡向結(jié)合城市點(diǎn)和支流位置評(píng)估城市上游匯水區(qū)的產(chǎn)流產(chǎn)沙潛力。這些分析放在以前往往要自己拼湊好幾個(gè)不同來(lái)源的圖層現(xiàn)在一份數(shù)據(jù)集基本上就能打通。6.2 黃土高原范圍線在分析中的輔助作用黃土高原地區(qū)界線單獨(dú)成層實(shí)際用途比想象中多。做統(tǒng)計(jì)報(bào)告時(shí)可以把黃土高原范圍作為掩膜單獨(dú)計(jì)算該區(qū)域內(nèi)的子流域數(shù)量、河網(wǎng)密度、高程起伏度這些指標(biāo)直接關(guān)系到水土流失治理重點(diǎn)區(qū)域的描述。也可以用來(lái)做城市點(diǎn)分布分析看看有多少城市落在黃土高原范圍內(nèi)配合土地利用數(shù)據(jù)進(jìn)一步評(píng)估裸露地表和植被覆蓋情況。我自己在做一次生態(tài)環(huán)境敏感性評(píng)估時(shí)就是拿這條界線把黃河流域切成兩個(gè)大區(qū)黃土高原區(qū)執(zhí)行一套指標(biāo)體系其他區(qū)域執(zhí)行另一套比一張圖打天下合理得多。6.3 可以繼續(xù)擴(kuò)展的分析方向如果熟用R或Python這套數(shù)據(jù)可以延伸出很多更深入的課題。比如用DEM重新劃分子流域再疊加上降雨格網(wǎng)數(shù)據(jù)做分布式水文模型的初步概化。也可以把湖泊、支流和流量累積柵格聯(lián)合起來(lái)識(shí)別潛在的洪水調(diào)蓄節(jié)點(diǎn)。還可以用DEM計(jì)算地形濕度指數(shù)與黃土高原范圍線疊加圈定優(yōu)先治理區(qū)域。數(shù)據(jù)本身是靜態(tài)的但分析思路一旦打開(kāi)那這份數(shù)據(jù)就能在不同項(xiàng)目中反復(fù)發(fā)揮價(jià)值。我自己的感受是這套黃河流域資料最討喜的地方不在于某一張DEM有多精細(xì)而是把干流、支流、湖泊、界線和DEM整合到了一個(gè)統(tǒng)一空間基準(zhǔn)里。它讓“拿到數(shù)據(jù)就能干活”這件事變成了現(xiàn)實(shí)。接下來(lái)要做的就是帶著自己的分析目標(biāo)好好用起來(lái)。本文還有配套的精品資源點(diǎn)擊獲取