简介：湖南省土地利用覆盖数据集是一套覆盖2000年、2005年、2010年、2015年和2020年五个关键时相的TIF格式地理空间数据，面向GIS、遥感与土地科学研究人员，主要解决土地利用变化分析与趋势研判问题。资源包共30个文件，含五个年份的遥感分类图层，并配备DBF地类属性表、XML元数据、OVR低分辨率快视图和TFW坐标配准文件，压缩包大小约157.3MB，整体结构清晰，可直接导入ArcGIS等平台开展空间分析。已有537人学习下载。借助这些图层，可对比不同年份耕地、林地、水域和建设用地等地类的分布变化，识别近二十年土地转换规律；也可结合社会经济数据做可持续性评价，或进行缓冲区、热点分析等GIS操作，为国土规划与生态保护提供参考。

1. 湖南省土地利用覆盖数据集：五期tif不是直接出图的，是要拿来算的

拿到《湖南省土地利用覆盖数据集（2000-2005-2010-2015-2020）五个年份数据集tif》这份数据的人，多数是冲着做长株潭扩张、洞庭湖湿地变化或者耕地红线评估去的。我最早用这类数据时也以为解压后拖进ArcMap就能出五张专题图，结果发现五个tif的分辨率、坐标系、地类编码并不完全一致，直接叠加做变化检测，出来的全是假变化。这个数据集真正值钱的不是那五张静态图，而是你把它处理成一套口径统一的时间序列后，能算出二十年里每一块地的去向——这才是土地利用数据集该有的用法。适合谁？做国土空间规划、生态评价、地理国情监测的从业者，以及需要拿现成分类结果做分析而不是自己跑分类算法的人。下文按“先核对数据、再裁剪、再算转移矩阵、最后排坑”的顺序展开，照着做基本能出可用成果。

2. 动手前先核对：五个年份tif的坐标系、分辨率与地类编码表

2.1 五个年份tif的数据构成：先分清是单景整幅还是分幅拼接

常见做法是，这类省级土地利用数据集解压后有两种形态：一种是一个完整的湖南省范围tif，直接覆盖全省；另一种是按标准分幅（如1:10万图幅）拆成多个tif，需要先拼接。我一般不会急着拼，先看文件列表里有没有“湖南省_LUCC_2000.tif”这类整幅文件。如果是分幅的，才考虑用栅格目录统一管理，或者用镶嵌到新栅格生成一个虚拟拼接层。

不管哪种形态，第一步都是用GDAL把五个tif的元数据全部列出来，确认它们是不是已经对齐。实际项目中最常见的坑，是2000年和2005年用的是Xian 80 / GK（高斯-克吕格）投影，2010年以后换成了CGCS2000，甚至有的年份直接给了WGS84经纬度。坐标系不一致，后面所有像元操作都是错的。

for f in hunan_2000.tif hunan_2005.tif hunan_2010.tif hunan_2015.tif hunan_2020.tif; do echo "== $f ==" gdalinfo "$f" | grep -E "Origin|Pixel Size|Coordinate System|Type=" done

这段循环把五个tif的起点坐标、像元大小、投影和数据类型一次列齐。重点看三处：一是Coordinate System是否完全一致，二是Pixel Size是否都是30米，三是GDAL Type是否相同。Data Type如果是Int16或者Byte，属性表才可能带地类代码；如果是Float32，多半是别人做过处理或重采样，属性表可能已经丢了。

2.2 属性表与地类编码：读不出代码，一切分析无从谈起

土地利用覆盖tif的每个像元存的是地类编码，不是RGB颜色。像元值1通常代表耕地、2代表林地、3代表草地，但不同批次的编码未必一样，有的用6类一级类，有的用25类二级类。拿到数据后没有项目文档时，我的习惯是先读属性表，没有属性表就统计唯一值和像元数量。

import rasterio from rasterio.plot import show import numpy as np for year in [2000, 2005, 2010, 2015, 2020]: path = f"hunan_{year}.tif" with rasterio.open(path) as src: band = src.read(1) vals, counts = np.unique(band, return_counts=True) nodata = src.nodata print(f"{year}: nodata={nodata}, 像元类别数={len(vals)}") for v, c in zip(vals[:12], counts[:12]): print(f" 代码 {v}: {c} 个像元")

这里把五个年份的像元值分布全部打出来，目的是确认年份之间的地类编码是否有漂移。如果2000年有代码0作为背景，而2020年用255作为NoData，不统一编码就去算转移矩阵，会出现整片“无中生有”。如果某个年份只有两三种像元值，那说明tif可能被拉伸成了RGB显示图，而不是原始分类栅格，这种文件不能直接用于计算，只能当底图看。

2.3 用地类代码重映射：统一到一套口径再动手

