简介：本资源为2020年广东省10米分辨率土地覆盖与土地利用数据包，面向地理信息、遥感分析、城市规划及生态环境研究等领域的从业者与学习者，可解决省级、市级尺度土地利用现状提取与空间分析的数据需求。数据基于10米哨兵影像，采用深度学习方法制作，涵盖耕地、林地、草地、灌木、湿地、水体、不透水面、裸地、雪冰等十类地物，并已由墨卡托投影转为WGS84地理坐标系，按最新省市级行政边界裁剪，每个地市均提供独立TIF成果。压缩包共147个文件，约75.8MB，包含21个tif主数据、21个xlsx属性表、21个tfw坐标文件、21个xml元数据、21个cpg编码文件、21个dbf属性库及21个png预览图，配套完整便于直接加载与制图。目前已有329人学习下载，适合用于土地利用变化监测、城市扩张分析、生态评估及GIS课程实践等场景。
1. 拿到「2020年10m精度广东省土地覆盖土地利用.rar」先别急着解压：这份数据到底能干什么
如果你手头正好有一份「2020年10m精度广东省土地覆盖土地利用.rar」，大概率是冲着两件事来的：一是想给广东某个市、某个流域、某个项目做一份像样的地表覆盖底图，二是想拿它当训练样本或验证底图，去跑自己的遥感分类模型。10m 分辨率意味着一个像素对应地面 100 平方米，这个尺度刚好卡在「能看清地块、又不至于让文件大到打不开」的甜点区，对省级尺度的耕地监测、城市扩张分析、生态红线核查都够用。但我要先泼一盆冷水：这类打包数据最常见的翻车点，不是精度不够，而是坐标系、分类体系、NoData 值三件事没对齐，直接拿去算面积能给你算出离谱结果。这篇就按「先搞懂它是什么 → 再动手跑通 → 最后避开那些血泪坑」的顺序讲清楚，新手能照着复现，熟手能直接跳到参数和边界那几节。
2. 先搞清楚这份广东土地覆盖数据的技术底子：分类体系、坐标系与文件组织
2.1 10m 分辨率在省级尺度上意味着什么
先把量级算清楚，不然后面所有操作都是玄学。广东省陆域面积约 17.97 万平方公里，换算成平方米是 1.797×10¹¹ 平方米。10m 分辨率下每个像素 100 平方米，那么整省满打满算约 1.797×10⁹ 个像素，也就是接近 18 亿个像元。这个数字决定了你后面所有工具的选择：用 Python 的 rasterio 分块读没问题，但想一次性read()进内存再转 numpy 数组做全局运算，普通 16G 内存的机器会直接给你脸色看。所以从第一步就要建立「分块处理」的意识，这不是优化技巧，是能不能跑通的前提。
10m 这个精度在国内土地覆盖产品里属于中高分辨率档位。常见的 30m 产品（比如早期的 GlobeLand30 系列思路）在省级尺度上做宏观统计够用，但一旦你要看珠三角那种村镇级别的建设用地碎片，30m 就会把很多小地块糊成一团。10m 能明显改善这一点，代价是数据量和计算量大约涨了 9 倍。理解这个取舍，你才知道为什么这份数据值得单独拿出来做，而不是随便找个 30m 的凑合。
2.2 分类体系：先确认它用的是哪套编码
这是最容易埋雷的地方。国内土地覆盖数据常见的分类体系有两套思路：一套是面向遥感解译的「一级大类 + 二级细分」，比如耕地、林地、草地、水体、建设用地、未利用地这六大类往下再分；另一套是面向国土调查的「三调」式分类，编码更细、更贴近管理口径。你拿到手的这份数据到底用哪套，不能猜，必须打开属性表或配套说明确认。
我一般的做法是先看栅格值的分布，用一段极短的代码把唯一值和计数打出来，一眼就能判断它是不是连续编码、有没有异常值：
import rasterio import numpy as np # 打开栅格，只读元数据，不读全部像素 with rasterio.open("gd_landcover_2020_10m.tif") as src: print("尺寸:", src.width, "x", src.height) print("波段数:", src.count) print("坐标系:", src.crs) print("像素大小:", src.res) print("NoData:", src.nodata) # 分块统计唯一值，避免一次性读入 uniq, counts = np.unique( src.read(1, out_shape=(1, src.height // 10, src.width // 10)), return_counts=True ) for v, c in zip(uniq, counts): print(f"值 {v}: {c} 个采样像元")这段代码的关键在out_shape参数，它让 rasterio 做降采样读取，把整幅图缩到十分之一再统计，速度快几十倍，用来快速摸清分类编码足够。逻辑说明：src.crs告诉你坐标系，src.res告诉你实际像素分辨率（有些数据标称 10m 但重采样过，实际 res 可能是 9.8 或 10.2），src.nodata是后面算面积时必须排除的值。参数说明：out_shape里的height // 10是降采样倍数，你想更精确可以改成// 5，但别直接去掉，否则大图会卡死。
2.3 坐标系与投影：面积计算前必须确认的事
如果这份数据是地理坐标系（比如 CGCS2000 经纬度，单位是度），那你直接算面积就是错的，因为经纬度下每个像素代表的实际面积随纬度变化。广东省跨纬度约 20°N 到 25°N，这个范围内 10m 像素的实际面积差异虽然不算巨大，但做省级统计时累积误差不能忽略。正确做法是先投影到等面积投影或高斯-克吕格投影，再算面积。
判断方法很简单，看src.crs.is_geographic是 True 还是 False。如果是 True，你需要先重投影。常见做法是用 CGCS2000 的 3 度带高斯投影，广东大致落在 38 带（中央经线 114°E）附近，但跨带的话要按实际范围选。这一步别偷懒，我见过太多人拿着经纬度数据直接像素数 × 100算面积，最后报出去的耕地面积偏差百分之几，被甲方一问就露馅。
3. 从 .rar 到可用栅格：解压、校验、裁剪到目标区域的完整流程
3.1 解压后先做文件清单和完整性校验
拿到 .rar 别急着双击解压到桌面。先确认里面到底是什么结构：是单个大 tif，还是按地级市切好的分幅，还是带了 shp 辅助文件。用命令行列一下清单最稳妥：
# 列出压缩包内容，不解压 unrar l "2020年10m精度广东省土地覆盖土地利用.rar" # 解压到指定目录，保留目录结构 unrar x "2020年10m精度广东省土地覆盖土地利用.rar" ./gd_lulc_2020/逻辑说明：l是 list，只查看不落地，先看清楚有几个文件、多大、什么格式，再决定解压策略。x是带完整路径解压，避免文件散落一地。参数说明：如果压缩包有密码，unrar会提示输入；如果文件名含中文，确保终端编码是 UTF-8，否则可能解压出乱码文件名。
解压完做一次校验：用gdalinfo看栅格是否完整、有没有损坏。如果gdalinfo报错说无法读取，那多半是解压不完整或文件本身有问题，这时候别硬着头皮往下做，先重新解压或找原始来源核对。
3.2 用 GDAL 做重投影和裁剪：一条命令解决
假设你已经确认数据是经纬度坐标系，现在要投影并裁剪到某个市，比如广州市。分两步走，先投影再裁剪，或者用gdalwarp一步到位：
# 一步完成：重投影到高斯投影 + 按广州边界裁剪 gdalwarp -t_srs "EPSG:4547" \ -cutline guangzhou_boundary.shp \ -crop_to_cutline \ -dstnodata 0 \ -co COMPRESS=LZW \ -co TILED=YES \ gd_landcover_2020_10m.tif \ gz_landcover_2020_10m.tif逻辑说明：-t_srs "EPSG:4547"指定目标投影，EPSG:4547 是 CGCS2000 3 度带 39 带（中央经线 117°E），广州大致在这个带附近，具体用哪个带要看你的数据范围，选错了会有投影变形。-cutline指定裁剪边界 shp，-crop_to_cutline让输出范围严格贴合边界。-dstnodata 0把裁剪外的区域设为 0，方便后续排除。-co COMPRESS=LZW用 LZW 压缩，能显著减小文件体积，-co TILED=YES让输出带内部瓦片结构，后续分块读取更快。
参数说明：如果你的边界 shp 和栅格坐标系不一致，gdalwarp会自动做动态投影，但前提是 shp 有正确的 .prj 文件。没有 .prj 的话，它会按栅格坐标系硬套，结果可能整体偏移，这是很隐蔽的坑。裁剪完记得再用gdalinfo确认一下输出范围、像素大小、NoData 值是否符合预期。
3.3 用 Python 做分块统计：算各类面积
裁剪完的小区域就可以用 Python 处理了。算面积的核心逻辑是：统计每个类别的像素数，乘以单个像素的实际面积。注意，投影后的像素面积是固定的（比如 10m × 10m = 100 平方米），但前提是你已经投影到等面积或局部投影，且像素是正方形。
import rasterio import numpy as np with rasterio.open("gz_landcover_2020_10m.tif") as src: nodata = src.nodata pixel_area = abs(src.res[0] * src.res[1]) # 单像素面积，平方米 # 分块读取并累计各类别像素数 counts = {} for _, window in src.block_windows(1): block = src.read(1, window=window) # 排除 NoData valid = block[block != nodata] if nodata is not None else block uniq, cnt = np.unique(valid, return_counts=True) for v, c in zip(uniq, cnt): counts[v] = counts.get(v, 0) + int(c) # 输出面积，单位：公顷 for cls, cnt in sorted(counts.items()): area_ha = cnt * pixel_area / 10000 print(f"类别 {cls}: {area_ha:.2f} 公顷")逻辑说明：block_windows(1)按栅格内部瓦片逐块读取，内存占用可控。src.res返回像素的 x、y 方向尺寸，投影后通常是 (10.0, 10.0)，相乘得单像素面积。np.unique统计每块内各类别数量，累加到全局字典。参数说明：pixel_area单位是平方米，除以 10000 转公顷。如果你的数据是经纬度且没投影，src.res会是度为单位的小数，这时候算出来的「面积」没有物理意义，必须回到上一步先投影。
4. 避坑与排查：这份广东土地覆盖数据最容易翻车的 5 个地方
4.1 现象：算出来的面积比官方统计大了一圈
原因：NoData 值没排除干净。很多土地覆盖数据用 0 表示 NoData，但 0 也可能被某些分类体系用作「其他」或「未分类」的有效值。如果你直接block != 0排除，可能把有效类别也排掉了，或者反过来，NoData 被当成有效值参与统计。
解决：先确认src.nodata的返回值。如果它是 None，说明数据没定义 NoData，你需要根据数据说明手动指定。如果它是 0，但分类编码里 0 也是有效类，那就得用掩膜文件或边界 shp 来限定统计范围，而不是靠值排除。
4.2 现象：裁剪后的图和边界对不上，整体偏移几百米
原因：边界 shp 缺少 .prj 投影文件，或者 shp 和栅格的坐标系定义不一致但被强行套用。
解决：用ogrinfo检查 shp 的坐标系，和gdalinfo输出的栅格坐标系对比。如果不一致，先用ogr2ogr把 shp 转到和栅格相同的坐标系，再做裁剪。别指望gdalwarp每次都猜对。
4.3 现象：分块统计时结果和一次性读取不一致
原因：分块读取时，块与块之间的边界像素可能被重复计算或遗漏，尤其是当block_windows的窗口有重叠时。
解决：rasterio 的block_windows默认不重叠，但如果你手动指定了带重叠的窗口，就要在累加时去重。更稳妥的做法是用src.read(1)配合out_shape降采样做快速估算，再用分块做精确统计，两者交叉验证。
4.4 现象：重投影后像素大小不是 10m 了
原因：gdalwarp默认会保持输出像素大小和输入一致，但如果输入是经纬度、输出是投影坐标系，它会自动计算一个近似值，可能变成 9.8m 或 10.3m。
解决：显式指定-tr 10 10强制输出 10m 像素。但要注意，强制指定可能导致重采样引入误差，如果原始数据精度要求高，建议用-r near最近邻重采样，保持类别值不变。
4.5 现象：打开 tif 一片黑或一片白
原因：分类栅格的值域很小（比如 1 到 6），但显示时按 0-255 拉伸，导致所有类别挤在一起看不出区别。
解决：在 QGIS 或 ArcGIS 里手动设置符号系统的值域范围，或者用gdaldem color-relief生成彩色渲染图。这不是数据问题，是显示设置问题，别误以为数据坏了。
5. 进阶用法：把这份数据接进自己的分类模型或时序分析
5.1 当训练标签用：采样与格式转换
如果你想拿这份 2020 年的土地覆盖数据当训练标签，去训练自己的遥感分类模型，第一步是采样。别全图采样，那样样本极度不均衡（林地可能占一半以上）。按类别分层采样，每个类别抽几千个点，导出成模型能吃的格式。
import rasterio import numpy as np from sklearn.model_selection import train_test_split with rasterio.open("gz_landcover_2020_10m.tif") as src: data = src.read(1) nodata = src.nodata # 展平并排除 NoData flat = data.flatten() if nodata is not None: flat = flat[flat != nodata] # 分层采样：每个类别最多抽 5000 个 samples = [] labels = [] for cls in np.unique(flat): cls_pixels = flat[flat == cls] n = min(len(cls_pixels), 5000) samples.extend(np.random.choice(cls_pixels, n, replace=False)) labels.extend([cls] * n) # 划分训练集和验证集 X_train, X_val, y_train, y_val = train_test_split( samples, labels, test_size=0.2, stratify=labels, random_state=42 ) print(f"训练集 {len(X_train)}，验证集 {len(X_val)}")逻辑说明：np.random.choice做无放回采样，stratify=labels保证训练集和验证集的类别比例一致。参数说明：5000是每个类别的采样上限，类别多的可以调大，类别少的会自动取全部。random_state=42固定随机种子，保证可复现。
5.2 做时序变化分析：和 2010、2015 数据对齐
如果你手头还有 2010、2015 的同类数据，想做广东十年土地利用变化，核心是对齐。不同年份的数据可能来自不同传感器、不同分类体系，直接相减会得到一堆无意义的「变化」。常见做法是先把各年份的分类体系统一到同一套编码，再做逐像素对比。
对齐时注意三点：一是空间范围要完全一致，用同一份边界 shp 裁剪；二是像素网格要对齐，用gdalwarp的-te和-tr参数强制统一；三是分类编码要建立映射表，比如把 2010 的「耕地」编码映射到 2020 的对应编码。做完这三步，再算转移矩阵才有意义。
5.3 一个我常用的验证习惯
每次处理完这类数据，我都会做一件事：随机抽 100 个点，用高分辨率影像（比如天地图或谷歌地球的历史影像）人工核对类别。不用多，100 个点就能大致判断这份数据的可信度。如果准确率低于 80%，那后面的分析结论都要打问号。这个习惯帮我避开了好几次「数据看着漂亮、结论完全错误」的尴尬。土地覆盖数据不是拿来就用的一次性消耗品，它的分类体系、坐标系、NoData 定义决定了你能做什么、不能做什么。把这几件事在前两章确认清楚，后面所有操作才有意义。希望帮到你。
本文还有配套的精品资源，点击获取
