簡介全國生態(tài)功能區(qū)劃修編版矢量數(shù)據(jù)資源包專為地理信息、生態(tài)評價、國土空間規(guī)劃、環(huán)境管理領(lǐng)域的科研人員和工程師準備幫助快速獲得標準化的全國生態(tài)功能分區(qū)邊界及屬性信息。壓縮包共含七個文件核心為shp矢量圖層并配套dbf屬性表、prj投影坐標、shx、sbn、sbx索引文件和xml元數(shù)據(jù)合計七點七二兆字節(jié)結(jié)構(gòu)緊湊可被ArcGIS、QGIS等主流軟件直接調(diào)用。該資源經(jīng)由作者整理分享當前已有三百一十三人學習或下載。在實際使用中用戶能夠直接開展疊加分析、屬性查詢和專題圖繪制用于生態(tài)紅線評估、環(huán)境承載力測算、區(qū)域開發(fā)適宜性評價等任務(wù)數(shù)據(jù)邊界與分類代碼源自修編版具備較好權(quán)威性免去自行找圖和矢量化流程大幅節(jié)省前期處理時間。對于需要全國尺度生態(tài)區(qū)劃數(shù)據(jù)支撐的研究、規(guī)劃與教學場景這是一份可直接入庫使用的基礎(chǔ)數(shù)據(jù)有助于快速搭建分析底圖、開展空間統(tǒng)計與決策輔助。1. 全國生態(tài)功能區(qū)劃 2015修編版shp tif 一份底圖生態(tài)評估與選線避讓的起點在哪里做環(huán)境影響識別、區(qū)域規(guī)劃里的生態(tài)紅線對接或者把一條擬建公路沿線的生態(tài)敏感類型拉出來分析手里沒有一份全國生態(tài)功能區(qū)劃 2015修編版底圖報告里很難給出讓人信服的空間依據(jù)。這份數(shù)據(jù)和常見路網(wǎng)、行政區(qū)劃不一樣發(fā)布形態(tài)一般是 shp 矢量邊界配合 tif 柵格shp 用來查功能區(qū)名稱、等級和代碼tif 像元值則可以直接參與柵格計算、面積統(tǒng)計以及和 DEM、土地利用這類數(shù)據(jù)做疊置分析。兩套格式看著互補實際用起來卻往往因為坐標系不統(tǒng)一、NoData 沒處理、柵格和面邊界對不齊在一開始就把人卡住。下面按我實際操作的順序把讀入、預(yù)處理、裁剪、轉(zhuǎn)換和避坑的完整流程梳理一遍照著走能少花很多前期摸索的時間。2. 數(shù)據(jù)讀入與預(yù)處理先搞清 shp 字段和 tif 值域再做坐標系對齊2.1 屬性表里查生態(tài)功能區(qū)常見字段、區(qū)劃代碼與名稱的對應(yīng)不同渠道拿到的數(shù)據(jù)屬性字段命名經(jīng)常不統(tǒng)一我用過的版本里既有叫 code、name、type 的也有直接用拼音縮寫如 xzq行政區(qū)、stgnq生態(tài)功能區(qū)的。不管字段長什么樣建議進 ArcMap 后用 Python 窗口先把字段名全部打出來再確認區(qū)劃代碼到底存在哪個字段里。# ArcMap Python 窗口 import arcpy fc rD:\eco\eco_function_area.shp fields [f.name for f in arcpy.ListFields(fc)] print(fields) with arcpy.da.SearchCursor(fc, fields[:6]) as cur: for i, row in enumerate(cur): if i 6: break print([str(v)[:30] for v in row])這段代碼用 arcpy.ListFields 把屬性字段全部列出來再用 SearchCursor 讀前 6 條記錄的字段值。這樣能快速判斷 shp 里到底有沒有“區(qū)劃名稱”“生態(tài)功能類型”“所屬省份”這類信息避免做圖例時到處找字段。打印時把字符串截到 30 位只是為了讓輸出整齊不影響源數(shù)據(jù)本身?？炊侄魏笾攸c看兩個關(guān)鍵項區(qū)劃代碼和功能類型。2015 修編版的區(qū)劃體系是“生態(tài)功能區(qū)—生態(tài)功能亞區(qū)—生態(tài)功能小區(qū)”三個層級但屬性表不一定每層都列全。大多數(shù)發(fā)布版本會有一個漢字字段寫一級功能類型比如水源涵養(yǎng)、水土保持、防風固沙、生物多樣性保護等而 tif 柵格里的像元值往往是針對這套類型做了數(shù)值映射。所以先把代碼和名稱的對照關(guān)系記下來后面做柵格轉(zhuǎn)面、按類型統(tǒng)計面積全靠這層映射不出錯?？臻g參考也建議在圖層屬性里看一眼。有些渠道下載的華東、華北分幅數(shù)據(jù)自帶投影坐標有些直接是 WGS84 地理坐標。拿北京地區(qū)舉例如果 shp 是 CGCS2000 下的高斯投影帶tif 卻是 WGS84兩者疊加會整體錯位幾米到幾十米必須先統(tǒng)一。常見做法是全部轉(zhuǎn)到 CGCS2000或按你最終制圖要求轉(zhuǎn)到所在區(qū)域的高斯分帶然后矢量做投影、柵格做重采樣一步到位。2.2 柵格 tif 的像元值、NoData 和空值處理工具怎么選打開 tif 后第一件事不是看顏色而是看屬性里的像元大小、波段數(shù)和 NoData 值。生態(tài)功能區(qū)劃 tif 如果是單波段整數(shù)型像元值通常直接對應(yīng)區(qū)劃代碼如果包含多波段通常第一波段是區(qū)劃類型后續(xù)波段可能存放生態(tài)敏感性等級等輔助信息。不少發(fā)布版本為了壓縮體積把 NoData 設(shè)成 -9999 或 0這兩個值一旦參與面積統(tǒng)計結(jié)果會多出一大片“未知類型”區(qū)域。提示tif 加載后整幅圖是同一個色多半不是數(shù)據(jù)壞了而是符號化用了連續(xù)拉伸不適合整數(shù)類型的分類數(shù)據(jù)。改成“唯一值”分類符號化區(qū)劃界線立刻就能顯示出來。處理 NoData我常用的工具是柵格計算器里的 SetNull 和 IsNull也可以在 ArcToolbox 里找“條件函數(shù)”完成。下面表達式等價于一個條件賦值# ArcGIS 柵格計算器表達式 out SetNull(IsNull(eco_tif) | (eco_tif -9999), eco_tif)這句的意思是當像元為 NoData或像元值等于 -9999 時把該位置設(shè)為空其余像元保留原值。IsNull 先生成一個 0/1 掩膜SetNull 再根據(jù)條件把對應(yīng)像元置空。之所以必須處理 -9999是因為很多軟件寫柵格時把“數(shù)據(jù)缺失”寫成了 -9999 而不是標準 NoData不轉(zhuǎn)換的話后面轉(zhuǎn)面、統(tǒng)計會把 -9999 當成一個真實的區(qū)劃類型生成整片無效圖斑。用 QGIS 或 GDAL 命令行等價操作是gdal_calc.py -A eco_2015.tif --outfileeco_clean.tif \ --calcA*(A0) --NoDataValue0這條命令把小于等于 0 的像元全部寫成 0 并標記為 NoData。執(zhí)行前先看一眼直方圖確認 tif 的值是連續(xù)的短代碼如 101、102、201還是大數(shù)字編號如 1000001避免條件判斷寫錯把有效類型一起抹掉。2.3 坐標系與像元對齊投影轉(zhuǎn)換、重采樣方法和 Snap Rastershp 和 tif 都加載后如果兩者“看起來”重疊但相交面積統(tǒng)計出來卻少了一截多半是投影基準和像元對齊的問題。常見做法是先把 tif 投影到和 shp 同一個坐標系隨后在環(huán)境設(shè)置里把 Snap Raster 設(shè)為生態(tài)區(qū)劃 tifCell Size 也固定成它的像元大小。投影柵格時重采樣方法必須選“最近鄰”NEAREST。生態(tài)功能區(qū)劃是離散類型數(shù)據(jù)雙線性或三次卷積會對邊界像元做插值比如說 101 和 102 之間插出 101.5這個值轉(zhuǎn)面后就成了“未知類型”。而最近鄰法只取原始值不會生成新類別。矢量 shp 的投影用 Project 工具目標坐標系選 CGCS2000 的地區(qū)分帶或 Albers 等面積投影都可以關(guān)鍵是 shp 和 tif 最終必須在同一套坐標系里。環(huán)境設(shè)置往往是新手最容易忽略的一步。在 ArcMap 的“環(huán)境設(shè)置”中把處理范圍設(shè)為生態(tài)區(qū)劃 shp 的范圍捕捉柵格設(shè)為生態(tài)區(qū)劃 tif像元大小填 tif 的像元尺寸。這樣后續(xù)任何柵格運算的輸出都和源 tif 網(wǎng)格嚴格對齊不會出現(xiàn)錯半格的情況。如果輸出偏差半個像元邊界會呈現(xiàn)鋸齒狀轉(zhuǎn)面后也很容易產(chǎn)生細碎窄條多邊形。對齊完成后可以用下面幾句檢查兩個數(shù)據(jù)的范圍是否一致import arcpy eco_tif rD:\eco_2015.tif eco_shp rD:\eco_function_area.shp r1 arcpy.Describe(eco_tif).extent r2 arcpy.Describe(eco_shp).extent print(r1.XMin, r1.YMin, r1.XMax, r1.YMax) print(r2.XMin, r2.YMin, r2.XMax, r2.YMax)前后兩行范圍值基本重合說明后續(xù)裁剪、提取都能正常工作差別明顯時以 shp 面范圍為準做裁剪不要反過來用 tif 范圍去約束 shp否則容易把邊界切掉一圈。3. 用生態(tài)區(qū)劃 shp 裁剪 DEM 柵格Clip 與掩膜提取的差異和參數(shù)選擇3.1 Clip 和 Extract by Mask 在輸出邊界上的真實差異很多人搜“arcmap 中依靠面圖層裁剪 dem 柵格 tif 文件”和“依靠面圖層掩碼提取”分不清兩個工具到底有何區(qū)別。ArcGIS 的柵格 Clip數(shù)據(jù)管理可以按矩形范圍裁剪也可以勾選“使用輸入要素裁剪幾何”選項使輸出范圍貼合面要素邊界。Extract by Mask Spatial Analyst則是把掩膜面柵格化后保留掩膜內(nèi)部的像元掩膜外一律寫成 NoData。兩者從結(jié)果看經(jīng)常很接近但有三個實際差異NoData 策略不同。Clip 在裁剪范圍內(nèi)的 NoData 原樣保留Extract by Mask 會把掩膜邊界外全部設(shè)為 NoData輸出的有效像元范圍看起來更“干凈”。表達式能力不同。Extract by Mask 配合柵格計算器可以做帶條件的提取比如只保留生態(tài)類型等于“水源涵養(yǎng)”的像元Clip 沒有表達式選項只能按幾何切。性能不同。數(shù)據(jù)量大的時候Clip 更快因為它不做掩膜重分類直接按范圍切塊Extract by Mask 多了掩膜柵格化這一步速度略慢。所以我的選擇策略是只做幾何范圍切割用 Clip簡單直接后續(xù)還要按類型疊加統(tǒng)計用 Extract by Mask 更順手。但這里有個高頻坑之前的工程在環(huán)境設(shè)置里掛了“分析掩膜”或“捕捉柵格”ArcMap 會默默繼承這些設(shè)置導(dǎo)致明明選的全國范圍輸出還是別處的矩形塊。跑之前先打開環(huán)境設(shè)置把不相關(guān)的掩膜清掉。3.2 跟著操作ArcMap 中用生態(tài)區(qū)劃面裁剪 DEM 的掩膜提取步驟下面這套是我在 ArcMap 里最常用的流程前提是數(shù)據(jù)已經(jīng)按第 2 章做完坐標系和 NoData 處理。我之前拿這套方案把全國生態(tài)功能區(qū)劃里的“水源涵養(yǎng)”區(qū)單獨提取出來再接上 30 米 DEM 算每個子區(qū)域的平均坡度。import arcpy arcpy.env.workspace rD:\eco_work arcpy.env.extent rD:\eco_work\eco_function_area.shp arcpy.env.snapRaster rD:\eco_work\eco_tif.tif arcpy.env.cellSize rD:\eco_work\eco_tif.tif arcpy.env.mask rD:\eco_work\eco_function_area.shp out arcpy.sa.ExtractByMask(rD:\dem_30m.tif, rD:\eco_work\eco_function_area.shp) out.save(rD:\eco_work\dem_eco_clip.tif)這段代碼把處理范圍、捕捉柵格、像元大小和分析掩膜統(tǒng)一設(shè)到生態(tài)功能區(qū) shp 上。ExtractByMask 的第二個參數(shù)可以是矢量面工具內(nèi)部會先把它柵格化。注意 Spatial Analyst 需要啟用擴展模塊否則會直接報“工具不可用”這也是個高頻入門問題。環(huán)境參數(shù)里最有玄機的是 snapRaster。它保證輸出像元和生態(tài)區(qū) tif 的網(wǎng)格完全重合否則會把柵格重新對齊到當前坐標系默認網(wǎng)格邊界像元錯開半個像元。對于生態(tài)功能區(qū)的面積統(tǒng)計這半個像元在邊界上累計起來能差出幾百上千平方米。在環(huán)評報告里“面積對不上”往往會被評審盯住所以這個參數(shù)我每次都會顯式設(shè)置。如果手頭沒有 DEM這一步也可以直接裁生態(tài)區(qū) tif 本身邏輯完全一樣。很多人下載了全國 tif 后先按省份邊界裁成省圖再按市、縣 shp 繼續(xù)切分這樣后續(xù)每次加載不用扛著全國范圍的大文件效率明顯更好。3.3 驗證裁剪結(jié)果像元數(shù)、唯一值和邊界貼合度裁剪完成后不要直接拿去制圖先做三個快速檢查確認輸出像元大小沒被重采樣改掉用唯一值符號化確認值集合是原 tif 的子集疊加 shp 看邊界是否貼合。對應(yīng)命令行檢查import arcpy from arcpy.sa import * r Raster(rD:\eco_work\dem_eco_clip.tif) print(r.width, r.height, r.cellSize) uc arcpy.UniqueValues(r) print(uc[:20])這段代碼打印輸出柵格的寬度、高度、像元大小以及前 20 個唯一值。看到 -9999 或 0 混在列表里就回頭看掩膜和 NoData 是否處理干凈。更直觀的辦法是把生態(tài)區(qū)劃 shp 和裁出來的 tif 疊加在 ArcMap 里沿著邊界放大到 1:5000 左右肉眼就能看出有沒有錯開。不要嫌這一步啰嗦數(shù)據(jù)問題大多出在源數(shù)據(jù)本身把驗證步驟固定下來每次跑都少踩很多坑。4. 生態(tài)區(qū)劃數(shù)據(jù)處理排查與避坑五個高頻問題的現(xiàn)象、原因、解決4.1 現(xiàn)象tif 加載后整幅圖灰黑一片看不到任何區(qū)劃圖斑原因出在符號化。生態(tài)功能區(qū)劃是離散分類柵格很多軟件默認按連續(xù)拉伸顯示把所有像元值壓到一個灰度區(qū)間看起來就是一塊灰。有些版本還把 NoData 區(qū)域和有效區(qū)域混在一起有效圖斑占比又低符號化后幾乎看不見。解決方法是先把符號化方式從“拉伸”改成“唯一值”并給不同功能類型分配不同色帶。如果改成唯一值后仍是一片黑用第 2 章的 SetNull 表達式把 NoData 和 -9999 轉(zhuǎn)成空再重新符號化。絕大多數(shù)“數(shù)據(jù)壞了”的錯覺到這里就好了。4.2 現(xiàn)象shp 與 tif 看似重疊邊界卻總有一圈錯位現(xiàn)象是疊加后大輪廓基本對得上但放大看 shp 邊界和 tif 的像元邊界差出幾個像元像一圈貼邊。原因通常有兩個一個坐標系不是同一套基準比如 WGS84 和 CGCS2000 之間的差異在高精度要求下會表現(xiàn)出來另一個是 tif 經(jīng)過某次重采樣像元網(wǎng)格和 shp 的投影網(wǎng)格不再對齊。解決方法是把兩者統(tǒng)一投影到 CGCS2000 的同一分帶并在環(huán)境設(shè)置里把 Snap Raster 指向 tif。檢查辦法是在 ArcMap 里把 shp 和 tif 都打開用“按屬性選擇”選中任意一個圖斑放大后對比邊界。如果錯位均勻且都在 10 米以內(nèi)通常就是投影差異如果錯位方向不一致還要懷疑 shp 本身是否做過局部編輯這種只能對照原始發(fā)布數(shù)據(jù)重新下載。4.3 現(xiàn)象裁剪輸出邊緣鋸齒明顯柵格轉(zhuǎn)面后出現(xiàn)大量細碎多邊形裁剪之后邊緣呈明顯階梯狀轉(zhuǎn)成 shp 后沿邊界出現(xiàn)一串長條狀多邊形這種問題多半是環(huán)境設(shè)置里沒有指定 Snap Raster 和 Cell Size系統(tǒng)按當前數(shù)據(jù)框的默認網(wǎng)格重新對齊了一次。生態(tài)區(qū)劃 tif 的像元原本是正南正北方向排列的重對齊后像元邊界就歪了。解決方法是重新執(zhí)行裁剪在環(huán)境設(shè)置里把捕捉柵格設(shè)為生態(tài)區(qū)劃 tifCell Size 填 tif 的原始像元尺寸。轉(zhuǎn)面時如果碎多邊形已經(jīng)產(chǎn)生會用“消除”工具合并但最穩(wěn)妥的辦法還是在裁剪階段就把網(wǎng)格對齊磨刀不誤砍柴工。4.4 現(xiàn)象柵格轉(zhuǎn)面后屬性表里的代碼全是 -9999類型字段為空柵格轉(zhuǎn)面工具Raster to Polygon在沒有先處理 NoData 的情況下會把 -9999 或 0 也當成一類像元轉(zhuǎn)出來的面屬性里出現(xiàn)大量“-9999”記錄真正的區(qū)劃代碼和名稱反而查不到。原因是 NoData 在轉(zhuǎn)面時被當作普通類值參與建面。解決方法是先按第 2 章做法把 -9999 設(shè)為 NoData或者用柵格計算器把非目標范圍的像元賦空值再執(zhí)行轉(zhuǎn)面。另外一個常見習慣是直接把 tif 用“按掩膜提取”先裁到研究區(qū)范圍內(nèi)邊角就不會有 -9999 混進來。轉(zhuǎn)面前順手跑一次唯一值統(tǒng)計看到列表里只有目標代碼再轉(zhuǎn)。4.5 現(xiàn)象用漁網(wǎng)工具分割 shp 后各部分面積總和比原始面小用創(chuàng)建漁網(wǎng)工具把生態(tài)功能區(qū)分成網(wǎng)格再逐格裁剪最后匯總面積卻發(fā)現(xiàn)少了。原因是漁網(wǎng)是獨立生成的規(guī)則格網(wǎng)它的邊界不可能和原始面邊界完全重合。裁剪時漁網(wǎng)格與面邊界重疊的窄條被丟掉共同產(chǎn)生的邊緣損失累計后相當可觀。解決方法是不要用漁網(wǎng)直接裁剪原始面而應(yīng)該用“相交”工具讓原始面與漁網(wǎng)線做空間相交相交后每個格網(wǎng)部分保留原面屬性面積由相交結(jié)果重新計算。這樣所有面積都遵循原始邊界不會出現(xiàn)丟邊。硬要用裁剪也是一條路但必須接受邊緣損失并在匯總時單獨說明。5. 讓生態(tài)區(qū)劃柵格快速復(fù)用柵格轉(zhuǎn)面、歸一化處理和導(dǎo)出 WKT、3D Tiles生態(tài)區(qū)劃 shp 和 tif 的價值往往不在全國圖本身而在把它轉(zhuǎn)成可復(fù)用的工程數(shù)據(jù)。我常用的三條路徑是柵格轉(zhuǎn)面、歸一化處理和轉(zhuǎn)成通用文本格式。柵格轉(zhuǎn)面用來把 tif 的離散類型變成矢量圖斑方便和其他面圖層做空間連接歸一化則解決多源柵格之間的數(shù)值尺度問題。歸一化處理我通常用柵格計算器表達式是(A - min) / (max - min)把生態(tài)區(qū)劃代碼壓到 0 到 1 之間。這對類型型數(shù)據(jù)而言只適合做顯示真正做疊加分析時還是保留原始代碼更可靠別為了統(tǒng)一尺度把分類語義丟掉。導(dǎo)出 WKT 可以用 GDAL 的 ogr2ogrogr2ogr -f CSV eco_wkt.csv eco_function_area.shp -dialect sqlite \ -sql SELECT code, name, ST_AsText(geometry) AS geom FROM eco_function_area這條命令把 shp 的每個面要素轉(zhuǎn)成一行 WKT 文本適合做數(shù)據(jù)庫入庫或程序?qū)?。如果目標?3D Tiles常見做法是先用轉(zhuǎn)面工具把 tif 轉(zhuǎn)成 shp再通過其他工具鏈做瓦片化這一步并非 ArcMap 自帶需要單獨搭建。最后說一句我的個人習慣每次拿到生態(tài)區(qū)劃數(shù)據(jù)第一件事永遠是復(fù)制一份原始文件然后在副本上進坐標系、NoData 和投影轉(zhuǎn)換絕不直接動下載源文件。柵格數(shù)據(jù)破壞性操作沒有后悔藥寧可多占一點硬盤也別在原始底圖上反復(fù)折騰。這套流程幫我避開了很多沒法回退的局面希望也能幫到你。本文還有配套的精品資源點擊獲取