检查完编码后，下一步就是把五个tif统一到同一套分类体系。我这里给出一个最少操作：如果原始代码是二级类，先归并成一级类；如果年份之间有0和255的NoData差异，先把NoData统一改成255，再丢进后续流程。

import rasterio import numpy as np src_to_l1 = {11:1, 12:1, 21:2, 22:2, 23:2, 31:3, 32:3, 41:4, 42:4, 51:5, 61:6} def normalize(path_in, path_out, nodata=255): with rasterio.open(path_in) as src: data = src.read(1).astype("int16") profile = src.profile.copy() profile.update(dtype="int16", nodata=nodata, count=1) out = np.full(data.shape, nodata, dtype="int16") for k, v in src_to_l1.items(): out[data == k] = v with rasterio.open(path_out, "w", **profile) as dst: dst.write(out, 1) normalize("hunan_2000.tif", "hunan_2000_l1.tif")

这段代码的逻辑是：逐像元把二级地类编码映射为一级编码，没有映射到的像元一律置为NoData。参数说明：src_to_l1中的11、12代表水田、旱地归为耕地（码1），21-23归为林地（码2），31、32归为草地（码3），41、42归为水域（码4），51归为建设用地（码5），61归为未利用地（码6）。如果你的数据是6类直接编码，这步可以跳过；这一步是后面做转移矩阵和面积统计的前提，不做的话五个年份的图例会对不上。

3. 从全省tif到研究区成果：裁剪、掩膜提取与批量处理的三种做法

3.1 ArcMap或QGIS里的面图层裁剪：常规Extract by Mask与Con的差别

热词里“arcmap中依靠面图层裁剪DEM栅格”指的就是按行政边界切栅格。在ArcMap中，工具箱里有两个外观相近的工具：按掩膜提取（Extract by Mask）和按矩形裁剪（Clip）。按掩膜提取会严格以面要素的边界生成不规则裁剪结果，边界外是NoData；栅格裁剪则输出面边界的外接矩形，边界外保留背景值。做地类面积统计时，我坚持用按掩膜提取，因为外接矩形会把省外像元算进去，统计面积直接虚高。

QGIS里对应的是“栅格 → 裁剪 → 按掩膜图层裁剪”。有一个细节必须注意：面图层和tif的坐标系必须一致，否则裁剪结果会产生轻微偏移。我一般会先在面图层上“导出要素（重投影）”到tif相同坐标系，再执行裁剪。

3.2 用rasterio批量裁剪五期tif：一套代码出五张成果图

如果只是裁一次，ArcMap点选就行；但五期tif都按同一个研究区裁剪，手工点选五次容易因为捕捉或字段选择不一致而出错。我更喜欢用rasterio写批量脚本，保证五期结果严格对齐。下面是按矢量边界裁剪并重采样的完整流程。

import rasterio import rasterio.mask import geopandas as gpd import numpy as np boundary = gpd.read_file("changzhutan.shp").to_crs("EPSG:32649") for year in [2000, 2005, 2010, 2015, 2020]: with rasterio.open(f"hunan_{year}_l1.tif") as src: geom = [boundary.geometry.unary_union] out_image, out_transform = rasterio.mask.mask( src, geom, crop=True, nodata=255 ) profile = src.profile.copy() profile.update( height=out_image.shape[1], width=out_image.shape[2], transform=out_transform, nodata=255, ) with rasterio.open(f"czt_{year}_l1.tif", "w", **profile) as dst: dst.write(out_image) print(f"{year} 裁剪完成，尺寸 {out_image.shape}")

参数说明：boundary先重投影到tif坐标系，避免两个图层基准不一致；rasterio.mask.mask的crop=True表示按边界最小外接矩形裁剪并压缩尺寸，nodata=255把边界外全部填成背景。输出的五张tif形状完全一致，像元一一对应，这是后一步做转移矩阵的前提。如果五个tif的分辨率不一致，裁剪时还要加resampling参数统一重采样为30米，否则后面的面积统计会失真。

3.3 裁剪后必须做的质量检查：边界黑边、错位和面积偏小

裁剪完别急着算面积，先做三件小事。第一，把裁剪结果叠加到影像或在线底图上，确认边界没有“黑边”——出现黑边的原因是面图层边界与tif边缘之间存在NoData条带，一般由坐标系不一致导致。第二，用“栅格唯一值”工具统计裁剪结果的像元数量，与原始tif按边界粗算的数量对比，偏差超过5%就要检查投影。第三，把五期裁剪结果的像元尺寸用gdalinfo再列一次，确保完全一致。

4. 算二十年地怎么变：转移矩阵、面积统计与变化检测的落地实现

4.1 数据准备好了，先解决“两期栅格怎么对到一起”

