據(jù)處理全流程:從7z解壓到坡度與流域提取)
簡(jiǎn)介這份資源面向GIS從業(yè)者、地理信息專(zhuān)業(yè)師生及區(qū)域規(guī)劃研究人員提供貴州省黔南布依族苗族自治州12.5米分辨率數(shù)字高程模型并附帶市級(jí)行政范圍矢量邊界可用于地形分析、環(huán)境模擬、城市規(guī)劃與災(zāi)害評(píng)估等場(chǎng)景。壓縮包共12個(gè)文件約219.69MB以tif柵格數(shù)據(jù)為核心配合dbf屬性表、prj坐標(biāo)系統(tǒng)、tfw世界文件、ovr與sbn/sbx金字塔緩存以及shp、shx等Shapefile組件兼顧高程讀取、坐標(biāo)定位與快速瀏覽。相比常見(jiàn)30米DEM12.5米精度能更細(xì)致刻畫(huà)山脊、山谷與坡度變化行政邊界則便于按區(qū)域裁剪統(tǒng)計(jì)。目前已有244人學(xué)習(xí)下載適合需要高精度地形底圖與邊界數(shù)據(jù)的中高級(jí)GIS用戶(hù)直接加載分析。1. 黔南州 12.5 米 DEM 拿到手之后先搞清楚這套數(shù)據(jù)能干什么黔南布依族苗族自治州位于貴州南部喀斯特地貌發(fā)育峰叢、洼地、河谷交錯(cuò)海拔落差從四百多米一直拉到近兩千米。這種地形條件下做水文分析、坡度坡向提取、工程選址或者三維場(chǎng)景搭建第一步都繞不開(kāi)一份靠譜的 DEM。12.5 米分辨率意味著每個(gè)柵格像元對(duì)應(yīng)地面約 156 平方米比常見(jiàn)的 30 米 DEM 精細(xì)一倍多能看清中小尺度的山脊線和溝谷走向又不像 5 米、2 米那樣動(dòng)輒幾十 GB 讓人望而卻步。標(biāo)題里這份數(shù)據(jù)是 7z 壓縮包里面除了 DEM 柵格還附帶了市級(jí)范圍的 shp 文件——這個(gè)組合很實(shí)用柵格管高程矢量管邊界裁剪拿到就能直接進(jìn) ArcGIS 或 QGIS 干活。適合做流域分析的水利從業(yè)者、做地形可視化的 GIS 工程師以及需要黔南州地形底圖做規(guī)劃或科研的人。但 7z 解壓、坐標(biāo)系確認(rèn)、DEM 與 shp 對(duì)齊這幾步每年都有大量人翻車(chē)下面把整條鏈路拆開(kāi)講。2. 7z 包解壓與數(shù)據(jù)完整性校驗(yàn)別讓第一步就卡住2.1 為什么這類(lèi)地理數(shù)據(jù)偏愛(ài) 7z 而不是 zip地理空間數(shù)據(jù)有個(gè)特點(diǎn)單個(gè) GeoTIFF 動(dòng)輒幾百 MB 到幾個(gè) GBshp 雖然小但往往配套 dbf、shx、prj 一堆文件。用 zip 壓縮壓縮率一般傳起來(lái)慢7z 用的是 LZMA2 算法對(duì)柵格這種有大量連續(xù)相同值的數(shù)據(jù)壓縮率能高出 30% 到 50%。所以國(guó)內(nèi)做數(shù)據(jù)分發(fā)的尤其是 DEM 這類(lèi)大柵格普遍選 7z。代價(jià)就是 Windows 自帶解壓不支持得裝 7-Zip 或者用支持 7z 的工具。另一個(gè)常見(jiàn)場(chǎng)景是分卷壓縮。黔南州全域 12.5 米 DEM 如果按整州拼一張?bào)w量不小分發(fā)時(shí)經(jīng)常切成.7z.001、.7z.002這種分卷。分卷包必須全部下載齊缺一個(gè)都解不開(kāi)而且要用 7-Zip 右鍵第一個(gè)分卷選「提取到當(dāng)前位置」不能單獨(dú)解壓中間某個(gè)卷。2.2 Windows 和 Linux 下的解壓命令Windows 圖形界面操作簡(jiǎn)單但做批量處理或者跑在服務(wù)器上命令行更靠譜。7-Zip 安裝后會(huì)把7z.exe加到路徑直接調(diào)用# Windows解壓到指定目錄-o 后面不能有空格 7z x 黔南州12.5米DEM.7z -oD:\qiannan_dem -y # 查看壓縮包內(nèi)容但不解壓先確認(rèn)里面有什么 7z l 黔南州12.5米DEM.7zLinux 下多數(shù)發(fā)行版?zhèn)}庫(kù)里有 p7zip裝完命令是7z或7za# Debian/Ubuntu 安裝 sudo apt install p7zip-full # 解壓-p 指定密碼如果有無(wú)密碼可省略 7z x qiannan_dem.7z -o./qiannan_dem -y # 分卷包只對(duì)第一個(gè)卷操作7z 會(huì)自動(dòng)找后續(xù)卷 7z x qiannan_dem.7z.001 -o./qiannan_dem參數(shù)說(shuō)明x是保留目錄結(jié)構(gòu)解壓e是把所有文件平鋪到當(dāng)前目錄地理數(shù)據(jù)千萬(wàn)別用e會(huì)把 tif 和 shp 混在一起-o指定輸出目錄注意-o和路徑之間沒(méi)有空格這是 7z 命令行的經(jīng)典坑-y表示全部確認(rèn)避免交互卡住腳本。2.3 解壓后先做三件事解壓完別急著往 GIS 里拖先確認(rèn)三件事。第一看文件清單里 DEM 是什么格式常見(jiàn)是.tif也可能是.img或.asc。第二看 shp 文件是否成套——一個(gè)完整的 shapefile 至少要有.shp、.shx、.dbf、.prj四個(gè)同名文件少一個(gè)都打不開(kāi)或者丟失坐標(biāo)系。第三看有沒(méi)有.prj文件它決定了坐標(biāo)系沒(méi)有它后面裁剪必然對(duì)不上。# 列出解壓目錄確認(rèn) shp 配套文件齊全 ls -lh ./qiannan_dem/ # 重點(diǎn)看有沒(méi)有 .shp .shx .dbf .prj 四件套如果 shp 缺了.prj可以用 ArcGIS 的「定義投影」補(bǔ)但前提是你知道原始坐標(biāo)系。國(guó)內(nèi)這類(lèi)行政邊界 shp 常見(jiàn)的是 CGCS2000 或 WGS84 地理坐標(biāo)系DEM 則可能是投影坐標(biāo)系比如 UTM 或高斯克呂格。兩者不一致時(shí)裁剪會(huì)直接報(bào)錯(cuò)或者結(jié)果錯(cuò)位這是后面避坑章要重點(diǎn)講的。3. DEM 與 shp 的坐標(biāo)系對(duì)齊裁剪能不能成全看這一步3.1 地理坐標(biāo)系和投影坐標(biāo)系到底差在哪很多人栽在坐標(biāo)系上是因?yàn)闆](méi)分清「地理坐標(biāo)系」和「投影坐標(biāo)系」。地理坐標(biāo)系GCS用經(jīng)緯度表示位置單位是度比如 CGCS2000、WGS84投影坐標(biāo)系PCS是把球面展平到平面上單位是米比如 UTM Zone 49N、CGCS2000 3 Degree GK Zone 35。DEM 做坡度、坡長(zhǎng)、匯流這些分析時(shí)必須用投影坐標(biāo)系因?yàn)椤付取箾](méi)法直接算距離和面積。而行政邊界 shp 經(jīng)常是地理坐標(biāo)系。判斷方法很簡(jiǎn)單打開(kāi)屬性看單位。如果顯示 Degree 就是地理坐標(biāo)系顯示 Meter 就是投影坐標(biāo)系。兩者不一致先統(tǒng)一再裁剪。3.2 用 QGIS 做坐標(biāo)系轉(zhuǎn)換和裁剪的最小流程QGIS 免費(fèi)、跨平臺(tái)處理這套數(shù)據(jù)足夠。假設(shè) DEM 是投影坐標(biāo)系shp 是地理坐標(biāo)系流程是先把 shp 轉(zhuǎn)成和 DEM 一致的投影再裁剪。# 用 GDAL 命令行做坐標(biāo)系轉(zhuǎn)換比開(kāi)圖形界面快 # 第一步查看 DEM 的坐標(biāo)系 gdalinfo qiannan_dem.tif | findstr AUTHORITY # 第二步把 shp 重投影到 DEM 的坐標(biāo)系 ogr2ogr -t_srs EPSG:32649 qiannan_boundary_proj.shp qiannan_boundary.shp # 第三步用重投影后的 shp 裁剪 DEM gdalwarp -cutline qiannan_boundary_proj.shp -crop_to_cutline \ -dstnodata -9999 qiannan_dem.tif qiannan_dem_clip.tif邏輯說(shuō)明ogr2ogr的-t_srs指定目標(biāo)坐標(biāo)系EPSG:32649 對(duì)應(yīng) UTM 49N覆蓋貴州南部gdalwarp的-cutline指定裁剪邊界-crop_to_cutline讓輸出范圍緊貼邊界而不是保留原圖幅-dstnodata把邊界外區(qū)域設(shè)為 -9999方便后續(xù)分析時(shí)排除。參數(shù)上-dstnodata的值要和 DEM 本身的無(wú)效值區(qū)分開(kāi)很多 DEM 默認(rèn)無(wú)效值是 -32768 或 -9999設(shè)重復(fù)了會(huì)混淆。3.3 裁剪后必須驗(yàn)證的兩件事裁剪完不是就完事了。第一檢查輸出范圍是否和 shp 邊界吻合可以在 QGIS 里疊加看也可以用gdalinfo看角點(diǎn)坐標(biāo)。第二檢查像元值有沒(méi)有異常比如邊界處出現(xiàn) 0 或者極大值那通常是 nodata 沒(méi)設(shè)對(duì)。# 查看裁剪結(jié)果的統(tǒng)計(jì)信息重點(diǎn)看最小最大值和 nodata gdalinfo -stats qiannan_dem_clip.tif如果最小最大值出現(xiàn) -9999 以外的異常值比如 0 或 65535說(shuō)明原始 DEM 的 nodata 定義和裁剪時(shí)設(shè)的不一致需要回去用gdal_translate -a_nodata重新指定。這一步不做后面做水文分析時(shí)洼地填充會(huì)把異常值當(dāng)真實(shí)高程結(jié)果全錯(cuò)。4. 避坑與排查黔南 DEM 處理里最容易翻車(chē)的 5 個(gè)點(diǎn)4.1 7z 解壓報(bào)錯(cuò)「文件損壞」但壓縮包明明沒(méi)問(wèn)題現(xiàn)象解壓到一半提示 CRC 校驗(yàn)失敗或者「不可預(yù)料的壓縮文件末端」。原因通常不是包壞了而是分卷沒(méi)下全或者下載過(guò)程中某個(gè)分卷字節(jié)數(shù)不對(duì)。黔南這種整州數(shù)據(jù)分卷常見(jiàn)少一個(gè).7z.00X就解不開(kāi)。解決核對(duì)分卷數(shù)量和文件大小重新下載缺失或異常的那個(gè)卷確保所有分卷在同一目錄下再對(duì)第一個(gè)卷操作。4.2 shp 拖進(jìn) ArcGIS 提示「缺少空間參考」現(xiàn)象shp 能打開(kāi)但坐標(biāo)系顯示 Unknown裁剪時(shí)提示無(wú)法投影。原因壓縮包里 shp 的.prj文件丟失或者分發(fā)時(shí)被過(guò)濾掉了。解決先確認(rèn)原始坐標(biāo)系——國(guó)內(nèi)行政邊界多為 CGCS2000 或 WGS84用 ArcGIS「定義投影」工具補(bǔ)上再和 DEM 統(tǒng)一。注意「定義投影」是給沒(méi)有坐標(biāo)系的數(shù)據(jù)貼標(biāo)簽「投影變換」才是真正轉(zhuǎn)換坐標(biāo)兩者別搞混。4.3 DEM 和 shp 疊加后邊界錯(cuò)位幾百米現(xiàn)象裁剪出來(lái)的 DEM 邊界和 shp 對(duì)不上整體偏移。原因兩者坐標(biāo)系不一致且偏移量正好是地理坐標(biāo)系和投影坐標(biāo)系之間的差異。解決用gdalinfo和ogrinfo分別查兩者坐標(biāo)系統(tǒng)一到同一個(gè)投影坐標(biāo)系再裁剪。如果偏移是幾十米級(jí)別可能是基準(zhǔn)面差異比如 WGS84 和 CGCS2000 之間需要做基準(zhǔn)面轉(zhuǎn)換而不是簡(jiǎn)單重投影。4.4 裁剪后 DEM 出現(xiàn)大片 0 值現(xiàn)象邊界內(nèi)出現(xiàn)成片 0 值看起來(lái)像空洞。原因原始 DEM 的 nodata 值被設(shè)成了 0裁剪時(shí)邊界外區(qū)域也用了 0兩者混在一起。解決用gdal_translate -a_nodata 0先把原始 DEM 的 nodata 改成 -9999再裁剪?;蛘卟眉魰r(shí)顯式指定-dstnodata -9999確保邊界外和原始無(wú)效值區(qū)分開(kāi)。4.5 Linux 下 7z 命令找不到現(xiàn)象7z: command not found。原因只裝了p7zip沒(méi)裝p7zip-full前者只有 7za 基礎(chǔ)功能不支持某些 7z 特性。解決sudo apt install p7zip-full裝完用7z而不是7za。如果還是找不到檢查 PATH或者用絕對(duì)路徑/usr/bin/7z。5. 從 DEM 到可用地形產(chǎn)品坡度、暈渲和流域提取的實(shí)操參數(shù)拿到裁剪好的黔南 DEM真正值錢(qián)的是從它派生出來(lái)的地形產(chǎn)品。這里挑三個(gè)最常用的坡度、山體暈渲、流域邊界把參數(shù)和驗(yàn)證方法講透。坡度提取在 ArcGIS 里用「坡度」工具輸出單位選 DegreeZ 因子是關(guān)鍵參數(shù)。黔南 DEM 如果是投影坐標(biāo)系單位是米Z 因子設(shè) 1如果是地理坐標(biāo)系單位是度Z 因子要設(shè)約 111320赤道處一度對(duì)應(yīng)的米數(shù)不設(shè)這個(gè)坡度會(huì)嚴(yán)重偏小。QGIS 里用「坡度」工具同理Z 因子在高級(jí)參數(shù)里。# 用 GDAL 生成坡度-p 表示輸出坡度而非坡向 gdaldem slope qiannan_dem_clip.tif qiannan_slope.tif -p -s 1.0 # 生成山體暈渲-az 方位角 -alt 高度角常用 315/45 gdaldem hillshade qiannan_dem_clip.tif qiannan_hillshade.tif \ -az 315 -alt 45 -z 1.0參數(shù)說(shuō)明-s 1.0是垂直 exaggeration 因子投影坐標(biāo)系下保持 1.0-az 315是光源方位角315 度西北方向是制圖慣例能讓山脊陰影自然-alt 45是光源高度角45 度適合大多數(shù)地形太低陰影過(guò)長(zhǎng)太高立體感弱。流域提取用 ArcGIS 的水文分析工具鏈填洼、流向、流量累積、柵格計(jì)算器提取河網(wǎng)、河流鏈接、流域分割。填洼閾值默認(rèn)就行黔南喀斯特地區(qū)洼地多填洼后能消除大量偽洼地。流量累積閾值決定河網(wǎng)密度12.5 米 DEM 建議從 500 開(kāi)始試太小河網(wǎng)碎太大主干河都提不出來(lái)。提取完用 shp 邊界裁剪流域再和實(shí)際水系圖對(duì)比驗(yàn)證。驗(yàn)證方法坡度圖疊加暈渲看紋理是否合理流域邊界疊加 DEM 看是否沿山脊線走。如果流域邊界橫切山脊說(shuō)明流向計(jì)算有問(wèn)題多半是填洼沒(méi)做或者 DEM 有異常值。我一般會(huì)先在小范圍試參數(shù)確認(rèn)無(wú)誤再跑全州避免跑幾小時(shí)出來(lái)發(fā)現(xiàn)參數(shù)錯(cuò)了。這套流程走通一次以后換任何地區(qū)的 DEM 都能復(fù)用希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取