简介本资源为2020年河南省10米精度土地覆盖与土地利用数据包面向地理信息、遥感、城市规划及生态环境研究等领域的科研人员与学生可解决省级尺度高分辨率地表覆被分析、土地利用变化监测等场景下的数据获取与预处理问题。压缩包共126个文件约65.7MB以18个tif栅格影像为核心配套tfw坐标文件、dbf属性表、cpg编码文件、xml元数据、xlsx统计表与png预览图覆盖河南省各地级市便于直接加载与属性查询。数据基于10米哨兵影像采用深度学习方法制作分为耕地、林地、草地、灌木、湿地、水体、不透水面、裸地、雪冰等类别并已由墨卡托投影转为WGS84地理坐标系按最新省市级行政边界裁剪省去了自行重投影与裁剪的繁琐步骤。目前已有353人学习下载适合需要快速开展河南土地利用制图、空间统计与变化分析的研究者使用。1. 拿到一份 2020 年 10m 分辨率河南省土地覆盖数据先别急着打开你手上如果有一份名为「2020年10m精度河南省土地覆盖土地利用.rar」的压缩包大概率是从公开地理数据平台、科研数据共享站点或者同行手里流转过来的。它本质上是一套栅格数据产品覆盖河南省全域空间分辨率 10 米时间截面是 2020 年每个像元记录一个土地覆盖/土地利用类别码。10 米意味着什么意味着一个像元在地面上是 10m×10m一平方公里有 1 万个像元河南约 16.7 万平方公里粗算下来是十几亿个像元。这个量级决定了它不可能用 Excel 打开也不适合直接扔进 QGIS 里靠鼠标点。这类数据真正解决的问题是你需要知道 2020 年河南每一块地是耕地、林地、草地、水体、建设用地还是未利用地而且要精确到地块级别。做国土空间规划、耕地非农化监测、城市扩张分析、生态红线评估、遥感变化检测的人都会盯上这种 10m 级别的省级产品。它比 30m 的全球产品细比 1m 的商业影像便宜甚至免费是性价比很高的中间档。但拿到 .rar 只是起点真正的坑在于解压之后你看到的可能是一堆 .tif、一个 .prj、一个说明文档甚至还有 .img 或 .dat坐标系、类别编码、无效值全都不一定写在显眼处。这篇就按我实际处理这类数据的顺序把解压、校验、裁剪、统计、踩坑一条线讲清楚。2. 解压之后先做三件事看结构、验坐标系、读类别码2.1 用命令行快速摸清压缩包里的文件构成很多人习惯双击解压然后被一堆同名不同后缀的文件搞晕。我一般先在 Linux 或 WSL 下用unzip -l看一眼列表再决定解压到哪个目录。Windows 上如果没有 unzip用 7-Zip 的命令行版本7z l也一样。# 先列出压缩包内容不实际解压确认文件数量和类型 unzip -l 2020年10m精度河南省土地覆盖土地利用.rar # 如果确认是分幅或分县存放解压到独立目录避免污染当前工作区 mkdir -p henan_lulc_2020 unzip 2020年10m精度河南省土地覆盖土地利用.rar -d henan_lulc_2020 # 进入目录后按扩展名统计快速判断数据组织方式 cd henan_lulc_2020 find . -type f | sed s/.*\.// | sort | uniq -c | sort -nr这段命令的逻辑很直接第一步只读列表防止压缩包里有大量小文件直接铺满当前目录第二步指定-d解压到独立文件夹第三步用find加sed提取扩展名并计数。如果输出里.tif占绝大多数说明是标准栅格分幅如果出现.img那是 ERDAS IMAGINE 格式GDAL 也能读但要注意可能附带.rrd金字塔文件如果出现.dat加.hdr那是 ENVI 格式需要成对保留。参数上唯一要留意的是-d后面不要带中文空格路径里如果有中文建议先改成英文目录再操作否则某些老版本 GDAL 会报编码错误。2.2 用 gdalinfo 确认坐标系、分辨率和无效值解压完别急着写代码先用gdalinfo把第一个 .tif 的元数据打出来。这一步能直接告诉你三件关键事投影是什么、像元大小是不是真的 10m、NoData 值是多少。# 查看单个栅格的完整元数据重点关注 Coordinate System 和 NoData Value gdalinfo henan_lulc_2020/henan_2020_lulc.tif # 如果文件很多批量提取坐标系和分辨率输出成表格方便比对 for f in henan_lulc_2020/*.tif; do echo $f gdalinfo $f | grep -E Pixel Size|Coordinate System|NoData Value donegdalinfo的输出里Pixel Size如果显示(10, 10)说明分辨率达标如果显示(0.0001, 0.0001)且坐标系是地理坐标WGS84那实际地面分辨率会随纬度变化河南大约在 10m 上下浮动需要重投影才能得到严格等面积统计。Coordinate System常见的是WGS 84 / UTM zone 49N或50N河南跨了 49 和 50 两个带如果数据是分幅的很可能一部分在 49N 一部分在 50N合并前必须统一。NoData Value如果是0而类别码里 0 又代表「无数据」那统计时就要显式排除否则会把背景值算成某一类。我见过最坑的情况是 NoData 写成-128但说明文档里没提结果统计出来多出一个莫名其妙的类别。2.3 类别码对照表没有它栅格只是一堆数字10m 土地覆盖产品的类别码通常沿用几套主流体系一套是 6 大类耕地、林地、草地、水体、建设用地、未利用地一套是 8 到 10 类的细分体系还有的会带二级码。压缩包里一般会有一个.txt或.pdf说明但有时候被漏掉。如果确实没有可以按像元值分布反推但更稳妥的是找同源产品的公开文档。常见类别码含义统计时注意事项1耕地含水田旱地部分产品会拆成 11/122林地有林地、灌木林可能合并3草地高覆盖、中覆盖、低覆盖可能合并4水体河流、湖泊、水库注意季节性水体5建设用地不透水面城市、村庄、工矿6未利用地裸地、沙地、盐碱地0 或 255无效值/背景必须排除不能参与面积统计拿到对照表后建议在代码里写成一个字典后续统计直接映射不要每次靠记忆。如果说明文档里的类别码和gdalinfo看到的实际值对不上以实际值分布为准先做一次直方图统计再下结论。3. 用 Python 把整省数据裁成分析单元并做面积统计3.1 环境准备与读取rasterio 比 gdal 更顺手命令行看元数据可以但真正做裁剪和统计我习惯用 Python 的rasterio加numpy。它比直接调 GDAL 的 Python 绑定干净读进来就是数组配合geopandas做矢量边界裁剪很顺。import rasterio from rasterio.mask import mask import geopandas as gpd import numpy as np # 打开整省栅格先看基本信息 src_path henan_lulc_2020/henan_2020_lulc.tif with rasterio.open(src_path) as src: print(CRS:, src.crs) print(分辨率:, src.res) print(波段数:, src.count) print(NoData:, src.nodata) # 读取第一波段为数组方便后续统计 data src.read(1) print(数组形状:, data.shape) print(唯一值:, np.unique(data)[:20])这段代码先打印元数据再把第一波段读成 numpy 数组。src.read(1)里的1是波段索引土地覆盖产品通常只有一个波段所以固定读 1。np.unique只打印前 20 个唯一值防止类别太多刷屏。如果src.nodata是None说明文件里没写无效值需要手动根据直方图判断比如 0 或 255 出现频率异常高那基本就是背景值。参数上唯一要注意的是内存整省 10m 数据读成数组后可能占几个 GB如果机器内存不够不要一次性read改用分块读取后面会讲。3.2 按行政区裁剪用矢量边界切出你要的市或县整省数据直接统计只能得到全省总面积但实际分析往往要落到市、县甚至乡镇。这时候需要一份行政边界矢量用rasterio.mask按几何裁剪。# 读取行政边界矢量假设字段名是 NAME要提取郑州市 admin gpd.read_file(henan_admin.shp) zhengzhou admin[admin[NAME] 郑州市] # 用郑州市边界裁剪栅格 with rasterio.open(src_path) as src: # 确保矢量坐标系和栅格一致不一致先转 if zhengzhou.crs ! src.crs: zhengzhou zhengzhou.to_crs(src.crs) # mask 返回裁剪后的数组和变换参数 out_image, out_transform mask(src, zhengzhou.geometry, cropTrue) out_meta src.meta.copy() # 更新元数据保持地理信息正确 out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) # 写出裁剪结果 with rasterio.open(zhengzhou_lulc_2020.tif, w, **out_meta) as dest: dest.write(out_image)mask函数的第一个参数是打开的栅格数据集第二个是几何对象列表cropTrue表示裁剪后去掉外围空白减小文件体积。out_meta复制原元数据后必须更新height、width和transform否则写出的文件地理参考会错位。坐标系一致性检查不能省我遇到过矢量是 CGCS2000 而栅格是 WGS84 UTM 的情况不转直接裁结果边界偏了几百米。如果裁剪多个县建议写成循环每个县输出一个文件命名带上行政区代码方便后续批量统计。3.3 面积统计从像元数到平方公里别忘了几何校正裁剪完就可以统计各类别面积了。核心逻辑是统计每个类别码的像元数乘以单个像元面积再换算成平方公里。如果数据是投影坐标系且分辨率严格 10m单个像元面积就是 100 平方米。# 读取裁剪后的栅格 with rasterio.open(zhengzhou_lulc_2020.tif) as src: data src.read(1) nodata src.nodata # 排除无效值 if nodata is not None: valid data[data ! nodata] else: valid data # 统计每个类别的像元数 classes, counts np.unique(valid, return_countsTrue) # 单个像元面积投影坐标系下 10m 分辨率即 100 平方米 pixel_area 100 # 平方米 # 换算成平方公里 for cls, cnt in zip(classes, counts): area_km2 cnt * pixel_area / 1e6 print(f类别 {cls}: {cnt} 像元, 面积 {area_km2:.2f} 平方公里)np.unique配合return_countsTrue直接返回类别和对应像元数比循环快得多。pixel_area这里写死 100 是因为确认了投影坐标系下分辨率为 10m如果数据是地理坐标系这个值会随纬度变化必须先用gdalwarp重投影到等面积投影或 UTM再统计。1e6是平方米到平方公里的换算系数。如果要做多年对比建议把结果存成 CSV字段包括行政区、年份、类别、面积后续直接读表做变化矩阵。我一般还会加一列百分比方便快速看结构。4. 避坑与排查10m 省级土地覆盖数据最常见的五个翻车点4.1 现象统计出来总面积比河南实际面积大很多原因通常有两个一是没排除 NoData背景值被算进了某一类二是数据分幅之间有重叠合并时没有去重。解决方法是先确认src.nodata如果为None用直方图看哪个值出现频率异常高手动设为无效值分幅数据用gdal_merge.py合并时加-n 0指定无效值或者用gdalbuildvrt先建虚拟镶嵌再转 tif避免重复像元被累加。4.2 现象裁剪后的栅格边界和矢量对不上整体偏移原因几乎都是坐标系不一致。栅格是 UTM 49N矢量是 WGS84 地理坐标直接裁就会偏。解决方法是统一坐标系用to_crs把矢量转到栅格坐标系或者用gdalwarp把栅格转到矢量坐标系。注意不要反复转转一次就固定下来反复转会引入重采样误差。4.3 现象类别码和说明文档对不上多出几个没见过的值原因可能是数据经过了二次处理或者说明文档对应的是另一版产品。解决方法是先做全图唯一值统计把实际出现的类别列出来再对照公开文档反推。如果多出的值集中在边界区域可能是混合像元或边缘效应可以考虑用众数滤波平滑但不要过度处理否则会损失真实细节。4.4 现象整省数据读进内存直接爆掉原因是一次性read了十几亿像元。解决方法是分块读取用rasterio的block_windows或者window参数每次只读一块统计完累加。也可以先用gdal_translate做降采样生成低分辨率概览用于快速预览确认无误后再对目标区域做全分辨率处理。4.5 现象面积统计结果和官方公布数据差几个百分点原因可能是分辨率导致的边界混合像元、类别定义差异、或者统计时包含了水域中的岛屿等细节。解决方法是对比时明确口径比如官方耕地面积是否包含临时种植园建设用地是否包含农村道路。如果差异在 5% 以内通常属于正常范围如果超过 10%要检查投影和无效值处理。我一般会保留一份统计脚本和中间结果方便回溯每一步。5. 进阶技巧用分块统计加变化矩阵把 10m 数据用出 30m 没有的细节5.1 分块统计内存不够时的标准操作整省 10m 数据在普通笔记本上确实吃力但分块统计可以解决。核心思路是按行块或窗口遍历每块统计类别像元数最后累加。from rasterio.windows import Window def block_statistics(tif_path, block_size2048): 分块统计各类别像元数适合大文件 counts {} with rasterio.open(tif_path) as src: nodata src.nodata height, width src.height, src.width # 按块遍历block_size 控制每次读取的行数 for row in range(0, height, block_size): for col in range(0, width, block_size): # 计算当前窗口防止越界 win Window(col, row, min(block_size, width - col), min(block_size, height - row)) data src.read(1, windowwin) if nodata is not None: data data[data ! nodata] classes, cnts np.unique(data, return_countsTrue) for c, n in zip(classes, cnts): counts[c] counts.get(c, 0) int(n) return counts # 调用并换算面积 result block_statistics(henan_lulc_2020/henan_2020_lulc.tif) for cls, cnt in sorted(result.items()): print(f类别 {cls}: {cnt * 100 / 1e6:.2f} 平方公里)block_size默认 2048意味着每次读 2048 行内存占用可控。Window的四个参数是列起点、行起点、宽度、高度用min防止最后一块越界。统计结果累加到字典里最后统一换算。这个方法比一次性读整图慢一点但内存占用从几个 GB 降到几十 MB普通办公机也能跑。如果要做多年变化矩阵把每年的统计结果存成字典再两两对比就能得到耕地转建设用地、林地转耕地等转移面积。5.2 变化矩阵两期数据叠加看谁在变有了分块统计的基础做变化矩阵就是多一步叠加。假设你有 2020 和 2015 两期同分辨率数据先确保坐标系和范围一致然后逐像元比较。# 假设两期数据已经对齐读取同一窗口 with rasterio.open(henan_2015_lulc.tif) as src15, \ rasterio.open(henan_2020_lulc.tif) as src20: # 读取同一区域这里以整图为例实际建议分块 old src15.read(1) new src20.read(1) # 排除无效值 mask_valid (old ! src15.nodata) (new ! src20.nodata) old_valid old[mask_valid] new_valid new[mask_valid] # 组合成转移对统计频次 transition np.zeros((10, 10), dtypenp.int64) for o, n in zip(old_valid, new_valid): transition[o, n] 1 # 打印主要转移方向 for i in range(10): for j in range(10): if transition[i, j] 0 and i ! j: print(f{i} - {j}: {transition[i, j] * 100 / 1e6:.2f} 平方公里)transition矩阵的行是旧类别列是新类别对角线是未变化部分。mask_valid确保两期都有效的像元才参与统计。循环部分如果数据量大可以用numpy的bincount加速但为了逻辑清晰这里用双重循环展示。实际跑的时候建议先分块再叠加否则两期整图同时读进内存容易爆。变化矩阵最有价值的地方是能看出「耕地转建设用地」这种关键转移10m 分辨率下连农村宅基地的扩张都能捕捉到这是 30m 数据做不到的。5.3 一个我常犯的错忘了记录处理链最后说个血泪经验。我早期处理这类数据时裁剪、重投影、统计各写一个脚本跑完就扔结果三个月后要复现某个市的面积发现忘了当时用的哪版矢量边界也忘了有没有做众数滤波。后来我养成了一个习惯每个项目目录下放一个README.md记录数据来源、处理步骤、关键参数、每一步的输出文件名。哪怕只是临时分析也写三行。这个习惯帮我省了无数次后悔药。10m 省级数据量不小处理一次不容易把处理链记清楚比多跑几次脚本重要得多。希望帮到你。本文还有配套的精品资源点击获取