簡介這份資源為2020年廣東省10米精度土地覆蓋土地利用數(shù)據(jù)面向地理信息、遙感、城市規(guī)劃及生態(tài)環(huán)境研究等方向的學(xué)習(xí)者與從業(yè)者。數(shù)據(jù)基于哨兵影像與深度學(xué)習(xí)方法制作經(jīng)墨卡托投影轉(zhuǎn)WGS84地理坐標(biāo)系并按最新省市級(jí)行政邊界裁剪覆蓋廣東省各地級(jí)市共分耕地、林地、草地、灌木、濕地、水體、不透水面、裸地、雪冰等十類。壓縮包共147個(gè)文件約75.8MB包含21個(gè)tif柵格主數(shù)據(jù)、21個(gè)xlsx屬性表、21個(gè)dbf與cpg等配套文件、21個(gè)tfw坐標(biāo)信息、21個(gè)xml元數(shù)據(jù)及21個(gè)png預(yù)覽圖便于直接加載與屬性查詢。已有329人學(xué)習(xí)下載。讀者可獲得統(tǒng)一坐標(biāo)系、邊界規(guī)范的廣東省市級(jí)土地利用數(shù)據(jù)用于空間分析、制圖與變化監(jiān)測省去投影轉(zhuǎn)換與裁剪步驟。1. 2020年10m精度廣東省土地覆蓋土地利用數(shù)據(jù)這份柵格包到底能干什么如果你手頭正好缺一份能直接扔進(jìn) GIS 里做分析的廣東省地表覆蓋底圖又不想花幾天時(shí)間去拼 Landsat 或者等哨兵影像慢慢下載那這份「2020年10m精度廣東省土地覆蓋土地利用.rar」大概率就是你要找的東西。它本質(zhì)上是一套已經(jīng)分類好的柵格數(shù)據(jù)產(chǎn)品空間分辨率 10 米覆蓋廣東全省時(shí)間截面鎖定在 2020 年。10 米這個(gè)精度意味著什么一個(gè)像素代表地面上 10 米 × 10 米的方塊廣州珠江新城一個(gè)標(biāo)準(zhǔn)地鐵站出入口的面積大概也就幾十個(gè)像素用來做市域、縣域尺度的國土空間分析、生態(tài)評(píng)估、城市擴(kuò)張監(jiān)測完全夠用。它適合做自然資源調(diào)查的技術(shù)人員、做城市研究的規(guī)劃師、寫論文需要土地利用底圖的研究生以及做遙感應(yīng)用開發(fā)需要一份現(xiàn)成分類結(jié)果做驗(yàn)證的工程師。但先別急著解壓這份數(shù)據(jù)能不能用、怎么用、坑在哪得先把它拆開看明白。2. 10米分辨率土地覆蓋柵格從分類體系到廣東地類的對(duì)應(yīng)關(guān)系2.1 為什么是10米而不是30米或1米土地覆蓋產(chǎn)品的分辨率選擇從來不是拍腦袋定的。30 米級(jí)別的產(chǎn)品比如早期基于 Landsat 的全球地表覆蓋在省級(jí)尺度上做宏觀趨勢沒問題但一旦落到珠三角這種城市蔓延劇烈、地塊破碎度高的區(qū)域30 米像素會(huì)把一個(gè)城中村和旁邊的工業(yè)園混成一個(gè)類別邊界模糊得讓人抓狂。1 米級(jí)別的高分影像當(dāng)然精細(xì)但廣東省全域覆蓋的高分?jǐn)?shù)據(jù)獲取成本高、處理量大而且分類精度受陰影和建筑紋理干擾反而可能下降。10 米恰好卡在一個(gè)甜點(diǎn)上哨兵二號(hào)Sentinel-2的多光譜波段原生就是 10 米、20 米和 60 米其中可見光和近紅外是 10 米用這些波段做分類既保留了城市內(nèi)部地塊的區(qū)分度又不至于讓數(shù)據(jù)量爆炸。這份 2020 年的產(chǎn)品大概率就是基于哨兵二號(hào)時(shí)序影像生產(chǎn)的這也是目前省級(jí) 10 米土地覆蓋產(chǎn)品的主流技術(shù)路線。2.2 分類體系怎么讀一級(jí)類和二級(jí)類的映射拿到柵格后第一個(gè)要搞清楚的是它的分類編碼。國內(nèi)常見的土地覆蓋產(chǎn)品通常參考《土地利用現(xiàn)狀分類》國標(biāo)或者類似 FROM-GLC、CLCD 等產(chǎn)品的分類體系。這份數(shù)據(jù)大概率包含耕地、林地、草地、灌木地、濕地、水體、不透水面、裸地這幾大類。你需要找到隨數(shù)據(jù)附帶的說明文件或者屬性表確認(rèn)每個(gè)像元值對(duì)應(yīng)什么類別。常見做法是值 1 代表耕地2 代表林地3 代表草地4 代表灌木5 代表濕地6 代表水體7 代表不透水面8 代表裸地。但不同生產(chǎn)單位編碼順序可能不同所以務(wù)必以實(shí)際說明為準(zhǔn)。如果壓縮包里沒有說明文件那就得靠目視比對(duì)打開珠三角區(qū)域看哪個(gè)值大面積出現(xiàn)在城市核心區(qū)那個(gè)值大概率就是不透水面。2.3 數(shù)據(jù)格式與坐標(biāo)系統(tǒng)確認(rèn)解壓后你大概率會(huì)看到一個(gè)或多個(gè) .tif 文件可能按地級(jí)市分幅也可能是一整幅廣東全省的鑲嵌影像。坐標(biāo)系統(tǒng)通常是 WGS84 地理坐標(biāo)系EPSG:4326或者 CGCS2000 投影坐標(biāo)系。如果是 4326單位是度做面積統(tǒng)計(jì)時(shí)需要先投影到合適的投影坐標(biāo)系比如 CGCS2000 3 度帶高斯克呂格投影否則算出來的面積是平方度毫無意義。檢查方法很簡單在 QGIS 或 ArcGIS 里加載后看圖層屬性里的 CRS 信息。如果是分幅的還需要先鑲嵌成一整幅再裁剪到研究區(qū)邊界。# 用 gdalinfo 查看柵格基本信息確認(rèn)坐標(biāo)系、行列數(shù)、像元大小和 NoData 值 gdalinfo gd_landcover_2020.tif # 關(guān)鍵輸出解讀 # Size is 120000, 85000 - 全省幅面行列數(shù)很大 # Pixel Size (0.0000898, -0.0000898) - 約10米分辨率地理坐標(biāo)系 # Coordinate System is WGS 84 - EPSG:4326 # NoData Value 0 - 背景值統(tǒng)計(jì)時(shí)要排除上面這段gdalinfo是拿到任何柵格數(shù)據(jù)后的第一步。重點(diǎn)看四個(gè)東西行列數(shù)判斷數(shù)據(jù)范圍是否符合預(yù)期像元大小確認(rèn)是不是真的 10 米坐標(biāo)系決定后續(xù)要不要投影轉(zhuǎn)換NoData 值決定統(tǒng)計(jì)時(shí)怎么設(shè)掩膜。如果 NoData 是 0 而你的分類編碼里 0 恰好代表某種地類那就得小心了需要跟數(shù)據(jù)說明核對(duì)清楚。2.4 用 Python 快速做面積統(tǒng)計(jì)確認(rèn)完基本信息后最常做的操作就是統(tǒng)計(jì)各地類面積。下面這段代碼用 rasterio 和 numpy 實(shí)現(xiàn)邏輯清晰適合直接抄。import rasterio import numpy as np # 打開柵格文件 with rasterio.open(gd_landcover_2020.tif) as src: data src.read(1) # 讀第一個(gè)波段 transform src.transform nodata src.nodata # 計(jì)算單個(gè)像元面積平方米 # transform[0] 是 x 方向像元大小度transform[4] 是 y 方向像元大小 # 地理坐標(biāo)系下需要換算成米粗略按緯度30度附近1度約111km和96km pixel_width_m transform[0] * 111000 * np.cos(np.radians(23.5)) # 廣東平均緯度約23.5 pixel_height_m abs(transform[4]) * 111000 pixel_area_m2 pixel_width_m * pixel_height_m # 排除 NoData valid data ! nodata unique, counts np.unique(data[valid], return_countsTrue) # 輸出各地類面積平方公里 for cls, cnt in zip(unique, counts): area_km2 cnt * pixel_area_m2 / 1e6 print(f類別 {cls}: {area_km2:.2f} 平方公里)這段代碼的核心邏輯是先算出一個(gè)像元代表多少平方米再統(tǒng)計(jì)每個(gè)類別有多少個(gè)像元相乘得到總面積。參數(shù)方面transform[0]和transform[4]是仿射變換里的像元尺寸np.cos(np.radians(23.5))是因?yàn)榈乩碜鴺?biāo)系下經(jīng)度方向的實(shí)際距離隨緯度變化廣東跨緯度約 20.5 到 25.5 度取中間值 23.5 度做近似。如果你追求更精確的面積建議先用gdalwarp投影到等面積投影再做統(tǒng)計(jì)那樣每個(gè)像元面積恒定不用做余弦校正。3. 把柵格用起來裁剪、重分類與疊加分析的操作鏈3.1 按行政區(qū)裁剪用掩膜提取研究區(qū)全省數(shù)據(jù)直接分析往往沒必要你通常只關(guān)心某個(gè)市或某個(gè)流域。裁剪有兩種常見做法用矢量邊界做掩膜提取或者按坐標(biāo)范圍做窗口讀取。前者更精準(zhǔn)后者更快。下面用gdalwarp做矢量裁剪。# 用廣州市行政邊界裁剪土地覆蓋柵格 # -cutline 指定矢量邊界文件 # -crop_to_cutline 讓輸出范圍嚴(yán)格貼合邊界 # -dstnodata 設(shè)置輸出背景值 gdalwarp -cutline guangzhou_boundary.shp \ -crop_to_cutline \ -dstnodata 0 \ -tr 0.0000898 0.0000898 \ gd_landcover_2020.tif \ gz_landcover_2020.tif-cutline后面跟矢量文件路徑支持 shp、geojson 等格式。-crop_to_cutline是關(guān)鍵參數(shù)不加的話輸出還是全省范圍只是邊界外被設(shè)為 NoData文件大小沒變。-tr指定輸出分辨率這里保持和原數(shù)據(jù)一致避免重采樣引入誤差。裁剪完之后再用前面的 Python 腳本統(tǒng)計(jì)廣州各地類面積就方便多了。3.2 重分類把二級(jí)類合并成你需要的一級(jí)類原始數(shù)據(jù)的分類可能很細(xì)比如把耕地分成水田和旱地把林地分成有林地和灌木林。但你的分析可能只需要「耕地、林地、水體、建設(shè)用地」這幾大類。這時(shí)候就需要重分類。用gdal_calc或者 Python 的 numpy 都能做。import rasterio import numpy as np with rasterio.open(gz_landcover_2020.tif) as src: data src.read(1) profile src.profile # 假設(shè)原始編碼1水田 2旱地 3有林地 4灌木 5水體 6建設(shè)用地 7裸地 # 重分類規(guī)則1,2-1(耕地) 3,4-2(林地) 5-3(水體) 6-4(建設(shè)用地) 7-5(其他) reclass_map {1:1, 2:1, 3:2, 4:2, 5:3, 6:4, 7:5} reclassified np.zeros_like(data) for old, new in reclass_map.items(): reclassified[data old] new # 保持 NoData reclassified[data 0] 0 profile.update(dtyperasterio.uint8, nodata0) with rasterio.open(gz_landcover_reclass.tif, w, **profile) as dst: dst.write(reclassified, 1)重分類的邏輯就是建立一個(gè)舊值到新值的映射字典然后遍歷每個(gè)舊值把對(duì)應(yīng)像元賦成新值。注意np.zeros_like(data)創(chuàng)建的全零數(shù)組如果原始數(shù)據(jù)里有 0 值代表某種地類這里就會(huì)混淆所以務(wù)必確認(rèn) NoData 和有效類別不重疊。profile.update里把數(shù)據(jù)類型改成 uint8因?yàn)橹胤诸惡箢悇e數(shù)變少了用不著原來的 uint16 或 float能省不少存儲(chǔ)空間。3.3 疊加分析土地覆蓋與人口格網(wǎng)的交叉統(tǒng)計(jì)土地覆蓋數(shù)據(jù)很少單獨(dú)用通常要和其他空間數(shù)據(jù)疊加。比如你想知道廣州市每個(gè)公里格網(wǎng)里建設(shè)用地占比和人口密度的關(guān)系就需要把土地覆蓋柵格和人口格網(wǎng)做分區(qū)統(tǒng)計(jì)。常見做法是先用gdalwarp把兩套數(shù)據(jù)重采樣到同一分辨率和坐標(biāo)系然后用 Python 做逐像元計(jì)算或者用 QGIS 的 Zonal Statistics 插件。import rasterio import numpy as np import pandas as pd # 讀取土地覆蓋和人口柵格確保兩者已經(jīng)對(duì)齊 with rasterio.open(gz_landcover_reclass.tif) as src: landcover src.read(1) with rasterio.open(gz_population_2020.tif) as src: population src.read(1) # 創(chuàng)建 1km 格網(wǎng)編號(hào)假設(shè)分辨率約10米100個(gè)像元約1km block_size 100 rows, cols landcover.shape results [] for i in range(0, rows, block_size): for j in range(0, cols, block_size): lc_block landcover[i:iblock_size, j:jblock_size] pop_block population[i:iblock_size, j:jblock_size] if lc_block.size 0: continue # 計(jì)算建設(shè)用地類別4占比 built_ratio np.sum(lc_block 4) / lc_block.size # 計(jì)算平均人口 mean_pop np.mean(pop_block) results.append({built_ratio: built_ratio, mean_pop: mean_pop}) df pd.DataFrame(results) print(df.corr()) # 輸出相關(guān)系數(shù)矩陣這段代碼把影像切成 100×100 像元的塊每塊約 1 公里見方然后統(tǒng)計(jì)每塊里建設(shè)用地的比例和平均人口最后算相關(guān)系數(shù)。block_size可以根據(jù)你的研究尺度調(diào)整想粗一點(diǎn)就設(shè) 200 或 500。注意人口柵格如果是 WorldPop 或 GPW 這類產(chǎn)品單位可能是每像元人數(shù)直接取平均沒問題如果是密度值那就要乘以面積換算成人數(shù)再統(tǒng)計(jì)。3.4 變化檢測的前置準(zhǔn)備和 2010 年數(shù)據(jù)對(duì)齊如果你手頭還有 2010 年的同類型數(shù)據(jù)想做十年變化檢測那第一步是確保兩期數(shù)據(jù)坐標(biāo)系、分辨率、范圍完全一致。常見做法是以 2020 年數(shù)據(jù)為基準(zhǔn)把 2010 年數(shù)據(jù)重采樣和裁剪到相同網(wǎng)格。# 將2010年數(shù)據(jù)對(duì)齊到2020年的網(wǎng)格 gdalwarp -t_srs EPSG:4326 \ -tr 0.0000898 0.0000898 \ -te 109.5 20.5 117.5 25.5 \ -r near \ gd_landcover_2010.tif \ gd_landcover_2010_aligned.tif-t_srs指定目標(biāo)坐標(biāo)系-tr指定目標(biāo)分辨率-te指定目標(biāo)范圍xmin ymin xmax ymax-r near表示用最近鄰重采樣這對(duì)分類數(shù)據(jù)是必須的不能用雙線性或立方卷積否則會(huì)出現(xiàn)不存在的類別值。對(duì)齊之后兩期數(shù)據(jù)逐像元相減非零的地方就是變化區(qū)域。4. 避坑與排查10米土地覆蓋數(shù)據(jù)常見的五個(gè)翻車點(diǎn)4.1 面積算出來偏大或偏小現(xiàn)象用地理坐標(biāo)系直接統(tǒng)計(jì)面積發(fā)現(xiàn)各地類加總跟官方公布的行政區(qū)面積對(duì)不上偏差能到百分之十幾。原因地理坐標(biāo)系下像元面積隨緯度變化高緯度地區(qū)像元實(shí)際面積比低緯度小用固定值乘會(huì)引入系統(tǒng)誤差。解決先投影到等面積投影或高斯克呂格投影再做統(tǒng)計(jì)。廣東常用 CGCS2000 3 度帶中央經(jīng)線根據(jù)研究區(qū)經(jīng)度選擇比如 114°E 對(duì)應(yīng) EPSG:4547 左右。投影后像元面積恒定統(tǒng)計(jì)結(jié)果才可靠。4.2 分類編碼對(duì)不上說明書現(xiàn)象按說明書里的編碼去提取水體結(jié)果提取出來的是大片山區(qū)。原因不同生產(chǎn)批次或不同來源的數(shù)據(jù)編碼可能調(diào)整過說明書沒同步更新。解決不要盲信文檔先做目視驗(yàn)證。在 QGIS 里加載柵格疊加天地圖或影像底圖找?guī)讉€(gè)已知地物點(diǎn)比如廣州塔附近應(yīng)該是建設(shè)用地珠江應(yīng)該是水體查看像元值反推編碼含義。確認(rèn)后再批量處理。4.3 分幅數(shù)據(jù)接邊處出現(xiàn)裂縫現(xiàn)象把各地級(jí)市分幅數(shù)據(jù)鑲嵌后發(fā)現(xiàn)市界處有一條明顯的類別突變線。原因分幅生產(chǎn)時(shí)各幅獨(dú)立分類接邊處分類結(jié)果不一致或者鑲嵌時(shí)沒有做羽化處理。解決如果數(shù)據(jù)是分幅的優(yōu)先找全省整幅版本。如果沒有鑲嵌后用眾數(shù)濾波做一下后處理或者在接受范圍內(nèi)忽略接邊處的少量誤差。做變化檢測時(shí)尤其要注意接邊處的假變化會(huì)干擾結(jié)果。4.4 NoData 值被當(dāng)成有效類別統(tǒng)計(jì)現(xiàn)象統(tǒng)計(jì)結(jié)果里多出一個(gè)面積巨大的未知類別或者面積加總遠(yuǎn)超全省面積。原因NoData 值通常是 0 或 255沒有排除被當(dāng)成一個(gè)地類參與了統(tǒng)計(jì)。解決統(tǒng)計(jì)前先確認(rèn) NoData 值用data ! nodata做掩膜。如果 NoData 是 0 而 0 又恰好是某個(gè)地類的編碼那就得回去找生產(chǎn)方確認(rèn)或者用柵格邊緣的 0 值區(qū)域判斷——通常 NoData 會(huì)出現(xiàn)在影像邊界外。4.5 重采樣方法選錯(cuò)導(dǎo)致類別污染現(xiàn)象把 10 米數(shù)據(jù)重采樣到 30 米后出現(xiàn)了一些原始數(shù)據(jù)里沒有的類別值。原因用了雙線性或立方卷積等連續(xù)型重采樣方法把類別值當(dāng)連續(xù)值插值了。解決分類柵格的重采樣必須用最近鄰法nearest neighbor。gdalwarp里加-r nearPython 里用rasterio的Resampling.nearest。如果確實(shí)需要降分辨率也可以先用眾數(shù)濾波再抽樣效果比直接最近鄰更好。5. 進(jìn)階技巧用土地覆蓋數(shù)據(jù)做生態(tài)質(zhì)量評(píng)價(jià)的完整鏈路5.1 從土地覆蓋到生態(tài)指數(shù)土地覆蓋數(shù)據(jù)本身是基礎(chǔ)底圖真正體現(xiàn)價(jià)值的是用它算出衍生指標(biāo)。一個(gè)經(jīng)典應(yīng)用是計(jì)算區(qū)域生態(tài)質(zhì)量指數(shù)比如基于土地利用的景觀生態(tài)風(fēng)險(xiǎn)指數(shù)或者生境質(zhì)量模型。這里給一個(gè)簡化但可復(fù)現(xiàn)的鏈路用土地覆蓋計(jì)算香農(nóng)多樣性指數(shù)再結(jié)合不透水面比例做生態(tài)質(zhì)量分級(jí)。import rasterio import numpy as np from scipy.stats import entropy with rasterio.open(gz_landcover_reclass.tif) as src: data src.read(1) profile src.profile # 滑動(dòng)窗口計(jì)算香農(nóng)多樣性指數(shù) window_size 50 # 500米窗口 rows, cols data.shape shannon np.zeros_like(data, dtypenp.float32) for i in range(0, rows, window_size): for j in range(0, cols, window_size): block data[i:iwindow_size, j:jwindow_size] if block.size 0: continue # 統(tǒng)計(jì)各類別出現(xiàn)頻率 unique, counts np.unique(block, return_countsTrue) # 排除 NoData mask unique ! 0 if np.sum(mask) 0: continue probs counts[mask] / np.sum(counts[mask]) shannon[i:iwindow_size, j:jwindow_size] entropy(probs) profile.update(dtyperasterio.float32, nodata-9999) with rasterio.open(gz_shannon.tif, w, **profile) as dst: dst.write(shannon, 1)這段代碼用 50×50 像元的滑動(dòng)窗口計(jì)算香農(nóng)多樣性指數(shù)窗口約 500 米見方。entropy函數(shù)來自 scipy輸入是各類別的概率分布。指數(shù)越高說明窗口內(nèi)土地覆蓋類型越豐富生態(tài)多樣性越好。計(jì)算完后可以按自然斷點(diǎn)法分成高、中、低幾檔再和不透水面比例做疊加識(shí)別出「高多樣性但高不透水面」的沖突區(qū)域這些往往是城市擴(kuò)張前沿值得重點(diǎn)關(guān)注。5.2 驗(yàn)證分類精度的一個(gè)土辦法如果你手頭沒有獨(dú)立的驗(yàn)證樣本但又想對(duì)數(shù)據(jù)質(zhì)量有個(gè)基本判斷可以用一個(gè)土辦法找?guī)追瑫r(shí)期的高分影像截圖在土地覆蓋柵格上隨機(jī)撒點(diǎn)目視比對(duì)。具體操作是在 QGIS 里加載土地覆蓋和影像底圖創(chuàng)建 100 個(gè)隨機(jī)點(diǎn)逐個(gè)查看點(diǎn)所在位置的影像地物和柵格類別是否一致。如果一致率低于 80%說明這份數(shù)據(jù)在你的研究區(qū)可能不太靠譜需要謹(jǐn)慎使用或者考慮自己重新分類。這個(gè)方法雖然粗糙但比盲目信任數(shù)據(jù)強(qiáng)。我一般會(huì)在正式分析前花半小時(shí)做這個(gè)抽查翻車過幾次之后就養(yǎng)成了習(xí)慣。5.3 數(shù)據(jù)融合的邊界最后說一個(gè)容易上頭的地方有人會(huì)想把這份 10 米土地覆蓋和夜間燈光、POI、路網(wǎng)數(shù)據(jù)全部融合起來做城市邊界提取。思路沒錯(cuò)但要注意尺度匹配。夜間燈光數(shù)據(jù)常見的是 500 米分辨率POI 是點(diǎn)數(shù)據(jù)路網(wǎng)是矢量線直接和 10 米柵格疊加會(huì)產(chǎn)生大量空值和尺度不一致問題。常見做法是先把所有數(shù)據(jù)統(tǒng)一到同一個(gè)格網(wǎng)比如 100 米或 250 米再做多要素加權(quán)。不要試圖在 10 米尺度上融合所有數(shù)據(jù)那樣計(jì)算量大不說精度提升也有限。從那以后我每次做多源融合之前都強(qiáng)制先跑一遍尺度一致性檢查確認(rèn)所有數(shù)據(jù)的空間支撐匹配了再往下走。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取