土地利用转移矩阵的本质是t时刻的类别i在t+1时刻变成类别j的面积统计。在Arcgis中有“栅格叠置分析→交叉制表”工具可以一键输出转移矩阵，但前提是两期tif已经像元对齐，且类别编码一致。如果你是按３.２节流程做的裁剪，这个前提已经成立。如果直接用原始五个tif去做，20年前的栅格和现在的栅格像元起点可能差半个像元，算出来的转移矩阵几乎全是噪声。

判断两期栅格是否对齐的快速方法是：用rasterio打开两张tif，读transform，比较左上角坐标和像元尺寸。

import rasterio def check_aligned(path_a, path_b): with rasterio.open(path_a) as a, rasterio.open(path_b) as b: print("A transform:", a.transform) print("B transform:", b.transform) print("形状一致:", a.shape == b.shape) print("仿射一致:", a.transform == b.transform) check_aligned("czt_2000_l1.tif", "czt_2020_l1.tif")

对齐检查完，如果输出全是False，说明不能直接做转移矩阵，需要重投影并重采样。使用rasterio.warp.reproject把晚年份数据重采样到早年份的网格上，这里强调一下重采样方法选择：地类栅格是离散分类数据，只能选nearest最近邻，不能选bilinear或cubic，否则地类边界会混出很多不存在的类别。

4.2 用numpy直接算转移矩阵：摆脱ArcGIS的依赖

如果你的ArcGIS许可在机构外不方便用，或者只是要快速验证某两期变化，可以用numpy一行统计转移矩阵。这个思路适合任何能读进Python的分类栅格，速度也快。

import numpy as np import rasterio def transition_matrix(path_a, path_b, n_classes=6, nodata=255): with rasterio.open(path_a) as a, rasterio.open(path_b) as b: da = a.read(1).astype("int16") db = b.read(1).astype("int16") valid_a = da != nodata valid_b = db != nodata valid = valid_a & valid_b matrix = np.zeros((n_classes, n_classes), dtype="int64") np.add.at(matrix, (da[valid], db[valid]), 1) return matrix mat_00_20 = transition_matrix("czt_2000_l1.tif", "czt_2020_l1.tif") print("2000→2020 转移矩阵（单位：像元）") print(mat_00_20)

逻辑说明：这个函数把两期栅格里同为有效值的像元找出来，以2000年的类为行、2020年的类为列，统计每个“来源类→目标类”组合的像元数量。n_classes=6对应前面统一后的一级类；如果保留二级类，这里要改成实际类别数。像元数乘以单像元面积（30米×30米=0.09公顷）就是转移面积，单位公顷。matrix[i][j]表示年份A的类i变成了年份B的类j的数量。

4.3 变化检测别只看总量：一张图找出“哪变了”

转移矩阵只能告诉你数量，不能告诉变化发生在哪。更常见的工作流是把两期tif按公式“变化后 = 旧值×100 + 新值”合成一个双时相变化编码栅格，再对特定组合赋色。这样做可以把“耕地转为建设用地”这一类变化单独提出来出图。

import numpy as np import rasterio with rasterio.open("czt_2000_l1.tif") as a, rasterio.open("czt_2020_l1.tif") as b: da = a.read(1) db = b.read(1) profile = a.profile.copy() change_code = np.where((da != 255) & (db != 255), da * 100 + db, 255) change_farm_to_built = np.where(change_code == 105, 1, 0) with rasterio.open("czt_change_2000_2020.tif", "w", **profile) as dst: dst.update(nodata=255, dtype="int16", count=1) dst.write(change_code, 1)

参数说明：da*100+db生成的编码中，105代表耕地（码1）变成建设用地（码5），205代表林地变成建设用地，依此类推。这样不用查转移矩阵就能快速定位哪些区域发生了“耕地流失”。后面如果要出图，再对change_code写一个颜色映射表。这个做法的好处是单波段栅格可以直接拖进ArcGIS做渲染，不需要连数据库。

5. 避坑指南：五个年份tif翻车现场与排查方法

5.1 现象：tif拖进ArcMap是黑色的，属性表也读不出地类代码

原因：五期tif里混入了经过拉伸或渲染的RGB预览图，原始分类波段被压缩成了三波段RGB，像元值不是地类编码而是一个颜色值。解决：检查Data Type是否为Byte且波段数为3，如果是，说明不是分类栅格，需要回溯原始文件。群里的数据如果只有这一版，只能退而求其次，用颜色映射表反推地类，但精度看运气，建议不要用于面积统计，当底图用即可。

5.2 现象：用面图层裁剪后，研究区边的耕地被裁掉一圈，面积比统计公报少一截

原因：面边界和栅格像元之间的配准误差，加上矢量边界比真实权属边界精度高，导致边界像元被判为NoData。解决：先做缓冲，对面图层做负缓冲拆掉边界毛刺，或者用“按掩膜提取”后再做一次“多数滤波”补边界像元。裁剪后如果非要逐级对比，必须说明统计口径是以像元归属为准，与统计公报的口径本身就存在差异。

