简介：这份广东Landcover数据面向GIS从业者、环境与城市规划研究者及高校师生，提供2020年ESRI发布的10米分辨率土地覆盖栅格成果，可用于地表分类制图、生态评估与空间叠加分析。资源包共16个文件，约163.29MB，以tif栅格与tfw坐标辅助文件为核心，辅以shp、shx、dbf等矢量边界文件，prj与xml记录投影和元数据，ovr、sbx、sbn等则服务于影像金字塔与空间索引，兼顾UTM投影与WGS地理两套坐标系。目前已有274人学习下载。数据覆盖森林、农田、城市、水域等多种地类，配合广东省行政边界矢量，可完成重分类、裁剪、缓冲区与叠加统计等操作，为区域生态研究、城市扩张监测和自然资源管理提供高精度底图支撑。
1. 广东Landcover数据（10m分辨率）：从一张地表覆盖图到能跑的分析底图
做广东区域的项目，绕不开一件事：手里得有一张能对得上地块、对得上影像、对得上行政边界的土地覆盖底图。广东Landcover数据（10m分辨率）就是干这个的——它把地表分成耕地、林地、草地、水体、建设用地等类型，每个像元代表地面10米×10米的范围。这个分辨率意味着什么？一条双向四车道的公路、一片几十亩的果园、一个镇级工业园，都能在图上被单独识别出来，而不是被糊成一团。对做国土空间规划、生态评估、城市扩张监测、农业估产的人来说，这张图是分析底图，不是装饰品。但拿到数据只是开始，真正决定成败的是：坐标系对不对、分类体系能不能和你的业务口径对齐、边界裁切有没有把海岸线和岛屿搞丢。这篇笔记就按我实际用这套数据的顺序，把选型理由、处理步骤、参数设置和踩过的坑讲清楚，让第一次接触的人能照着跑通，也让用过的人看到边界在哪。
2. 广东Landcover数据到底怎么来的：分类体系与精度边界
2.1 10m分辨率背后的数据源与分类逻辑
广东Landcover数据（10m分辨率）的常见生产路径，是以 Sentinel-2 多光谱影像为主要数据源，结合时序特征和地形辅助数据，用监督分类或深度学习模型生成。Sentinel-2 的 10m 空间分辨率对应四个波段：蓝、绿、红、近红外。这四波段是分类的核心输入，因为植被在近红外波段的反射特征和建设用地、水体的差异非常明显。分类体系通常参照国内土地利用现状分类标准，但不同生产方会有细微差别，常见的一级类包括：耕地、林地、草地、灌木、水体、建设用地、裸地、湿地等。
这里有一个容易被忽略的点：10m 分辨率不等于 10m 精度。分辨率是像元大小，精度是分类正确的比例。广东地形复杂，珠三角城市群建筑密集、粤北山区地形阴影重、粤东粤西海岸带水体边界破碎，这些区域的实际分类精度会明显低于平原农田区。我一般会先看数据生产方提供的精度评价报告，如果没有，就自己抽验证点做混淆矩阵。常见做法是每类抽 50 到 100 个点，用高分辨率影像或实地调查数据做参考，算总体精度和 Kappa 系数。总体精度低于 80% 的数据，在做精细地块分析时要谨慎。
2.2 坐标系与投影：为什么你拿到的数据可能对不上
广东Landcover数据常见的坐标系有两种：WGS84 地理坐标系（EPSG:4326）和 CGCS2000 投影坐标系（如 EPSG:4547，CGCS2000 / 3-degree Gauss-Kruger CM 114E）。如果你要做面积统计，必须用投影坐标系，因为地理坐标系下算面积会随纬度变化产生形变。广东跨了约 3 个经度带，用 EPSG:4547 能覆盖大部分区域，但粤东部分可能需要换带。
我遇到过最典型的问题：把 WGS84 的数据直接和 CGCS2000 的行政边界叠加，发现整体偏移了几十米。这不是数据错了，是坐标系没统一。WGS84 和 CGCS2000 在广东区域的差异通常在厘米到米级，但如果你用的边界数据是更老的北京54或西安80，偏移可能达到上百米。处理前先确认两件事：数据的坐标系是什么，你要叠加的其它数据坐标系是什么。不一致就先做投影转换，再裁切。
import rasterio from rasterio.warp import calculate_default_transform, reproject, Resampling # 查看原始数据坐标系和范围 with rasterio.open('guangdong_landcover_10m.tif') as src: print('CRS:', src.crs) print('Bounds:', src.bounds) print('Shape:', src.shape) print('Resolution:', src.res) # 将数据重投影到 CGCS2000 3度带 EPSG:4547 dst_crs = 'EPSG:4547' with rasterio.open('guangdong_landcover_10m.tif') as src: transform, width, height = calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds) kwargs = src.meta.copy() kwargs.update({ 'crs': dst_crs, 'transform': transform, 'width': width, 'height': height }) with rasterio.open('guangdong_landcover_4547.tif', 'w', **kwargs) as dst: for i in range(1, src.count + 1): reproject( source=rasterio.band(src, i), destination=rasterio.band(dst, i), src_transform=src.transform, src_crs=src.crs, dst_transform=transform, dst_crs=dst_crs, resampling=Resampling.nearest # 分类数据必须用最近邻 )这段代码做了三件事：读取原始数据的元信息、计算目标投影下的变换参数、逐波段重投影。关键参数是Resampling.nearest，分类数据只能用最近邻重采样，用双线性或三次卷积会把类别值插值成无意义的数字。重投影后要检查边界是否完整，特别是海岸线和岛屿区域，有时候重投影会把边缘像元裁掉。
2.3 分类体系对齐：你的业务口径和数据口径差在哪
不同来源的 Landcover 数据分类体系不一样。有的把果园归到耕地，有的归到林地；有的把养殖水面单独列出来，有的并到水体。如果你的业务需要统计“耕地面积”，而数据里果园算林地，结果就会偏小。我一般会先做分类体系映射表，把数据的原始类别码映射到业务需要的类别。
| 原始类别码 | 原始类别名 | 业务类别 | 处理方式 |
|---|---|---|---|
| 1 | 水田 | 耕地 | 保留 |
| 2 | 旱地 | 耕地 | 保留 |
| 3 | 有林地 | 林地 | 保留 |
| 4 | 灌木林地 | 林地 | 保留 |
| 5 | 果园 | 园地 | 单独统计 |
| 6 | 水体 | 水体 | 保留 |
| 7 | 建设用地 | 建设用地 | 保留 |
映射完成后，用栅格计算器或 Python 做重分类。注意：重分类后的数据要保留原始数据的空间参考和像元对齐关系，否则后续做变化检测时会对不上。
3. 用 Python 把广东Landcover数据跑成可用底图：裁切、重分类与统计
3.1 按行政边界裁切：别把海岸线和岛屿裁丢了
拿到全省数据后，第一步通常是裁切到研究区。常见做法是用行政边界矢量做掩膜提取。这里有个坑：广东海岸线曲折，岛屿多，如果边界矢量本身精度不够，或者裁切时用了简单的矩形范围，会把沿海的滩涂、红树林、岛屿直接切掉。我一般会先检查边界矢量的拓扑质量，确保没有自相交和缝隙，然后用rasterio.mask做精确裁切。
import geopandas as gpd import rasterio from rasterio.mask import mask import numpy as np # 读取行政边界 gdf = gpd.read_file('guangdong_boundary.shp') # 确保边界和栅格坐标系一致 with rasterio.open('guangdong_landcover_4547.tif') as src: if gdf.crs != src.crs: gdf = gdf.to_crs(src.crs) # 合并所有几何体作为掩膜 geoms = gdf.geometry.values out_image, out_transform = mask(src, geoms, crop=True, nodata=0) out_meta = src.meta.copy() out_meta.update({ 'height': out_image.shape[1], 'width': out_image.shape[2], 'transform': out_transform, 'nodata': 0 }) with rasterio.open('guangdong_clipped.tif', 'w', **out_meta) as dst: dst.write(out_image)crop=True会把输出范围收紧到边界的最小外接矩形，减少文件体积。nodata=0把边界外的区域设为 0，后续统计时排除。裁切后一定要打开看一眼，特别是珠江口、雷州半岛、南澳岛这些区域，确认没有异常缺失。
3.2 重分类与面积统计：从像元数到平方公里
裁切完成后，按业务需求做重分类，然后统计各类面积。面积统计的核心公式：面积 = 像元数 × 像元面积。10m 分辨率的像元面积是 100 平方米，即 0.0001 平方公里。但要注意，如果数据经过了重投影，像元大小可能不是精确的 10m，要以实际元数据为准。
import rasterio import numpy as np import pandas as pd # 读取裁切后的数据 with rasterio.open('guangdong_clipped.tif') as src: data = src.read(1) pixel_size = src.res[0] # 像元大小，单位与投影一致（米） pixel_area_km2 = (pixel_size * pixel_size) / 1e6 # 转换为平方公里 # 定义重分类映射 reclass_map = { 1: 1, 2: 1, # 水田、旱地 -> 耕地 3: 2, 4: 2, # 有林地、灌木林地 -> 林地 5: 3, # 果园 -> 园地 6: 4, # 水体 -> 水体 7: 5 # 建设用地 -> 建设用地 } # 应用重分类 reclassed = np.zeros_like(data) for old_val, new_val in reclass_map.items(): reclassed[data == old_val] = new_val # 统计各类面积 class_names = {1: '耕地', 2: '林地', 3: '园地', 4: '水体', 5: '建设用地'} results = [] for class_val, class_name in class_names.items(): count = np.sum(reclassed == class_val) area = count * pixel_area_km2 results.append({'类别': class_name, '像元数': count, '面积(km2)': round(area, 2)}) df = pd.DataFrame(results) print(df)这段代码先读取数据，然后按映射表重分类，最后统计各类像元数和面积。关键参数是pixel_area_km2，它决定了面积统计的基准。如果数据是地理坐标系（度为单位），不能直接用这个公式，必须先投影到米制坐标系。统计结果建议导出为 CSV，方便后续和统计年鉴或其它数据源做交叉验证。
3.3 用 GDAL 命令行做批量处理：适合多县区并行
如果研究区涉及多个县区，逐个用 Python 脚本跑效率低。我一般用 GDAL 命令行做批量裁切和重投影，配合 shell 脚本并行。GDAL 的gdalwarp支持按矢量裁切，gdal_calc.py支持栅格计算。
# 按掩膜裁切，-cutline 指定边界矢量，-crop_to_cutline 收紧范围 gdalwarp -cutline guangdong_boundary.shp \ -crop_to_cutline \ -dstnodata 0 \ -tr 10 10 \ -r near \ guangdong_landcover_4547.tif \ guangdong_clipped_gdal.tif # 用 gdal_calc.py 做重分类，将原始类别映射为新类别 gdal_calc.py -A guangdong_clipped_gdal.tif \ --outfile=guangdong_reclass.tif \ --calc="(A==1)*1+(A==2)*1+(A==3)*2+(A==4)*2+(A==5)*3+(A==6)*4+(A==7)*5" \ --NoDataValue=0 \ --type=Byte-tr 10 10强制输出 10m 像元，-r near指定最近邻重采样。gdal_calc.py的--calc参数用布尔运算实现重分类，逻辑是：如果像元值等于 1，则输出 1；等于 2，也输出 1；以此类推。这种写法比写 Python 循环快得多，适合大批量处理。注意--type=Byte限制输出为 8 位无符号整数，因为重分类后的类别码不会超过 255。
4. 广东Landcover数据避坑与排查：五个让我返工的血泪教训
4.1 坑一：像元对齐导致叠加分析错位
现象：把 Landcover 数据和土地利用矢量叠加，发现矢量地块边界和栅格类别边界错开半个像元，统计面积时出现“一地两类”。
原因：两个数据的栅格原点或像元对齐方式不一致。常见情况是 Landcover 数据以左上角为原点，而你的矢量数据或其它栅格以中心为原点。
解决：用gdalinfo查看数据的 GeoTransform，确认原点坐标。如果不对齐，用gdalwarp的-tap参数强制对齐到目标分辨率网格。或者在 Python 里用rasterio的align功能统一网格。
4.2 坑二：NoData 值被当成类别参与统计
现象：统计面积时，发现某一类面积异常大，检查发现是边界外的 NoData 区域被算进去了。
原因：NoData 值设置不当，或者统计时没有排除 NoData。有些数据用 0 表示 NoData，但 0 也可能是某个类别的编码。
解决：先确认数据的 NoData 值，统计时用data[data != nodata]排除。如果 0 既是 NoData 又是类别码，需要先做掩膜处理，把 NoData 区域单独标记。
4.3 坑三：重投影后面积统计偏差
现象：同一区域，用地理坐标系和投影坐标系分别统计面积，结果差了几个百分点。
原因：地理坐标系下像元大小随纬度变化，直接用度数算面积会引入误差。广东跨纬度约 4 度，误差不可忽略。
解决：面积统计必须在投影坐标系下进行。如果原始数据是地理坐标系，先重投影到 CGCS2000 或 UTM 投影，再做统计。重投影时用最近邻重采样，避免类别值被插值改变。
4.4 坑四：分类体系映射错误导致业务口径偏差
现象：业务需要“林地”面积，但数据里“果园”被归到林地，导致林地面积虚高。
原因：不同数据的分类体系定义不同，直接使用原始类别码统计会偏离业务口径。
解决：先建立分类映射表，明确每个原始类别对应到业务类别的规则。映射表要经过业务方确认，不能自己拍脑袋。映射完成后，用混淆矩阵或抽样验证检查映射结果的合理性。
4.5 坑五：大文件处理内存溢出
现象：用 Python 读取全省 10m 数据时，内存直接爆掉，程序崩溃。
原因：全省 10m 数据量巨大，广东面积约 17.97 万平方公里，10m 像元约 1.8 亿个，单波段 float32 数据超过 700MB，加上中间变量很容易超过内存限制。
解决：用分块读取（rasterio的block_windows）或窗口读取，每次只处理一小块。或者用 GDAL 命令行工具，它们底层做了流式处理，内存占用低。如果必须用 Python，考虑用dask或xarray做分块计算。
5. 把广东Landcover数据用出进阶价值：变化检测与精度验证的实操技巧
5.1 用两期数据做变化检测：从类别变化到空间格局
单期 Landcover 数据能告诉你“现在是什么”，两期数据叠加才能告诉你“变成了什么”。我一般会做两件事：一是类别转移矩阵，二是变化热点图。转移矩阵用 pandas 的crosstab就能算，变化热点图则需要把变化像元单独提取出来，做核密度分析。
import rasterio import numpy as np import pandas as pd import seaborn as sns import matplotlib.pyplot as plt # 读取两期数据 with rasterio.open('guangdong_landcover_2015.tif') as src: data_2015 = src.read(1) with rasterio.open('guangdong_landcover_2020.tif') as src: data_2020 = src.read(1) # 确保两期数据形状一致 assert data_2015.shape == data_2020.shape, '两期数据形状不一致，需要先对齐' # 排除 NoData mask = (data_2015 > 0) & (data_2020 > 0) t1 = data_2015[mask].flatten() t2 = data_2020[mask].flatten() # 计算转移矩阵 transition = pd.crosstab(t1, t2, rownames=['2015'], colnames=['2020']) print(transition) # 可视化 plt.figure(figsize=(10, 8)) sns.heatmap(transition, annot=True, fmt='d', cmap='YlOrRd') plt.title('2015-2020 广东Landcover类别转移矩阵') plt.tight_layout() plt.savefig('transition_matrix.png', dpi=300)这段代码的核心是pd.crosstab，它统计了每个类别从 2015 到 2020 的转移数量。对角线上的数字表示未变化的像元，非对角线表示发生了转移。重点关注“耕地转建设用地”“林地转园地”这类业务敏感的变化。转移矩阵算完后，可以进一步算变化率：变化率 = 转移像元数 / 基期该类像元总数。
5.2 精度验证：没有验证点的替代方案
理想情况下，精度验证需要独立的地面验证点。但实际项目中，验证点往往不够。我一般用两种替代方案：一是用高分辨率影像（如 Google Earth 历史影像）做目视解译抽样，二是用统计数据进行总量控制。目视解译抽样时，每类至少抽 50 个点，用随机分层抽样保证空间分布均匀。总量控制则是把 Landcover 统计的耕地面积和统计年鉴的耕地面积对比，如果偏差超过 10%，就要检查分类体系或数据质量。
| 验证方法 | 所需数据 | 适用场景 | 局限性 |
|---|---|---|---|
| 地面验证点 | GPS 实测点 | 精度要求高、有实地条件 | 成本高、覆盖有限 |
| 高分辨率影像目视解译 | Google Earth 等 | 快速验证、历史数据 | 影像时相不一致 |
| 统计数据总量控制 | 统计年鉴 | 宏观校验 | 统计口径可能不同 |
| 交叉验证 | 其它 Landcover 产品 | 多源对比 | 产品间误差传递 |
5.3 一个我常用的技巧：用众数滤波消除“椒盐噪声”
10m 分辨率的分类结果，在建筑和植被交界处容易出现零星的错分像元，看起来像椒盐噪声。我一般用众数滤波（majority filter）做后处理，窗口大小设为 3×3 或 5×5。GDAL 的gdal_sieve.py可以做类似的事，但它针对的是小图斑。众数滤波更直接：对每个像元，取窗口内出现次数最多的类别值作为输出。
from scipy.ndimage import generic_filter import rasterio import numpy as np def majority_filter(values): """返回窗口内出现次数最多的值""" counts = np.bincount(values.astype(int)) return np.argmax(counts) with rasterio.open('guangdong_reclass.tif') as src: data = src.read(1) profile = src.profile # 应用 3x3 众数滤波 filtered = generic_filter(data, majority_filter, size=3, mode='nearest') # 保存结果 with rasterio.open('guangdong_reclass_filtered.tif', 'w', **profile) as dst: dst.write(filtered, 1)generic_filter的size=3表示 3×3 窗口，mode='nearest'处理边界像元。滤波后要对比滤波前后的面积变化，如果某一类面积变化超过 5%，说明滤波窗口可能太大，把真实的小图斑也滤掉了。我一般先用 3×3 试，效果不够再考虑 5×5，但很少用更大的窗口，因为 10m 数据的小图斑本身就有意义。
做这套数据这些年，我最大的习惯是：拿到数据先不看分类结果，先看元数据、坐标系和 NoData 值。这三样对了，后面的事才顺。希望帮到你。
本文还有配套的精品资源，点击获取
