拓冰网站定制 
简介本资源为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_countsTrue ) for v, c in zip(uniq, counts): print(f值 {v}: {c} 个采样像元)这段代码的关键在out_shape参数它让 rasterio 做降采样读取把整幅图缩到十分之一再统计速度快几十倍用来快速摸清分类编码足够。逻辑说明src.crs告诉你坐标系src.res告诉你实际像素分辨率有些数据标称 10m 但重采样过实际 res 可能是 9.8 或 10.2src.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 COMPRESSLZW \ -co TILEDYES \ 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 COMPRESSLZW用 LZW 压缩能显著减小文件体积-co TILEDYES让输出带内部瓦片结构后续分块读取更快。参数说明如果你的边界 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, windowwindow) # 排除 NoData valid block[block ! nodata] if nodata is not None else block uniq, cnt np.unique(valid, return_countsTrue) 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, replaceFalse)) labels.extend([cls] * n) # 划分训练集和验证集 X_train, X_val, y_train, y_val train_test_split( samples, labels, test_size0.2, stratifylabels, random_state42 ) print(f训练集 {len(X_train)}验证集 {len(X_val)})逻辑说明np.random.choice做无放回采样stratifylabels保证训练集和验证集的类别比例一致。参数说明5000是每个类别的采样上限类别多的可以调大类别少的会自动取全部。random_state42固定随机种子保证可复现。5.2 做时序变化分析和 2010、2015 数据对齐如果你手头还有 2010、2015 的同类数据想做广东十年土地利用变化核心是对齐。不同年份的数据可能来自不同传感器、不同分类体系直接相减会得到一堆无意义的「变化」。常见做法是先把各年份的分类体系统一到同一套编码再做逐像素对比。对齐时注意三点一是空间范围要完全一致用同一份边界 shp 裁剪二是像素网格要对齐用gdalwarp的-te和-tr参数强制统一三是分类编码要建立映射表比如把 2010 的「耕地」编码映射到 2020 的对应编码。做完这三步再算转移矩阵才有意义。5.3 一个我常用的验证习惯每次处理完这类数据我都会做一件事随机抽 100 个点用高分辨率影像比如天地图或谷歌地球的历史影像人工核对类别。不用多100 个点就能大致判断这份数据的可信度。如果准确率低于 80%那后面的分析结论都要打问号。这个习惯帮我避开了好几次「数据看着漂亮、结论完全错误」的尴尬。土地覆盖数据不是拿来就用的一次性消耗品它的分类体系、坐标系、NoData 定义决定了你能做什么、不能做什么。把这几件事在前两章确认清楚后面所有操作才有意义。希望帮到你。本文还有配套的精品资源点击获取 2020年10m精度广东省土地覆盖数据：从解压到模型训练全流程避坑指南
简介本资源为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_countsTrue ) for v, c in zip(uniq, counts): print(f值 {v}: {c} 个采样像元)这段代码的关键在out_shape参数它让 rasterio 做降采样读取把整幅图缩到十分之一再统计速度快几十倍用来快速摸清分类编码足够。逻辑说明src.crs告诉你坐标系src.res告诉你实际像素分辨率有些数据标称 10m 但重采样过实际 res 可能是 9.8 或 10.2src.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 COMPRESSLZW \ -co TILEDYES \ 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 COMPRESSLZW用 LZW 压缩能显著减小文件体积-co TILEDYES让输出带内部瓦片结构后续分块读取更快。参数说明如果你的边界 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, windowwindow) # 排除 NoData valid block[block ! nodata] if nodata is not None else block uniq, cnt np.unique(valid, return_countsTrue) 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, replaceFalse)) labels.extend([cls] * n) # 划分训练集和验证集 X_train, X_val, y_train, y_val train_test_split( samples, labels, test_size0.2, stratifylabels, random_state42 ) print(f训练集 {len(X_train)}验证集 {len(X_val)})逻辑说明np.random.choice做无放回采样stratifylabels保证训练集和验证集的类别比例一致。参数说明5000是每个类别的采样上限类别多的可以调大类别少的会自动取全部。random_state42固定随机种子保证可复现。5.2 做时序变化分析和 2010、2015 数据对齐如果你手头还有 2010、2015 的同类数据想做广东十年土地利用变化核心是对齐。不同年份的数据可能来自不同传感器、不同分类体系直接相减会得到一堆无意义的「变化」。常见做法是先把各年份的分类体系统一到同一套编码再做逐像素对比。对齐时注意三点一是空间范围要完全一致用同一份边界 shp 裁剪二是像素网格要对齐用gdalwarp的-te和-tr参数强制统一三是分类编码要建立映射表比如把 2010 的「耕地」编码映射到 2020 的对应编码。做完这三步再算转移矩阵才有意义。5.3 一个我常用的验证习惯每次处理完这类数据我都会做一件事随机抽 100 个点用高分辨率影像比如天地图或谷歌地球的历史影像人工核对类别。不用多100 个点就能大致判断这份数据的可信度。如果准确率低于 80%那后面的分析结论都要打问号。这个习惯帮我避开了好几次「数据看着漂亮、结论完全错误」的尴尬。土地覆盖数据不是拿来就用的一次性消耗品它的分类体系、坐标系、NoData 定义决定了你能做什么、不能做什么。把这几件事在前两章确认清楚后面所有操作才有意义。希望帮到你。本文还有配套的精品资源点击获取