简介面向土地利用变化分析与遥感应用研究压缩包内含1980、1990、1995、2000、2005、2010、2015及2020年八期中国土地利用现状遥感监测数据数据来源于资源环境科学与数据中心其中2020年数据基于2015年成果与Landsat 8遥感影像人工目视解译完成。土地利用类型划分详细包括耕地、林地、草地、水域、居民地和未利用土地六个一级类型并延伸出二十五个二级类型可直接服务于国土空间规划、生态环境评价、耕地保护及长时间序列土地利用/覆被变化研究。资源包共163个文件压缩后约26.25MB主要文件类型为adf栅格数据、dat与nit数据文件、doc说明文档及分卷rar包adf格式支持ArcGIS等主流GIS平台直接读取适合开展区域制图与变化检测分析。目前已有4499人学习下载尤其适合GIS、遥感专业的高校师生与科研人员作为基础底图或参考数据集使用可大幅节省数据预处理时间支撑多期对比研究与成果产出的快速实现。1. 拿到中国土地利用现状遥感监测数据.rar先别急着双击解压做遥感应用和土地利用研究的同行对这套数据应该都不陌生中国土地利用现状遥感监测数据是基于 Landsat 等中分辨率遥感影像结合人工目视解译与自动分类生成的栅格分类产品分辨率通常是 30 米覆盖 1990 年到 2020 年之间多个时相。数据发布方一般打包成 .rar 压缩包里面是按年份组织的 tif 栅格文件单个文件几十到几百兆不等。很多第一次用的人第一步就被卡住解压报错、tif 打开全黑、面积统计对不上、投影乱掉。这篇文章不聊虚的直接讲我从拿到压缩包到完成统计、变化检测的完整落地流程包括参数怎么设和血泪踩坑。2. 解压之前先认识数据30 米分辨率、分类体系与文件组织2.1 数据规格30 米分辨率、多期覆盖、栅格格式这套数据产品的基本规格是以 Landsat 系列影像为主要数据源经过几何纠正、辐射校正后用人工解译加自动分类的方式生成土地覆盖分类图。空间分辨率 30 米意味着一个像元对应地面 30m×30m 的范围也就是 900 平方米。对一个县城来说一张 30 米分辨率分类图的行列数通常在几千乘几千数据量并不小。时相覆盖上大多数公开版本包含 1990 年、1995 年、2000 年、2005 年、2010 年、2015 年、2020 年等年份每期都对应独立的文件。年份之间的对比就是土地利用变化分析的基础。拿到压缩包以后建议先看下载页面的说明文档确认你手里的压缩包是只有一年还是多期合集避免后面做时序分析时才发现缺期。数据格式方面常见的是 GeoTIFF 栅格理论上应该自带投影和地理参考信息。但实际下载的数据里坐标参考丢失、投影定义缺失的情况并不罕见这块我会在第四章节专门讲。还有一个容易忽略的点这套数据的分类值是字节型整数不是浮点型读取时不需要做任何辐射定标直接读像元值就行。2.2 分类体系一级类与二级类的编码规则中国土地利用现状遥感监测数据使用一套比较固定的分类系统一级类共六类二级类若干。这是做任何统计之前必须搞清楚的如果不知道像元值对应的地类后面算面积、做转移矩阵都是空谈。一级类和二级类的对应关系如下一级编码一级类名称二级编码二级类名称1耕地11 / 12水田 / 旱地2林地21 / 22 / 23 / 24有林地 / 灌木林 / 疏林地 / 其他林地3草地31 / 32 / 33高覆盖度草地 / 中覆盖度草地 / 低覆盖度草地4水域41 / 42 / 43 / 44 / 45 / 46河渠 / 湖泊 / 水库坑塘 / 永久性冰川雪地 / 滩涂 / 滩地5城乡、工矿、居民用地51 / 52 / 53城镇用地 / 农村居民点 / 其他建设用地6未利用土地61 / 62 / 63 / 64 / 65 / 66 / 67沙地 / 戈壁 / 盐碱地 / 沼泽地 / 裸土地 / 裸岩石质地 / 其他配套的还有三个扩增类型常常被单独标注包括 0 值作为背景或无效区域、以及将某些特定区域的“海洋”等以外区域标记为 0。常见的做法是先对栅格中不等于 0 的像元做统计再单独检查 0 值是否为有效的滩涂/海洋区域。我一般会在拿到数据后先把整张栅格的唯一值列表打印出来跟上面的编码表对一下。如果出现数字 7、8、9 之类的编码说明你拿到的可能是扩展分类版本或者经过重分类的成果需要另外查对应的数据字典。这种数字对不上分类表的情况通常不是数据损坏而是之前有人对原始产品做过处理。2.3 压缩包内的文件命名与目录组织一个多期压缩包解压以后常见的目录组织方式是年份平铺文件夹里直接放着 2000 年、2005 年、2010 年、2015 年、2020 年几个 tif或者以日期格式命名比如landuse_2000.tif、landuse_2005.tif。有可能还附带一个landuse_cls.txt分类说明文件或者一个README文档。这里要提醒的是不同发布渠道的目录命名规范差异很大。有的压缩包内部是一个单独文件夹解压路径直接带一层目录有的则是多个文件平铺在包根解压时如果不注意几十个 tif 全部散落在同一目录。先列出压缩包内容清单确认目录结构再整体解压能省掉后面整理文件的时间。另外我也见过发布方在 rar 里同时放了矢量边界 shp 和行政区划代码表下载后建议保存好后续按区域统计时会用到。3. 用 7-Zip 解开 rar解压命令、完整性校验与只读属性3.1 7-Zip 能解压 rar 文件吗可以而且比带广告的 rar 软件干净打开 .rar 文件最常用的解压工具是 WinRAR 和 7-Zip。WinRAR 有官方中文版和试用版没有购买授权的试用期会反复弹窗现在不少 rar 解压软件在安装时还会捆绑广告组件解压到一半弹个推荐安装的页面数据还没出来系统先多了一堆不需要的软件。7-Zip 是开源免费工具也是比较干净的选择它支持解压 rar、rar5以及 zip、7z、tar 等格式。对“7zip 可以解压 rar 文件吗”这个问题结论是可以。需要注意的是RAR5 是较新的压缩格式旧版 7-Zip 可能无法识别建议下载安装 21.00 以上版本WinRAR 5.0 以上也能处理 RAR5。如果你发现自己电脑上 7-Zip 解压 rar 报“Unsupported Method”先升级版本不要急着换工具。安装完成后建议把 7z.exe 所在目录加入系统 PATH这样能在命令行直接调用后续做批量解压和完整性校验方便很多。3.2 先列出压缩包内容再决定解压策略数据压缩包动辄几百兆甚至上 GB解压前先看清单是个好习惯。7-Zip 的命令行提供了完整的文件列表视图7z l 中国土地利用现状遥感监测数据.rar这个命令会列出压缩包内所有文件的名称、大小、压缩后大小和日期。执行后重点看三样东西第一文件总数多少个有没有明显的目录层级第二tif 文件单个多大估算解压后磁盘占用第三有没有同名的说明文档。如果列表显示某个年份的 tif 大小为 0 字节或者年份缺失说明压缩包本身不完整解压前就要返工。我这里用7z l而不是直接双击打开是因为命令行能看到被图形界面藏起来的细节比如文件数是否和下载页面标注一致。对做数据分析的人来说这一步顺手而且有效。关于文件名中的中文问题如果压缩包或者内部文件名含中文在 Windows 命令行下 7-Zip 通常能正常显示但在某些设置了非 UTF-8 代码页的系统上会乱码。遇到这种情况最简单的解决方式是先把压缩包复制到纯英文路径下再操作不要在中文字符串上纠结浪费时间。3.3 执行解压7z x 的参数说明与单文件释放确认清单无误后开始正式解压。用以下命令7z x 中国土地利用现状遥感监测数据.rar -oD:\landuse_data -y参数含义拆开说x表示按完整路径解压保留压缩包内的目录结构。如果只想把文件全部平铺到一个目录可以用e但建议用x避免不同年份同名文件互相覆盖。-oD:\landuse_data指定输出目录。注意-o后面直接跟路径中间不要加空格写成-o D:\landuse_data的话 7-Zip 会认为D是另一个参数导致解压失败或语法报错。-y表示遇到询问时全部自动确认比如目标目录不存在时自动创建或者覆盖已存在文件时不再逐个问。解压完成后检查目录里的文件数量和大小和7z l列出的结果一致才算完整。我通常在解压后马上执行下面这条校验命令防止解压过程本身被中断或文件损坏7z t 中国土地利用现状遥感监测数据.rart是测试模式会重新计算压缩包内每个文件的校验值和打包时记录的对比。如果输出里有Err或者Failed说明压缩包已经损坏解压出来的 tif 即使能打开也可能存在局部数据错误。这种损坏源文件导致的错误重解压多少次都一样需要回原站点重新下载完整压缩包。3.4 只读属性、密码保护与“强制解压”工具的坑解压出来的 tif 文件有时会带只读属性特别是在发布方用过的压缩工具比较老、或者打包时文件本身设置了只读的情况下。只读属性对直接读取影响不大但当你需要在 GIS 软件里给 tif 构建金字塔、写侧车文件.aux.xml、.ovr时就会报“无法创建”的错误。处理也简单批量去掉只读属性attrib -r -s -h D:\landuse_data\*.tif /s /d参数说明-r去掉只读-s去掉系统属性-h去掉隐藏属性/s处理所有子目录/d同时处理文件夹本身。执行完再用attrib查看确认。关于密码保护个别渠道下载的数据会带 rar 密码密码通常在下载页面的说明文字里。如果你确实把密码忘了正确做法是回到原发布页面查说明而不是去下载所谓的“密码移除”或“强制解压”工具。这类工具绝大多数是捆绑软件解压成功率极低还可能顺手装上推广程序不划算。4. 坐标参考与属性表处理把数据放到正确的投影上再统计4.1 先检查 tif 自带的空间参考信息别等统计时才发现解压完先别急着打开先用 GDAL 看一眼空间参考信息。Python 环境里如果没有 GDAL用gdalinfo命令行也可以gdalinfo D:\landuse_data\landuse2020.tif输出里重点看Driver、Size、Coordinate System is和GeoTransform这几行。如果Coordinate System is后面是空说明 tif 丢了投影定义。如果GeoTransform全为 0 或者第一、第五个参数是异常值说明地理参考也有问题。用 Python 检查更直接from osgeo import gdal ds gdal.Open(rD:\landuse_data\landuse2020.tif) print(ds.GetProjection()) print(行列数:, ds.RasterXSize, ds.RasterYSize) print(仿射变换:, ds.GetGeoTransform())GetProjection()返回空字符串时就定义投影。全国范围的土地利用数据常用的投影是 Albers 等积圆锥投影典型参数类似中央经线 105°E、两条标准纬线 25°N 和 47°N椭球体为 Krasovsky 或 CGCS2000。如果你不确定原数据投影可以靠两种方式推断一是看数据发布页的技术文档里写没写投影二是用原始影像或行政边界 shp 做参照在 GIS 里人眼比对。最忌讳的是不做任何检查就计算面积得到的数字一定是不对的。4.2 投影转换为什么要用 etc. 参数和最近邻重采样拿到数据后我一般先确定目标投影再统一转换。目标投影的选择取决于你要做什么应用场景推荐投影理由全国面积统计Albers 等积圆锥投影面积不变形是国家标准县域或市域分析对应分带的 Gauss-Kruger 投影精度高和当地测绘成果对齐与裸 Landsat 影像叠加WGS84 UTM 对应分带影像自带投影一致无需重投影出图或底图叠加与底图保持一致减少人为误差分类栅格的转换有一个铁律重采样方法必须用最近邻nearest。双线性bilinear和三次卷积cubic会为分类值插值出中间值比如耕地编码 1 和林地编码 2 之间可能生成 1.3、1.7 这种怪物不仅无法解释还会让后续统计全部报废。用 GDAL 转换的示例from osgeo import gdal gdal.Warp( rD:\landuse\landuse2020_albers.tif, rD:\landuse\landuse2020.tif, dstSRSprojaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 ellpskrass unitsm no_defs, resampleAlgnear ) print(重投影完成)参数说明dstSRS指定目标坐标系resampleAlgnear即最近邻重采样ellpskrass是 Krasovsky 椭球体参数。如果你的数据已经带投影只是和你的底图不一致直接从原坐标转到目标坐标即可。投影转换有个容易忽略的小坑转换前后统计结果可能略有差异原因在于重采样时栅格网格变化导致的像元边界重分配。要确保前后一致在转换前记录原始像元总数转换后再核对一次差异超过 0.5% 就要检查参数是不是设错了。4.3 分类编码与属性表把数字翻译成面积坐标问题解决后核心就是把编码变成统计结果。读取整张栅格并统计各分类的像元数import numpy as np from osgeo import gdal ds gdal.Open(rD:\landuse\landuse2020_albers.tif) arr ds.ReadAsArray() # 像元值是字节型直接统计出现次数 unique, counts np.unique(arr, return_countsTrue) for cls, cnt in zip(unique, counts): if cls 0: continue area_km2 cnt * 900 / 1e6 # 30m分辨率每像元900平换算平方千米 print(f分类编码 {cls}: {cnt} 像元, {area_km2:.2f} km²)逻辑说明np.unique返回唯一值和对应计数arr是二维 NumPy 数组。30 米分辨率下每个像元面积是 900 平方米除以 100 万得到平方千米。跳过 0 值是因为它表示背景或无效区域不能算作面积。这一步结果拿到后对照第 2 章的分类编码表就能直接得到各地类面积。如果你需要的是二级类面积将 11/12 分别统计即可如果想要一级类的总耕地面积就把 11 和 12 的计数相加。从这往下的常见做法是直接按一级类重分类。比如把 11、12 归为耕地21-24 归为林地31-33 归为草地41-46 归为水域51-53 归为城乡工矿居民用地61-67 归为未利用地。重分类以后再做面积统计输出更直观也更方便后续的变化检测。5. 避坑排查解压、投影、面积统计中的常见翻车现场5.1 解压到一半提示“文件末端”重解压还是一样现象用 7-Zip 解压到 60% 左右突然报错 “Unexpected end of data”重新下载后再解压仍然在同一进度附近失败。原因压缩包本身在下载过程中被截断。浏览器断点续传失效、网盘中转文件损坏都会导致压缩包数据不完整。重新下载的压缩包如果校验值未变仍然损坏。解决不要反复重试解压。先用7z t 压缩包.rar测试完整性如果报错直接回源站重新下载下载工具开启强制校验或对比发布方提供的 SHA-256 值。压缩包损坏时只有重新获取完整源文件这一条路没有其他捷径。5.2 tif 在 ArcGIS / QGIS 里打开全黑或者全绿现象在 ArcGIS 里打开 tif图层全黑一片或者用默认拉伸显示成异常的单色在 QGIS 里打开显示全绿。原因分类栅格是单波段整数数据值域在 0 到 67 之间软件默认按 RGB 三波段渲染单波段灰度拉伸的默认范围不对0 值背景占满色带看起来就是全黑。解决在 ArcGIS 的图层属性里把渲染方式改为“唯一值”Unique Values给每个分类编码配一个颜色在 QGIS 中右键图层属性渲染类型选“单波段伪彩色”或者“分类”手动给编码赋值颜色。不要拉直方图分类数据拉伸直方图没有任何意义。5.3 按经纬度坐标直接计算面积结果大得离谱现象用 ArcGIS 字段计算器对矢量化的土地利用数据直接计算面积得到的结果和统计年鉴差出好几倍。原因把未投影的经纬度坐标GCS_WGS_1984当成投影坐标系来计算面积。经度方向每度对应的地面距离随纬度变化直接算出来的“平方度”根本不是平方米数字自然不对。解决先投影到 Albers 等积圆锥投影或对应分带的高斯投影再算面积。投影操作用 ArcGIS 的 Project Raster 或者 GDAL 的 Warp 都行关键是选择等积性质好的投影。这是面积统计翻车的头号原因没有之一。5.4 打开 tif 提示“缺少空间参考”或者导出后坐标系丢失现象数据在 ArcGIS 里打开时提示 “Missing Spatial Reference”强制加载后导出成新 tif 坐标信息还是为空。原因原始 tif 的投影信息在生成或者压缩转换时被剔除。常见于某些第三方转换工具把 GeoTIFF 转成了普通 TIFF丢失了地理标签。解决先用 gdalinfo 检查原始文件是否有投影如果是空使用定义投影工具手动指定。如果发布方文档说明了坐标系就用文档参数没说就用周围地物比对推断。定义完成后再 gdalinfo 复核一次确认投影字符串非空再往下走。5.5 面积统计结果和实际明显不符怀疑数据精度现象经过正确处理和投影后某县耕地面积和统计年鉴差异仍然很大甚至出现目视明显是湖泊的水域统计面积偏小。原因数据本身的几何精度和分类精度都存在误差。30 米分辨率下小水体、小带状地物容易被错分或漏分边界类地物天然有像元化误差。另外部分数据是经过概括和化简的版本不是原始分类结果。解决把统计结果当作趋势参考而不是真值。验证方法是在原始影像上随机抽几个区域目视解译比对分类标签或者和同地区更高分辨率土地覆盖产品交叉验证。如果只是做趋势分析和变化检测这类系统误差在可接受范围内。6. 进阶用法多期土地利用数据的变化检测与面积转移矩阵6.1 两期栅格差值运算快速提取变化图斑2000 年和 2020 年两期数据在手第一步就是做差值。相同位置像元的分类编码不同说明这个像元的地类发生了变化。import numpy as np from osgeo import gdal def read_band(path): ds gdal.Open(path) return ds.ReadAsArray(), ds arr2000, ds2000 read_band(rD:\landuse\landuse2000_albers.tif) arr2020, _ read_band(rD:\landuse\landuse2020_albers.tif) # 两期栅格必须行列数一致 assert arr2000.shape arr2020.shape, 两期栅格大小不一致先统一网格 # 变化像元标记2000和2020都有有效分类且编码不同 change np.where( (arr2000 0) (arr2020 0) (arr2000 ! arr2020), 1, 0 ) # 写入变化图斑 driver gdal.GetDriverByName(GTiff) out_ds driver.Create( rD:\landuse\change_2000_2020.tif, ds2000.RasterXSize, ds2000.RasterYSize, 1, gdal.GDT_Byte ) out_ds.SetGeoTransform(ds2000.GetGeoTransform()) out_ds.SetProjection(ds2000.GetProjection()) out_band out_ds.GetRasterBand(1) out_band.WriteArray(change) out_ds.FlushCache() print(变化图斑已输出)逻辑说明np.where在三元条件成立时填 1否则填 0得到一个值只有 0 和 1 的变化掩膜。写 tif 时用原始数据的仿射变换和投影参数保证输出空间参考一致。创建的栅格是单波段字节型占用空间很小。一个重要提醒在做差值前必须先确认两期栅格的行列数和投影完全一致。不一致的话先做重投影和网格对齐再计算差值否则结果全是无效数据。6.2 用转移矩阵量化“从什么变成什么”差值能告诉你哪里变了但回答不了“耕地变成建设用地有多少、林地变成草地有多少”这类问题。转移矩阵把两期分类编码做交叉统计import numpy as np arr2000 read_band(rD:\landuse\landuse2000_albers.tif)[0] arr2020 read_band(rD:\landuse\landuse2020_albers.tif)[0] classes [1, 2, 3, 4, 5, 6] idx {v: i for i, v in enumerate(classes)} cls_idx {i: v for i, v in enumerate(classes)} # 只统计两期都落在有效一级类里的像元 mask np.isin(arr2000, classes) np.isin(arr2020, classes) a arr2000[mask] b arr2020[mask] # 用一种成对编码技巧高效统计转移量 pair_code a * 100 b unique, counts np.unique(pair_code, return_countsTrue) matrix np.zeros((len(classes), len(classes)), dtypenp.int64) for code, cnt in zip(unique, counts): from_cls code // 100 to_cls code % 100 matrix[idx[from_cls], idx[to_cls]] cnt print(转移矩阵 (行: 2000年, 列: 2020年)) print( .join(cls_idx.values())) for i, row in enumerate(matrix): print(f{cls_idx[i]} .join(f{v:8d} for v in row))逻辑说明转移矩阵的行代表 2000 年的一级类列代表 2020 年的一级类。对角线上是两期保持不变的面积非对角线就是转移量。用a * 100 b把“从类 A 到类 B”的组合压缩成一个整数再用np.unique一次统计全部分组计数避免写逐像元循环处理大型栅格时性能差距明显。拿到转移矩阵后一个常用操作是计算变化率比如耕地转建设用地的面积占耕地总面积的百分比。我一般在矩阵基础上加上行和列占比列导出成 Excel 做报告。另一个验证方式是把最大变化类型单独提取出来叠到当年原始影像上做目视抽查——变化检测说到底是为了缩小嫌疑范围最终结论还要回到影像上确认。做这套流程我有个习惯两期差值前先数一遍各自像元总数确认行列一致面积统计结果出来再和统计年鉴交叉比对一次差太多就先查投影参数再怀疑数据。这套检查看起来繁琐但每次都能拦截到至少一个低级错误倒不是技术玄学纯粹是数据工作中交过的学费。希望帮到你。本文还有配套的精品资源点击获取