簡介這份資源為2020年廣東省10米精度土地覆蓋土地利用數(shù)據(jù)面向地理信息、遙感、城市規(guī)劃及生態(tài)環(huán)境研究等方向的學(xué)習(xí)者與從業(yè)者。數(shù)據(jù)基于哨兵影像與深度學(xué)習(xí)方法制作經(jīng)墨卡托投影轉(zhuǎn)WGS84地理坐標(biāo)系并按最新省市級行政邊界裁剪覆蓋廣東省各地級市共分耕地、林地、草地、灌木、濕地、水體、不透水面、裸地、雪冰等十類。壓縮包共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ī)范的廣東省市級土地利用數(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)評估、城市擴(kuò)張監(jiān)測完全夠用。它適合做自然資源調(diào)查的技術(shù)人員、做城市研究的規(guī)劃師、寫論文需要土地利用底圖的研究生以及做遙感應(yīng)用開發(fā)需要一份現(xiàn)成分類結(jié)果做驗(yàn)證的工程師。但先別急著解壓這份數(shù)據(jù)能不能用、怎么用、坑在哪得先把它拆開看明白。2. 10米分辨率土地覆蓋柵格從分類體系到廣東地類的對應(yīng)關(guān)系2.1 為什么是10米而不是30米或1米土地覆蓋產(chǎn)品的分辨率選擇從來不是拍腦袋定的。30 米級別的產(chǎn)品比如早期基于 Landsat 的全球地表覆蓋在省級尺度上做宏觀趨勢沒問題但一旦落到珠三角這種城市蔓延劇烈、地塊破碎度高的區(qū)域30 米像素會把一個(gè)城中村和旁邊的工業(yè)園混成一個(gè)類別邊界模糊得讓人抓狂。1 米級別的高分影像當(dāng)然精細(xì)但廣東省全域覆蓋的高分?jǐn)?shù)據(jù)獲取成本高、處理量大而且分類精度受陰影和建筑紋理干擾反而可能下降。10 米恰好卡在一個(gè)甜點(diǎn)上哨兵二號Sentinel-2的多光譜波段原生就是 10 米、20 米和 60 米其中可見光和近紅外是 10 米用這些波段做分類既保留了城市內(nèi)部地塊的區(qū)分度又不至于讓數(shù)據(jù)量爆炸。這份 2020 年的產(chǎn)品大概率就是基于哨兵二號時(shí)序影像生產(chǎn)的這也是目前省級 10 米土地覆蓋產(chǎn)品的主流技術(shù)路線。2.2 分類體系怎么讀一級類和二級類的映射拿到柵格后第一個(gè)要搞清楚的是它的分類編碼。國內(nèi)常見的土地覆蓋產(chǎn)品通常參考《土地利用現(xiàn)狀分類》國標(biāo)或者類似 FROM-GLC、CLCD 等產(chǎn)品的分類體系。這份數(shù)據(jù)大概率包含耕地、林地、草地、灌木地、濕地、水體、不透水面、裸地這幾大類。你需要找到隨數(shù)據(jù)附帶的說明文件或者屬性表確認(rèn)每個(gè)像元值對應(yīng)什么類別。常見做法是值 1 代表耕地2 代表林地3 代表草地4 代表灌木5 代表濕地6 代表水體7 代表不透水面8 代表裸地。但不同生產(chǎn)單位編碼順序可能不同所以務(wù)必以實(shí)際說明為準(zhǔn)。如果壓縮包里沒有說明文件那就得靠目視比對打開珠三角區(qū)域看哪個(gè)值大面積出現(xiàn)在城市核心區(qū)那個(gè)值大概率就是不透水面。2.3 數(shù)據(jù)格式與坐標(biāo)系統(tǒng)確認(rèn)解壓后你大概率會看到一個(gè)或多個(gè) .tif 文件可能按地級市分幅也可能是一整幅廣東全省的鑲嵌影像。坐標(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ù)說明核對清楚。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 重分類把二級類合并成你需要的一級類原始數(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è)舊值把對應(yīng)像元賦成新值。注意np.zeros_like(data)創(chuàng)建的全零數(shù)組如果原始數(shù)據(jù)里有 0 值代表某種地類這里就會混淆所以務(wù)必確認(rèn) NoData 和有效類別不重疊。profile.update里把數(shù)據(jù)類型改成 uint8因?yàn)橹胤诸惡箢悇e數(shù)變少了用不著原來的 uint16 或 float能省不少存儲空間。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)對齊 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)編號假設(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ù)對齊如果你手頭還有 2010 年的同類型數(shù)據(jù)想做十年變化檢測那第一步是確保兩期數(shù)據(jù)坐標(biāo)系、分辨率、范圍完全一致。常見做法是以 2020 年數(shù)據(jù)為基準(zhǔn)把 2010 年數(shù)據(jù)重采樣和裁剪到相同網(wǎng)格。# 將2010年數(shù)據(jù)對齊到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表示用最近鄰重采樣這對分類數(shù)據(jù)是必須的不能用雙線性或立方卷積否則會出現(xiàn)不存在的類別值。對齊之后兩期數(shù)據(jù)逐像元相減非零的地方就是變化區(qū)域。4. 避坑與排查10米土地覆蓋數(shù)據(jù)常見的五個(gè)翻車點(diǎn)4.1 面積算出來偏大或偏小現(xiàn)象用地理坐標(biāo)系直接統(tǒng)計(jì)面積發(fā)現(xiàn)各地類加總跟官方公布的行政區(qū)面積對不上偏差能到百分之十幾。原因地理坐標(biāo)系下像元面積隨緯度變化高緯度地區(qū)像元實(shí)際面積比低緯度小用固定值乘會引入系統(tǒng)誤差。解決先投影到等面積投影或高斯克呂格投影再做統(tǒng)計(jì)。廣東常用 CGCS2000 3 度帶中央經(jīng)線根據(jù)研究區(qū)經(jīng)度選擇比如 114°E 對應(yīng) EPSG:4547 左右。投影后像元面積恒定統(tǒng)計(jì)結(jié)果才可靠。4.2 分類編碼對不上說明書現(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)象把各地級市分幅數(shù)據(jù)鑲嵌后發(fā)現(xiàn)市界處有一條明顯的類別突變線。原因分幅生產(chǎn)時(shí)各幅獨(dú)立分類接邊處分類結(jié)果不一致或者鑲嵌時(shí)沒有做羽化處理。解決如果數(shù)據(jù)是分幅的優(yōu)先找全省整幅版本。如果沒有鑲嵌后用眾數(shù)濾波做一下后處理或者在接受范圍內(nèi)忽略接邊處的少量誤差。做變化檢測時(shí)尤其要注意接邊處的假變化會干擾結(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 會出現(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ì)量評價(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ì)量分級。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)法分成高、中、低幾檔再和不透水面比例做疊加識別出「高多樣性但高不透水面」的沖突區(qū)域這些往往是城市擴(kuò)張前沿值得重點(diǎn)關(guān)注。5.2 驗(yàn)證分類精度的一個(gè)土辦法如果你手頭沒有獨(dú)立的驗(yàn)證樣本但又想對數(shù)據(jù)質(zhì)量有個(gè)基本判斷可以用一個(gè)土辦法找?guī)追瑫r(shí)期的高分影像截圖在土地覆蓋柵格上隨機(jī)撒點(diǎn)目視比對。具體操作是在 QGIS 里加載土地覆蓋和影像底圖創(chuàng)建 100 個(gè)隨機(jī)點(diǎn)逐個(gè)查看點(diǎn)所在位置的影像地物和柵格類別是否一致。如果一致率低于 80%說明這份數(shù)據(jù)在你的研究區(qū)可能不太靠譜需要謹(jǐn)慎使用或者考慮自己重新分類。這個(gè)方法雖然粗糙但比盲目信任數(shù)據(jù)強(qiáng)。我一般會在正式分析前花半小時(shí)做這個(gè)抽查翻車過幾次之后就養(yǎng)成了習(xí)慣。5.3 數(shù)據(jù)融合的邊界最后說一個(gè)容易上頭的地方有人會想把這份 10 米土地覆蓋和夜間燈光、POI、路網(wǎng)數(shù)據(jù)全部融合起來做城市邊界提取。思路沒錯(cuò)但要注意尺度匹配。夜間燈光數(shù)據(jù)常見的是 500 米分辨率POI 是點(diǎn)數(shù)據(jù)路網(wǎng)是矢量線直接和 10 米柵格疊加會產(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)擊獲取