5.3 现象：做转移矩阵时，矩阵对角线特别大，但总有一部分“1→1”明显不合理，相邻像元全变成同类

原因：两期tif的像元错位半个到几个像元，变化检测时出现系统性误配。解决：回到第4章的check_aligned检查transform，用rasterio.warp.reproject把晚一期重采样到早一期网格，重采样方法必须用nearest。如果重采样后仍有大量“1→1”伪变化，再用3×3众数滤波压掉孤立变化像元。

5.4 现象：五个年份的面积加起来对不上，2010年总像元数少了一万多个

原因：某个tif在传输或压缩过程中丢过NoData，也可能原始分类把云区和阴影直接标成了0，0被当成有效地类统计。解决：在统一编码时把所有不在类别字典里的值（0、255、NaN）全部映射为NoData，统计面积时排除NoData，并在成果表中给出每个年份的“有效面积”字段，方便后续对比。

5.5 现象：QGIS里打开五期tif，颜色显示对不上，同一块林地在不同年份颜色不同

原因：QGIS默认按像元值的直方图拉伸渲染，五期直方图不同，渲染出来颜色当然不一致。解决：在图层样式里手动设置唯一值渲染（Paletted/Unique Values），把六类地类固定成同一套颜色，比如林地绿色、建设用地红色。保存成QML样式文件后批量应用到五期tif，保证制图口径统一，这个坑在给项目做汇报图时特别常见。

6. 进阶：把五期tif合成时间序列数据集，再用变化检测验证成果

6.1 合成多波段时间序列GeoTIFF：一张tif里放五期数据

做完整套流程后，我会把五个年份的一级类栅格合并成一个五波段tif，后续做趋势分析、训练分类模型、做时序分割都从这个文件读，不用每次开五个文件。合并时波段顺序按年份排列，波段名写进tif的tags。

import rasterio import numpy as np years = [2000, 2005, 2010, 2015, 2020] paths = [f"czt_{y}_l1.tif" for y in years] profile = None stack = [] for p in paths: with rasterio.open(p) as src: if profile is None: profile = src.profile.copy() profile.update(count=len(years), dtype="int16", nodata=255) stack.append(src.read(1).astype("int16")) with rasterio.open("czt_lucc_2000_2020_stack.tif", "w", **profile) as dst: dst.write(np.stack(stack), [1, 2, 3, 4, 5]) dst.update_tags(ns="LUCC", years="2000,2005,2010,2015,2020")

这段代码把五期像元严格叠加成五波段文件，波段索引1到5对应2000到2020。注意profile.update(count=5)必须在第一个波段读取后执行，否则写入时波段数不对。叠加后的tif可以用rasterio直接读某一个波段，也可以传给xarray做时间维度分析，后续跑随机森林或者时序分割都从这一个文件读数据，比维护五个路径省心得多。

6.2 用像元级变化频率验证成果：一个简单的合理性检查

合成文件后，有必要做一次合理性验证，防止前面每一步掩盖了错误。方法：统计每个像元在这五期里的类别变化次数。耕地、林地这类稳定地类变化次数应为0或1，建设用地一旦从耕地转来后不应再变成林地，变化次数超过2的像元要重点检查。

import rasterio import numpy as np with rasterio.open("czt_lucc_2000_2020_stack.tif") as src: data = src.read() data[data == 255] = -1 change_counts = np.zeros(data.shape[1:], dtype="int16") for i in range(1, data.shape[0]): change_counts += (data[i] != data[i - 1]) & (data[i] != -1) & (data[i - 1] != -1) print("变化0次:", np.sum(change_counts == 0)) print("变化1次:", np.sum(change_counts == 1)) print("变化≥2次:", np.sum(change_counts >= 2))

逻辑说明：逐波段比较前后两期，累计每个像元的类别切换次数。变化次数过多，通常意味着归一化编码环节出了问题，常见原因是2005年和2010年的地类编码体系不一致，导致大量“假变化”。我自己的习惯是，任何变化检测成果发布前都要跑一遍这个频率统计，如果“变化≥2次”的像元超过总面积5%，先回头查编码表，而不是直接调颜色出图。

说一个个人习惯：做这类长时序土地利用数据，我会把每一步加工后的中间结果保留下来，比如原始tif、统一编码后的tif、裁剪后的tif、合成stack，各存一份。表面看占了三倍磁盘空间，但实际上任何环节出问题都能回退重做，不用从头解压原始数据。这算是吃了几次“成果做完发现2005年编码漂移”的亏之后养成的习惯。希望帮到你。

本文还有配套的精品资源，点击获取