據(jù)吉林省10米分辨率裁剪投影與面積統(tǒng)計(jì)避坑指南)
簡介2021年土地利用ESRI吉林省10m精度數(shù)據(jù)集由全球GIS領(lǐng)域頭部廠商ESRI官方制作采用WGS84坐標(biāo)系統(tǒng)與10米柵格分辨率完整呈現(xiàn)吉林省土地覆蓋狀況。數(shù)據(jù)源自海量全球原始柵格經(jīng)地理配準(zhǔn)、裁剪、轉(zhuǎn)換等處理按長春、吉林、四平、遼源、通化、白山、松原、白城、延邊9個地級行政區(qū)劃分獨(dú)立目錄壓縮包共63個文件大小約44MB核心為TIF格式土地覆蓋柵格圖配套TFW坐標(biāo)定位文件、DBF與CPG屬性表、XLSX分類統(tǒng)計(jì)表、XML元數(shù)據(jù)及PNG快速預(yù)覽圖。借助主流GIS軟件可提取農(nóng)田、森林、水域、城市建筑、草地等用地類型完成面積統(tǒng)計(jì)、變化檢測與空間制圖。數(shù)據(jù)適用于土地覆蓋變化監(jiān)測、生態(tài)環(huán)境評估、城鄉(xiāng)規(guī)劃、農(nóng)業(yè)管理、自然災(zāi)害風(fēng)險評估及氣候變化研究等多類任務(wù)。該數(shù)據(jù)集已有258人學(xué)習(xí)下載是GIS從業(yè)人員、規(guī)劃決策者與高校師生開展省級尺度空間分析的實(shí)用基礎(chǔ)數(shù)據(jù)。1. 2021 年 ESRI 土地利用數(shù)據(jù)吉林省 10 米分辨率到底能用來干什么很多做東北項(xiàng)目的同事電腦里還躺著 GlobeLand30 那套 30 米地表覆蓋理由就一句話“省級分析夠用了”。但 2021 年 ESRI 發(fā)布的這套全球 10 米土地利用數(shù)據(jù)把空間分辨率抬到了 Sentinel-2 能做到的極限吉林省作為長白山區(qū)和松嫩平原的交錯帶最吃這套數(shù)據(jù)的紅利山區(qū)能撕開溝谷和耕地邊界平原能大致分清村落、農(nóng)田和鹽堿地。不過它原始坐標(biāo)系是 WGS84 地理坐標(biāo)直接拿原圖算面積會翻車而且分類體系只有 9 類跟國內(nèi)常用的土地利用二級類對不上。這篇文章就圍繞“ESRI 2021 土地利用、吉林省、10 米精度”講清楚它是什么、怎么裁剪投影、怎么統(tǒng)計(jì)面積以及五個實(shí)際踩過的坑。適合 GIS 從業(yè)者、做國土空間規(guī)劃和農(nóng)業(yè)遙感監(jiān)測的人直接照著做。2. 分類體系與數(shù)據(jù)原理9 類不夠用但每一類都有明確含義2.1 九個類別字段每類對應(yīng)什么地表特征這套數(shù)據(jù)在 ArcGIS Living Atlas 里的全稱是 Sentinel-2 10m Land Use/Land Cover2021 年度合成版本使用 9 類體系。先看編碼后面所有面積統(tǒng)計(jì)都依賴這張表編碼英文類名中文對應(yīng)在吉林省的典型表現(xiàn)1Water水體松花江、查干湖、大型水庫2Trees樹木長白山針闊混交林、次生林3Grass草地西部草場、林間隙地4Flooded Vegetation被淹植被河漫灘、濕地、洪泛區(qū)植被5Crops作物玉米、水稻、大豆田塊6Built Area建成區(qū)長春、吉林市區(qū)及鄉(xiāng)鎮(zhèn)聚落7Bare Ground裸地白城鹽堿地、裸土、采石場8Snow/Ice冰雪冬季積雪也常造成誤分9Clouds云無意義類別統(tǒng)計(jì)前應(yīng)剔除注意第 4 類“被淹植被”在官方文檔里對應(yīng)的是紅樹林、洪泛區(qū)、沼澤這類長期或季節(jié)性淹水的植被并不是水稻田。但吉林省水稻移栽期田間有明水光譜特征和 Flooded Vegetation 很像后面避坑章會專門講這個事。2.2 物候合成與深度學(xué)習(xí)的生產(chǎn)邏輯ESRI 這套產(chǎn)品不是拿單景影像做的分類而是用 2021 年全年 Sentinel-2 L2A 數(shù)據(jù)按月合成每個像元保留最能代表該月地表狀態(tài)的反射率值再把 12 個月的合成影像輸入到一個深度學(xué)習(xí)語義分割模型里逐像元輸出類別標(biāo)簽。這樣做的直接結(jié)果是常綠針葉林和落葉闊葉林在年度合成里會被歸到同一個“Trees”類這在東北是能接受的畢竟林業(yè)二類調(diào)查用的是亞類而土地利用現(xiàn)狀調(diào)查要的是喬木林地這個一級類口徑。和 30 米 GlobeLand30 相比10 米產(chǎn)品在吉林省的實(shí)際差異非常直觀。第一是線性地物鄉(xiāng)間道路、河渠、防護(hù)林帶在 30 米像元里基本糊成混合像元10 米能單獨(dú)成行。第二是聚落邊界長春周邊城中村、開發(fā)區(qū)在建工地10 米分辨率下建成區(qū)和裸地的邊界干凈很多。第三是耕地破碎度吉林東部半山區(qū)的小塊農(nóng)田30 米產(chǎn)品經(jīng)常漏分10 米分出來的地塊更接近真實(shí)田塊邊界。2.3 為什么說它不能直接對標(biāo)國土三調(diào)國內(nèi)用戶最常問的問題這套數(shù)據(jù)和國土變更調(diào)查里的地類圖斑能不能直接互相替換答案是不能。三調(diào)的“水田”要求結(jié)合水利設(shè)施和種植制度認(rèn)定ESRI 的 Crops 類只反映 2021 年影像上“像作物的植被”兩者口徑不同。這套數(shù)據(jù)真正適合的場景是宏觀的土地利用變化監(jiān)測、生態(tài)評估、國土空間規(guī)劃的前期現(xiàn)狀摸底而不是地籍級的確權(quán)應(yīng)用。把它定位成“省級宏觀現(xiàn)狀速覽圖”最合適精度優(yōu)勢體現(xiàn)在空間表達(dá)上不在分類語義深度上。3. 從 Living Atlas 到吉林省 GeoTIFF三種獲取路徑與導(dǎo)出參數(shù)3.1 ArcGIS Pro 動態(tài)影像服務(wù)導(dǎo)出最穩(wěn)妥的獲取方式是直接在 ArcGIS Pro 里連 Living Atlas 的影像服務(wù)。打開目錄窗格在 Portal 選項(xiàng)卡下瀏覽 Living Atlas搜索 Sentinel-2 10m Land Use/Land Cover把它拖到地圖視圖。注意要選 2021 年份的圖層這個產(chǎn)品按年份拆成了多條服務(wù)混用年份會導(dǎo)致后續(xù)分析出錯。確認(rèn)圖層加載后通過“數(shù)據(jù)管理 → 導(dǎo)出柵格”或在圖層右鍵菜單里選“導(dǎo)出 → 為不同格式”。關(guān)鍵參數(shù)如下Extent范圍選擇“繪制范圍”并加載吉林省邊界矢量或用“環(huán)境 → 處理范圍”指定邊界要素類。Cell Size設(shè)為 10 米。Resampling必須選 Nearest Neighbor這是分類柵格雙線性插值會在類別邊緣產(chǎn)生不存在的地類編碼。FormatTIFF。CompressionLZW 或 DEFLATE10 米整省數(shù)據(jù)壓縮后體積能小很多。整省范圍直接在動態(tài)服務(wù)上導(dǎo)出經(jīng)常遇到服務(wù)端超時。我一般會先把吉林省邊界外擴(kuò) 1 公里作為導(dǎo)出范圍然后把目標(biāo)拆成東部山區(qū)、中部平原、西部草原三塊分別導(dǎo)出最后再用鑲嵌工具合并。動態(tài)服務(wù)的好處是不用關(guān)心原始分塊和你有沒有 1T 硬盤壞處是受服務(wù)性能限制大范圍導(dǎo)出容易斷。3.2 官方分塊 GeoTIFF 與 QGIS 加載ESRI 在 Living Atlas 的數(shù)據(jù)詳情頁還提供了按 100km × 100km 分塊的年度 GeoTIFF 直接下載入口每個分塊文件名包含它的經(jīng)緯度網(wǎng)格編號。吉林省橫跨大約 122°E 到 131°E、41°N 到 46°N需要匹配 4 到 6 個分塊文件下載后先鑲嵌成一張全圖再做裁剪。這個路徑適合網(wǎng)絡(luò)條件好、想離線分析的用戶。如果你用 QGIS動態(tài)影像服務(wù)同樣可以加載圖層 → 添加圖層 → ArcGIS REST Server粘貼服務(wù)的 REST URL。但沒有 ArcGIS Pro 的“導(dǎo)出柵格”那么順手需要先加載到畫布再用“柵格 → 轉(zhuǎn)換 → 裁剪”按邊界矢量導(dǎo)出速度比 Pro 慢不少內(nèi)存也吃得更緊。3.3 導(dǎo)出后先做三個檢查下載完成先別急著投影和統(tǒng)計(jì)花兩分鐘做三個檢查能省掉后面很多排查時間。第一檢查 NoData 設(shè)置。動態(tài)服務(wù)導(dǎo)出的 GeoTIFF 在省界外通常是有值的只是顯示為黑色背景要確認(rèn) NoData 有沒有正確寫入。第二檢查類別直方圖用 ArcGIS 的“查看屬性表”或 QGIS 的“唯一值”統(tǒng)計(jì)看看 9Clouds類的像元數(shù)量是否異常多如果超過總像元數(shù)的 5%說明 2021 年該區(qū)域云影響嚴(yán)重需要重新評估這個年份的質(zhì)量。第三檢查空間參考確認(rèn)文件仍是 WGS84 地理坐標(biāo)系如果已經(jīng)被自動轉(zhuǎn)成了 Web Mercator先退回原坐標(biāo)系再走下一步。4. 裁剪、投影與面積統(tǒng)計(jì)可抄作業(yè)的 GDAL 參數(shù)4.1 用邊界矢量做精確裁剪拿到吉林省范圍的 GeoTIFF 后第一步是裁剪。推薦用 GDAL 的 gdalwarp 一次性完成裁剪加范圍收緊命令行如下gdalwarp -overwrite \ -cutline jilin_2021.shp \ -crop_to_cutline \ -dstnodata 255 \ -r near \ -of GTiff \ LULC_2021_jilin_raw.tif \ LULC_2021_jilin_wgs84.tif-cutline指定吉林省邊界矢量-crop_to_cutline讓輸出柵格的范圍嚴(yán)格貼合邊界省界外不再保留多余的矩形區(qū)域。-dstnodata 255把邊界外的像元設(shè)為 255這樣后面統(tǒng)計(jì)時可以直接排除。-r near保持最近鄰重采樣分類柵格不能做平滑插值。如果你手里的原始數(shù)據(jù)是多個分塊先執(zhí)行g(shù)dalbuildvrt LULC_2021_jilin.vrt block1.tif block2.tif block3.tif block4.tif gdal_translate -of GTiff LULC_2021_jilin.vrt LULC_2021_jilin_raw.tifVRT 虛擬柵格的好處是零拷貝合并不產(chǎn)生中間文件適合分塊較多的情況。合并完成后還有一步很關(guān)鍵用 gdalinfo 查看輸出文件的尺寸確認(rèn)沒有因?yàn)榉謮K重疊產(chǎn)生坐標(biāo)偏移。4.2 為什么用 Albers 等積投影而不是三度帶吉林省跨越東經(jīng) 121 度到 131 度橫跨多個三度帶。很多習(xí)慣用國家 2000 三度帶投影的同行會按帶號分帶投影最后再拼接。這個做法在縣一級沒問題但整省范圍會出現(xiàn)兩個麻煩一是分帶接邊處地物錯位二是兩個帶分別做面積統(tǒng)計(jì)后再合并類別面積在拼接帶附近會有系統(tǒng)性偏差。做土地利用面積統(tǒng)計(jì)第一原則是投影必須等積。我一般用 Albers 等面積圓錐投影參數(shù)按東北地區(qū)常用設(shè)置gdalwarp -overwrite \ -t_srs projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 ellpsGRS80 unitsm no_defs \ -r near \ -tr 10 10 \ LULC_2021_jilin_wgs84.tif \ LULC_2021_jilin_albers.tif投影參數(shù)里lat_125和lat_247是雙標(biāo)準(zhǔn)緯線吉林省主體緯度 41 度到 46 度正好被夾在中間變形最小。lon_0105是中央經(jīng)線這是全國標(biāo)準(zhǔn) Albers 的常用值如果你只做吉林省改成 126 度會讓圖面中央變形更小但面積統(tǒng)計(jì)結(jié)果和 105 度幾乎沒有差別。-tr 10 10強(qiáng)制輸出像元為 10 米乘 10 米的正方形這一步直接決定后續(xù)面積統(tǒng)計(jì)的精度——投影后每個像元對應(yīng) 100 平方米。4.3 分塊統(tǒng)計(jì)各類面積避免一次性讀入內(nèi)存面積統(tǒng)計(jì)的思路很簡單統(tǒng)計(jì)每個類別編碼的像元數(shù)乘以單個像元面積。但吉林省全境 10 米分辨率柵格寬度約 74000 個像元高度約 62000 個像元總像元數(shù)超過 45 億。直接用 GDAL 的 ReadAsArray 把全圖讀進(jìn)內(nèi)存uint8 數(shù)據(jù)也要超過 4.5GB普通工作站會直接內(nèi)存溢出。正確做法是分塊讀取from osgeo import gdal from collections import Counter import numpy as np ds gdal.Open(LULC_2021_jilin_albers.tif, 0) band ds.GetRasterBand(1) xsize band.XSize ysize band.YSize # 讀取仿射變換得到單像元的寬度和高度乘起來就是單像元面積 gt ds.GetGeoTransform() px_area abs(gt[1] * gt[5]) print(f單像元面積: {px_area:.1f} 平方米) # 分塊大小取 2048 行寬度保持全幅單次內(nèi)存約 74000*2048*1 字節(jié) ≈ 150MB block_size 2048 counter Counter() for row in range(0, ysize, block_size): rows min(block_size, ysize - row) data band.ReadAsArray(0, row, xsize, rows) vals, counts np.unique(data, return_countsTrue) for v, c in zip(vals, counts): if v not in (255, 9): # 剔除邊界 NoData 和云類 counter[v] c print(f已處理第 {row}/{ysize} 行) # 輸出各面積平方米轉(zhuǎn)公頃 class_names {1: 水體, 2: 樹木, 3: 草地, 4: 被淹植被, 5: 作物, 6: 建成區(qū), 7: 裸地, 8: 冰雪} for cls_id in sorted(class_names): ha counter[cls_id] * px_area / 10000 print(f{class_names[cls_id]}: {ha:.1f} 公頃)px_area的值在 Albers 投影下就是 100.0 平方米它的正確性是面積統(tǒng)計(jì)可信的基礎(chǔ)。block_size可以按機(jī)器內(nèi)存調(diào)整內(nèi)存 16GB 以上的機(jī)器可以設(shè) 40968GB 的機(jī)器建議保持 2048。輸出結(jié)果后再在 Excel 里手動合并類別比如把 4 和 1 合并為水域濕地把 7 單獨(dú)拎出來看鹽堿地不需要重新生成柵格效率最高。統(tǒng)計(jì)完可以和已知面積做一次粗校驗(yàn)吉林省陸地面積約 18.74 萬平方公里把上面腳本輸出的所有類面積加起來如果總誤差超過 2%說明裁剪或投影環(huán)節(jié)出了問題優(yōu)先檢查投影后的像元分辨率是否真的是 10 米。5. 吉林省應(yīng)用避坑五個真實(shí)踩坑記錄5.1 坑一WGS84 原始圖直接算面積水體面積憑空多了幾十萬畝現(xiàn)象拿到下載好的 GeoTIFF沒做投影就直接在 ArcGIS 里用屬性表統(tǒng)計(jì)各類像元數(shù)再按“10 米乘 10 米等于 100 平方米”換算面積結(jié)果查干湖水體面積比實(shí)際大出近三成。原因WGS84 地理坐標(biāo)系下像元尺寸用度表示緯度越高經(jīng)度方向的實(shí)際距離越短北緯 45 度附近一個 10 米標(biāo)稱的像元實(shí)際東西向長度只有赤道處的七成左右像元根本不是正方形。解決必須先轉(zhuǎn) Albers 等積投影讓像元變?yōu)?10 米乘 10 米的標(biāo)準(zhǔn)方格再統(tǒng)計(jì)像元數(shù)。5.2 坑二水稻田被分到 Flooded Vegetation作物面積嚴(yán)重偏低現(xiàn)象統(tǒng)計(jì)完 Crops 類別面積只有 400 多萬公頃和統(tǒng)計(jì)公報(bào)里的全省糧食播種面積對不上翻看影像發(fā)現(xiàn)大片稻田被分到了第 4 類“被淹植被”。原因吉林省水稻在移栽和返青期田間保持水層Sentinel-2 影像上反映的是“水 植被”混合光譜恰好命中 Flooded Vegetation 的訓(xùn)練特征。解決不要急著改柵格把第 4 類的空間分布疊加到 2021 年 5 月到 6 月的 Sentinel-2 NDVI 時序圖上凡是水層期出現(xiàn)又在水稻成熟期轉(zhuǎn)為高 NDVI 的像元手動歸并到 Crops 類統(tǒng)計(jì)口徑。5.3 坑三長白山冬季積雪造成樹木面積高估現(xiàn)象統(tǒng)計(jì)結(jié)果里 Trees 類面積達(dá)到 800 萬公頃比吉林省林業(yè)數(shù)據(jù)偏大而且 Snow/Ice 類也有相當(dāng)數(shù)量分布在林區(qū)。原因年度合成影像里如果有冬季月份被保留下來積雪覆蓋在林冠上分類模型把“雪 樹冠”的混合信號分成了冰雪或樹木造成雙向誤差。解決不要在年度合成圖上做逐類絕對面積先按月層瀏覽 2021 年 12 月、1 月、2 月的合成影像把雪蓋嚴(yán)重的月份找出來用夏季月份的合成層重新做年度統(tǒng)計(jì)或者直接用 NDVI 峰值合成剔除雪的干擾。5.4 坑四白城鹽堿地 Bare Ground 與旱地邊界非常碎現(xiàn)象西部鹽堿地區(qū)域的裸地類別呈椒鹽狀和周邊旱作農(nóng)田犬牙交錯統(tǒng)計(jì)面積時同一塊地每年結(jié)果波動極大。原因鹽堿地反射率隨土壤含水量和鹽分變化明顯部分鹽堿地的光譜和裸土幾乎一樣單時相特征難以區(qū)分。解決統(tǒng)計(jì)時先做一次 3×3 像元眾數(shù)濾波把孤立單像元的類別噪聲壓掉?；蛘甙衙娣e小于 0.1 公頃的零碎裸地圖斑直接合并到相鄰主導(dǎo)類別再重新統(tǒng)計(jì)。5.5 坑五分帶投影拼接造成邊界錯位現(xiàn)象整省按三度帶分四個帶分別投影后再鑲嵌結(jié)果帶與帶交界處的線性地物出現(xiàn)百米級錯位面積統(tǒng)計(jì)在接邊帶異常偏大。原因每個帶單獨(dú)做投影變換時遠(yuǎn)離中央經(jīng)線的區(qū)域變形方向不同接邊處的像元重采樣又用了各自帶內(nèi)的參數(shù)拼接后位置對不齊。解決整省范圍必須一次性投影到同一個 Albers 坐標(biāo)系不要分帶裁剪、分帶投影再拼接。如果柵格太大寧可分塊投影后用 gdalmerge 合并也不要按投影帶切分。6. 怎么驗(yàn)證這套 10 米結(jié)果是可信的與 GlobeLand30 做交叉驗(yàn)證驗(yàn)證方法比想象中簡單不需要地面樣點(diǎn)只需要另一套公開數(shù)據(jù)做參照。用 QGIS 或 ArcGIS 在吉林省范圍內(nèi)生成 300 到 500 個隨機(jī)點(diǎn)分別提取 ESRI 2021 10 米數(shù)據(jù)和 GlobeLand30 數(shù)據(jù)的類別值做一個混淆矩陣。重點(diǎn)關(guān)注 Trees、Crops、Bare Ground 這三類的互相混淆程度。總體精度在 85% 以上、Kappa 系數(shù)在 0.75 以上基本可以認(rèn)為這套 10 米數(shù)據(jù)在吉林省是可用的。如果 Trees 和 Crops 的混淆比例超過 15%多半是 2021 年影像質(zhì)量問題需要回到月度合成層看看哪些月份被云污染。面積驗(yàn)證用統(tǒng)計(jì)公報(bào)口徑把腳本統(tǒng)計(jì)出的 Crops 類面積換算成公頃與省統(tǒng)計(jì)年鑒里的耕地和播種面積對比。不同定義下差 10% 以內(nèi)都算正常因?yàn)?ESRI 的 Crops 類不包含未耕種但具有耕地屬性的土地。如果差幅超過 25%優(yōu)先懷疑第 5 章里的水稻田誤分問題其次檢查投影是否真的轉(zhuǎn)成了等積。還有一個更輕量的驗(yàn)證技巧選一個自己最熟悉的縣比如榆樹市或農(nóng)安縣把 10 米分類結(jié)果按 5 公里網(wǎng)格打上格網(wǎng)每個格網(wǎng)里計(jì)算 Crops 類占比再把占比圖和該縣統(tǒng)計(jì)年鑒里的鄉(xiāng)鎮(zhèn)播種面積表做排序相關(guān)性分析。柵格數(shù)據(jù)和統(tǒng)計(jì)數(shù)據(jù)能對上序就說明這套數(shù)據(jù)在空間格局上是可信的。從那以后我每次拿到任何年度土地覆蓋柵格第一件事不是打開看顏色而是先跑一遍像元數(shù)校驗(yàn)和與統(tǒng)計(jì)公報(bào)的面積對比三行腳本就能攔住大多數(shù)翻車。這套 2021 年吉林省 10 米數(shù)據(jù)的裁剪、投影、統(tǒng)計(jì)流程跑通之后你對“10 米到底比 30 米好在哪”會有非常具體的體感希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取 ESRI 2021土地利用數(shù)據(jù)吉林省10米分辨率裁剪投影與面積統(tǒng)計(jì)避坑指南
簡介2021年土地利用ESRI吉林省10m精度數(shù)據(jù)集由全球GIS領(lǐng)域頭部廠商ESRI官方制作采用WGS84坐標(biāo)系統(tǒng)與10米柵格分辨率完整呈現(xiàn)吉林省土地覆蓋狀況。數(shù)據(jù)源自海量全球原始柵格經(jīng)地理配準(zhǔn)、裁剪、轉(zhuǎn)換等處理按長春、吉林、四平、遼源、通化、白山、松原、白城、延邊9個地級行政區(qū)劃分獨(dú)立目錄壓縮包共63個文件大小約44MB核心為TIF格式土地覆蓋柵格圖配套TFW坐標(biāo)定位文件、DBF與CPG屬性表、XLSX分類統(tǒng)計(jì)表、XML元數(shù)據(jù)及PNG快速預(yù)覽圖。借助主流GIS軟件可提取農(nóng)田、森林、水域、城市建筑、草地等用地類型完成面積統(tǒng)計(jì)、變化檢測與空間制圖。數(shù)據(jù)適用于土地覆蓋變化監(jiān)測、生態(tài)環(huán)境評估、城鄉(xiāng)規(guī)劃、農(nóng)業(yè)管理、自然災(zāi)害風(fēng)險評估及氣候變化研究等多類任務(wù)。該數(shù)據(jù)集已有258人學(xué)習(xí)下載是GIS從業(yè)人員、規(guī)劃決策者與高校師生開展省級尺度空間分析的實(shí)用基礎(chǔ)數(shù)據(jù)。1. 2021 年 ESRI 土地利用數(shù)據(jù)吉林省 10 米分辨率到底能用來干什么很多做東北項(xiàng)目的同事電腦里還躺著 GlobeLand30 那套 30 米地表覆蓋理由就一句話“省級分析夠用了”。但 2021 年 ESRI 發(fā)布的這套全球 10 米土地利用數(shù)據(jù)把空間分辨率抬到了 Sentinel-2 能做到的極限吉林省作為長白山區(qū)和松嫩平原的交錯帶最吃這套數(shù)據(jù)的紅利山區(qū)能撕開溝谷和耕地邊界平原能大致分清村落、農(nóng)田和鹽堿地。不過它原始坐標(biāo)系是 WGS84 地理坐標(biāo)直接拿原圖算面積會翻車而且分類體系只有 9 類跟國內(nèi)常用的土地利用二級類對不上。這篇文章就圍繞“ESRI 2021 土地利用、吉林省、10 米精度”講清楚它是什么、怎么裁剪投影、怎么統(tǒng)計(jì)面積以及五個實(shí)際踩過的坑。適合 GIS 從業(yè)者、做國土空間規(guī)劃和農(nóng)業(yè)遙感監(jiān)測的人直接照著做。2. 分類體系與數(shù)據(jù)原理9 類不夠用但每一類都有明確含義2.1 九個類別字段每類對應(yīng)什么地表特征這套數(shù)據(jù)在 ArcGIS Living Atlas 里的全稱是 Sentinel-2 10m Land Use/Land Cover2021 年度合成版本使用 9 類體系。先看編碼后面所有面積統(tǒng)計(jì)都依賴這張表編碼英文類名中文對應(yīng)在吉林省的典型表現(xiàn)1Water水體松花江、查干湖、大型水庫2Trees樹木長白山針闊混交林、次生林3Grass草地西部草場、林間隙地4Flooded Vegetation被淹植被河漫灘、濕地、洪泛區(qū)植被5Crops作物玉米、水稻、大豆田塊6Built Area建成區(qū)長春、吉林市區(qū)及鄉(xiāng)鎮(zhèn)聚落7Bare Ground裸地白城鹽堿地、裸土、采石場8Snow/Ice冰雪冬季積雪也常造成誤分9Clouds云無意義類別統(tǒng)計(jì)前應(yīng)剔除注意第 4 類“被淹植被”在官方文檔里對應(yīng)的是紅樹林、洪泛區(qū)、沼澤這類長期或季節(jié)性淹水的植被并不是水稻田。但吉林省水稻移栽期田間有明水光譜特征和 Flooded Vegetation 很像后面避坑章會專門講這個事。2.2 物候合成與深度學(xué)習(xí)的生產(chǎn)邏輯ESRI 這套產(chǎn)品不是拿單景影像做的分類而是用 2021 年全年 Sentinel-2 L2A 數(shù)據(jù)按月合成每個像元保留最能代表該月地表狀態(tài)的反射率值再把 12 個月的合成影像輸入到一個深度學(xué)習(xí)語義分割模型里逐像元輸出類別標(biāo)簽。這樣做的直接結(jié)果是常綠針葉林和落葉闊葉林在年度合成里會被歸到同一個“Trees”類這在東北是能接受的畢竟林業(yè)二類調(diào)查用的是亞類而土地利用現(xiàn)狀調(diào)查要的是喬木林地這個一級類口徑。和 30 米 GlobeLand30 相比10 米產(chǎn)品在吉林省的實(shí)際差異非常直觀。第一是線性地物鄉(xiāng)間道路、河渠、防護(hù)林帶在 30 米像元里基本糊成混合像元10 米能單獨(dú)成行。第二是聚落邊界長春周邊城中村、開發(fā)區(qū)在建工地10 米分辨率下建成區(qū)和裸地的邊界干凈很多。第三是耕地破碎度吉林東部半山區(qū)的小塊農(nóng)田30 米產(chǎn)品經(jīng)常漏分10 米分出來的地塊更接近真實(shí)田塊邊界。2.3 為什么說它不能直接對標(biāo)國土三調(diào)國內(nèi)用戶最常問的問題這套數(shù)據(jù)和國土變更調(diào)查里的地類圖斑能不能直接互相替換答案是不能。三調(diào)的“水田”要求結(jié)合水利設(shè)施和種植制度認(rèn)定ESRI 的 Crops 類只反映 2021 年影像上“像作物的植被”兩者口徑不同。這套數(shù)據(jù)真正適合的場景是宏觀的土地利用變化監(jiān)測、生態(tài)評估、國土空間規(guī)劃的前期現(xiàn)狀摸底而不是地籍級的確權(quán)應(yīng)用。把它定位成“省級宏觀現(xiàn)狀速覽圖”最合適精度優(yōu)勢體現(xiàn)在空間表達(dá)上不在分類語義深度上。3. 從 Living Atlas 到吉林省 GeoTIFF三種獲取路徑與導(dǎo)出參數(shù)3.1 ArcGIS Pro 動態(tài)影像服務(wù)導(dǎo)出最穩(wěn)妥的獲取方式是直接在 ArcGIS Pro 里連 Living Atlas 的影像服務(wù)。打開目錄窗格在 Portal 選項(xiàng)卡下瀏覽 Living Atlas搜索 Sentinel-2 10m Land Use/Land Cover把它拖到地圖視圖。注意要選 2021 年份的圖層這個產(chǎn)品按年份拆成了多條服務(wù)混用年份會導(dǎo)致后續(xù)分析出錯。確認(rèn)圖層加載后通過“數(shù)據(jù)管理 → 導(dǎo)出柵格”或在圖層右鍵菜單里選“導(dǎo)出 → 為不同格式”。關(guān)鍵參數(shù)如下Extent范圍選擇“繪制范圍”并加載吉林省邊界矢量或用“環(huán)境 → 處理范圍”指定邊界要素類。Cell Size設(shè)為 10 米。Resampling必須選 Nearest Neighbor這是分類柵格雙線性插值會在類別邊緣產(chǎn)生不存在的地類編碼。FormatTIFF。CompressionLZW 或 DEFLATE10 米整省數(shù)據(jù)壓縮后體積能小很多。整省范圍直接在動態(tài)服務(wù)上導(dǎo)出經(jīng)常遇到服務(wù)端超時。我一般會先把吉林省邊界外擴(kuò) 1 公里作為導(dǎo)出范圍然后把目標(biāo)拆成東部山區(qū)、中部平原、西部草原三塊分別導(dǎo)出最后再用鑲嵌工具合并。動態(tài)服務(wù)的好處是不用關(guān)心原始分塊和你有沒有 1T 硬盤壞處是受服務(wù)性能限制大范圍導(dǎo)出容易斷。3.2 官方分塊 GeoTIFF 與 QGIS 加載ESRI 在 Living Atlas 的數(shù)據(jù)詳情頁還提供了按 100km × 100km 分塊的年度 GeoTIFF 直接下載入口每個分塊文件名包含它的經(jīng)緯度網(wǎng)格編號。吉林省橫跨大約 122°E 到 131°E、41°N 到 46°N需要匹配 4 到 6 個分塊文件下載后先鑲嵌成一張全圖再做裁剪。這個路徑適合網(wǎng)絡(luò)條件好、想離線分析的用戶。如果你用 QGIS動態(tài)影像服務(wù)同樣可以加載圖層 → 添加圖層 → ArcGIS REST Server粘貼服務(wù)的 REST URL。但沒有 ArcGIS Pro 的“導(dǎo)出柵格”那么順手需要先加載到畫布再用“柵格 → 轉(zhuǎn)換 → 裁剪”按邊界矢量導(dǎo)出速度比 Pro 慢不少內(nèi)存也吃得更緊。3.3 導(dǎo)出后先做三個檢查下載完成先別急著投影和統(tǒng)計(jì)花兩分鐘做三個檢查能省掉后面很多排查時間。第一檢查 NoData 設(shè)置。動態(tài)服務(wù)導(dǎo)出的 GeoTIFF 在省界外通常是有值的只是顯示為黑色背景要確認(rèn) NoData 有沒有正確寫入。第二檢查類別直方圖用 ArcGIS 的“查看屬性表”或 QGIS 的“唯一值”統(tǒng)計(jì)看看 9Clouds類的像元數(shù)量是否異常多如果超過總像元數(shù)的 5%說明 2021 年該區(qū)域云影響嚴(yán)重需要重新評估這個年份的質(zhì)量。第三檢查空間參考確認(rèn)文件仍是 WGS84 地理坐標(biāo)系如果已經(jīng)被自動轉(zhuǎn)成了 Web Mercator先退回原坐標(biāo)系再走下一步。4. 裁剪、投影與面積統(tǒng)計(jì)可抄作業(yè)的 GDAL 參數(shù)4.1 用邊界矢量做精確裁剪拿到吉林省范圍的 GeoTIFF 后第一步是裁剪。推薦用 GDAL 的 gdalwarp 一次性完成裁剪加范圍收緊命令行如下gdalwarp -overwrite \ -cutline jilin_2021.shp \ -crop_to_cutline \ -dstnodata 255 \ -r near \ -of GTiff \ LULC_2021_jilin_raw.tif \ LULC_2021_jilin_wgs84.tif-cutline指定吉林省邊界矢量-crop_to_cutline讓輸出柵格的范圍嚴(yán)格貼合邊界省界外不再保留多余的矩形區(qū)域。-dstnodata 255把邊界外的像元設(shè)為 255這樣后面統(tǒng)計(jì)時可以直接排除。-r near保持最近鄰重采樣分類柵格不能做平滑插值。如果你手里的原始數(shù)據(jù)是多個分塊先執(zhí)行g(shù)dalbuildvrt LULC_2021_jilin.vrt block1.tif block2.tif block3.tif block4.tif gdal_translate -of GTiff LULC_2021_jilin.vrt LULC_2021_jilin_raw.tifVRT 虛擬柵格的好處是零拷貝合并不產(chǎn)生中間文件適合分塊較多的情況。合并完成后還有一步很關(guān)鍵用 gdalinfo 查看輸出文件的尺寸確認(rèn)沒有因?yàn)榉謮K重疊產(chǎn)生坐標(biāo)偏移。4.2 為什么用 Albers 等積投影而不是三度帶吉林省跨越東經(jīng) 121 度到 131 度橫跨多個三度帶。很多習(xí)慣用國家 2000 三度帶投影的同行會按帶號分帶投影最后再拼接。這個做法在縣一級沒問題但整省范圍會出現(xiàn)兩個麻煩一是分帶接邊處地物錯位二是兩個帶分別做面積統(tǒng)計(jì)后再合并類別面積在拼接帶附近會有系統(tǒng)性偏差。做土地利用面積統(tǒng)計(jì)第一原則是投影必須等積。我一般用 Albers 等面積圓錐投影參數(shù)按東北地區(qū)常用設(shè)置gdalwarp -overwrite \ -t_srs projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 ellpsGRS80 unitsm no_defs \ -r near \ -tr 10 10 \ LULC_2021_jilin_wgs84.tif \ LULC_2021_jilin_albers.tif投影參數(shù)里lat_125和lat_247是雙標(biāo)準(zhǔn)緯線吉林省主體緯度 41 度到 46 度正好被夾在中間變形最小。lon_0105是中央經(jīng)線這是全國標(biāo)準(zhǔn) Albers 的常用值如果你只做吉林省改成 126 度會讓圖面中央變形更小但面積統(tǒng)計(jì)結(jié)果和 105 度幾乎沒有差別。-tr 10 10強(qiáng)制輸出像元為 10 米乘 10 米的正方形這一步直接決定后續(xù)面積統(tǒng)計(jì)的精度——投影后每個像元對應(yīng) 100 平方米。4.3 分塊統(tǒng)計(jì)各類面積避免一次性讀入內(nèi)存面積統(tǒng)計(jì)的思路很簡單統(tǒng)計(jì)每個類別編碼的像元數(shù)乘以單個像元面積。但吉林省全境 10 米分辨率柵格寬度約 74000 個像元高度約 62000 個像元總像元數(shù)超過 45 億。直接用 GDAL 的 ReadAsArray 把全圖讀進(jìn)內(nèi)存uint8 數(shù)據(jù)也要超過 4.5GB普通工作站會直接內(nèi)存溢出。正確做法是分塊讀取from osgeo import gdal from collections import Counter import numpy as np ds gdal.Open(LULC_2021_jilin_albers.tif, 0) band ds.GetRasterBand(1) xsize band.XSize ysize band.YSize # 讀取仿射變換得到單像元的寬度和高度乘起來就是單像元面積 gt ds.GetGeoTransform() px_area abs(gt[1] * gt[5]) print(f單像元面積: {px_area:.1f} 平方米) # 分塊大小取 2048 行寬度保持全幅單次內(nèi)存約 74000*2048*1 字節(jié) ≈ 150MB block_size 2048 counter Counter() for row in range(0, ysize, block_size): rows min(block_size, ysize - row) data band.ReadAsArray(0, row, xsize, rows) vals, counts np.unique(data, return_countsTrue) for v, c in zip(vals, counts): if v not in (255, 9): # 剔除邊界 NoData 和云類 counter[v] c print(f已處理第 {row}/{ysize} 行) # 輸出各面積平方米轉(zhuǎn)公頃 class_names {1: 水體, 2: 樹木, 3: 草地, 4: 被淹植被, 5: 作物, 6: 建成區(qū), 7: 裸地, 8: 冰雪} for cls_id in sorted(class_names): ha counter[cls_id] * px_area / 10000 print(f{class_names[cls_id]}: {ha:.1f} 公頃)px_area的值在 Albers 投影下就是 100.0 平方米它的正確性是面積統(tǒng)計(jì)可信的基礎(chǔ)。block_size可以按機(jī)器內(nèi)存調(diào)整內(nèi)存 16GB 以上的機(jī)器可以設(shè) 40968GB 的機(jī)器建議保持 2048。輸出結(jié)果后再在 Excel 里手動合并類別比如把 4 和 1 合并為水域濕地把 7 單獨(dú)拎出來看鹽堿地不需要重新生成柵格效率最高。統(tǒng)計(jì)完可以和已知面積做一次粗校驗(yàn)吉林省陸地面積約 18.74 萬平方公里把上面腳本輸出的所有類面積加起來如果總誤差超過 2%說明裁剪或投影環(huán)節(jié)出了問題優(yōu)先檢查投影后的像元分辨率是否真的是 10 米。5. 吉林省應(yīng)用避坑五個真實(shí)踩坑記錄5.1 坑一WGS84 原始圖直接算面積水體面積憑空多了幾十萬畝現(xiàn)象拿到下載好的 GeoTIFF沒做投影就直接在 ArcGIS 里用屬性表統(tǒng)計(jì)各類像元數(shù)再按“10 米乘 10 米等于 100 平方米”換算面積結(jié)果查干湖水體面積比實(shí)際大出近三成。原因WGS84 地理坐標(biāo)系下像元尺寸用度表示緯度越高經(jīng)度方向的實(shí)際距離越短北緯 45 度附近一個 10 米標(biāo)稱的像元實(shí)際東西向長度只有赤道處的七成左右像元根本不是正方形。解決必須先轉(zhuǎn) Albers 等積投影讓像元變?yōu)?10 米乘 10 米的標(biāo)準(zhǔn)方格再統(tǒng)計(jì)像元數(shù)。5.2 坑二水稻田被分到 Flooded Vegetation作物面積嚴(yán)重偏低現(xiàn)象統(tǒng)計(jì)完 Crops 類別面積只有 400 多萬公頃和統(tǒng)計(jì)公報(bào)里的全省糧食播種面積對不上翻看影像發(fā)現(xiàn)大片稻田被分到了第 4 類“被淹植被”。原因吉林省水稻在移栽和返青期田間保持水層Sentinel-2 影像上反映的是“水 植被”混合光譜恰好命中 Flooded Vegetation 的訓(xùn)練特征。解決不要急著改柵格把第 4 類的空間分布疊加到 2021 年 5 月到 6 月的 Sentinel-2 NDVI 時序圖上凡是水層期出現(xiàn)又在水稻成熟期轉(zhuǎn)為高 NDVI 的像元手動歸并到 Crops 類統(tǒng)計(jì)口徑。5.3 坑三長白山冬季積雪造成樹木面積高估現(xiàn)象統(tǒng)計(jì)結(jié)果里 Trees 類面積達(dá)到 800 萬公頃比吉林省林業(yè)數(shù)據(jù)偏大而且 Snow/Ice 類也有相當(dāng)數(shù)量分布在林區(qū)。原因年度合成影像里如果有冬季月份被保留下來積雪覆蓋在林冠上分類模型把“雪 樹冠”的混合信號分成了冰雪或樹木造成雙向誤差。解決不要在年度合成圖上做逐類絕對面積先按月層瀏覽 2021 年 12 月、1 月、2 月的合成影像把雪蓋嚴(yán)重的月份找出來用夏季月份的合成層重新做年度統(tǒng)計(jì)或者直接用 NDVI 峰值合成剔除雪的干擾。5.4 坑四白城鹽堿地 Bare Ground 與旱地邊界非常碎現(xiàn)象西部鹽堿地區(qū)域的裸地類別呈椒鹽狀和周邊旱作農(nóng)田犬牙交錯統(tǒng)計(jì)面積時同一塊地每年結(jié)果波動極大。原因鹽堿地反射率隨土壤含水量和鹽分變化明顯部分鹽堿地的光譜和裸土幾乎一樣單時相特征難以區(qū)分。解決統(tǒng)計(jì)時先做一次 3×3 像元眾數(shù)濾波把孤立單像元的類別噪聲壓掉?；蛘甙衙娣e小于 0.1 公頃的零碎裸地圖斑直接合并到相鄰主導(dǎo)類別再重新統(tǒng)計(jì)。5.5 坑五分帶投影拼接造成邊界錯位現(xiàn)象整省按三度帶分四個帶分別投影后再鑲嵌結(jié)果帶與帶交界處的線性地物出現(xiàn)百米級錯位面積統(tǒng)計(jì)在接邊帶異常偏大。原因每個帶單獨(dú)做投影變換時遠(yuǎn)離中央經(jīng)線的區(qū)域變形方向不同接邊處的像元重采樣又用了各自帶內(nèi)的參數(shù)拼接后位置對不齊。解決整省范圍必須一次性投影到同一個 Albers 坐標(biāo)系不要分帶裁剪、分帶投影再拼接。如果柵格太大寧可分塊投影后用 gdalmerge 合并也不要按投影帶切分。6. 怎么驗(yàn)證這套 10 米結(jié)果是可信的與 GlobeLand30 做交叉驗(yàn)證驗(yàn)證方法比想象中簡單不需要地面樣點(diǎn)只需要另一套公開數(shù)據(jù)做參照。用 QGIS 或 ArcGIS 在吉林省范圍內(nèi)生成 300 到 500 個隨機(jī)點(diǎn)分別提取 ESRI 2021 10 米數(shù)據(jù)和 GlobeLand30 數(shù)據(jù)的類別值做一個混淆矩陣。重點(diǎn)關(guān)注 Trees、Crops、Bare Ground 這三類的互相混淆程度。總體精度在 85% 以上、Kappa 系數(shù)在 0.75 以上基本可以認(rèn)為這套 10 米數(shù)據(jù)在吉林省是可用的。如果 Trees 和 Crops 的混淆比例超過 15%多半是 2021 年影像質(zhì)量問題需要回到月度合成層看看哪些月份被云污染。面積驗(yàn)證用統(tǒng)計(jì)公報(bào)口徑把腳本統(tǒng)計(jì)出的 Crops 類面積換算成公頃與省統(tǒng)計(jì)年鑒里的耕地和播種面積對比。不同定義下差 10% 以內(nèi)都算正常因?yàn)?ESRI 的 Crops 類不包含未耕種但具有耕地屬性的土地。如果差幅超過 25%優(yōu)先懷疑第 5 章里的水稻田誤分問題其次檢查投影是否真的轉(zhuǎn)成了等積。還有一個更輕量的驗(yàn)證技巧選一個自己最熟悉的縣比如榆樹市或農(nóng)安縣把 10 米分類結(jié)果按 5 公里網(wǎng)格打上格網(wǎng)每個格網(wǎng)里計(jì)算 Crops 類占比再把占比圖和該縣統(tǒng)計(jì)年鑒里的鄉(xiāng)鎮(zhèn)播種面積表做排序相關(guān)性分析。柵格數(shù)據(jù)和統(tǒng)計(jì)數(shù)據(jù)能對上序就說明這套數(shù)據(jù)在空間格局上是可信的。從那以后我每次拿到任何年度土地覆蓋柵格第一件事不是打開看顏色而是先跑一遍像元數(shù)校驗(yàn)和與統(tǒng)計(jì)公報(bào)的面積對比三行腳本就能攔住大多數(shù)翻車。這套 2021 年吉林省 10 米數(shù)據(jù)的裁剪、投影、統(tǒng)計(jì)流程跑通之后你對“10 米到底比 30 米好在哪”會有非常具體的體感希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取