简介湖南省土地利用覆盖数据集是一套覆盖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 fhunan_{year}.tif with rasterio.open(path) as src: band src.read(1) vals, counts np.unique(band, return_countsTrue) 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, nodata255): with rasterio.open(path_in) as src: data src.read(1).astype(int16) profile src.profile.copy() profile.update(dtypeint16, nodatanodata, count1) out np.full(data.shape, nodata, dtypeint16) 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代表水田、旱地归为耕地码121-23归为林地码231、32归为草地码341、42归为水域码451归为建设用地码561归为未利用地码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(fhunan_{year}_l1.tif) as src: geom [boundary.geometry.unary_union] out_image, out_transform rasterio.mask.mask( src, geom, cropTrue, nodata255 ) profile src.profile.copy() profile.update( heightout_image.shape[1], widthout_image.shape[2], transformout_transform, nodata255, ) with rasterio.open(fczt_{year}_l1.tif, w, **profile) as dst: dst.write(out_image) print(f{year} 裁剪完成尺寸 {out_image.shape})参数说明boundary先重投影到tif坐标系避免两个图层基准不一致rasterio.mask.mask的cropTrue表示按边界最小外接矩形裁剪并压缩尺寸nodata255把边界外全部填成背景。输出的五张tif形状完全一致像元一一对应这是后一步做转移矩阵的前提。如果五个tif的分辨率不一致裁剪时还要加resampling参数统一重采样为30米否则后面的面积统计会失真。3.3 裁剪后必须做的质量检查边界黑边、错位和面积偏小裁剪完别急着算面积先做三件小事。第一把裁剪结果叠加到影像或在线底图上确认边界没有“黑边”——出现黑边的原因是面图层边界与tif边缘之间存在NoData条带一般由坐标系不一致导致。第二用“栅格唯一值”工具统计裁剪结果的像元数量与原始tif按边界粗算的数量对比偏差超过5%就要检查投影。第三把五期裁剪结果的像元尺寸用gdalinfo再列一次确保完全一致。4. 算二十年地怎么变转移矩阵、面积统计与变化检测的落地实现4.1 数据准备好了先解决“两期栅格怎么对到一起”土地利用转移矩阵的本质是t时刻的类别i在t1时刻变成类别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_classes6, nodata255): 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), dtypeint64) 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_classes6对应前面统一后的一级类如果保留二级类这里要改成实际类别数。像元数乘以单像元面积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(nodata255, dtypeint16, count1) dst.write(change_code, 1)参数说明da*100db生成的编码中105代表耕地码1变成建设用地码5205代表林地变成建设用地依此类推。这样不用查转移矩阵就能快速定位哪些区域发生了“耕地流失”。后面如果要出图再对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也可能原始分类把云区和阴影直接标成了00被当成有效地类统计。解决在统一编码时把所有不在类别字典里的值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 [fczt_{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(countlen(years), dtypeint16, nodata255) 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(nsLUCC, years2000,2005,2010,2015,2020)这段代码把五期像元严格叠加成五波段文件波段索引1到5对应2000到2020。注意profile.update(count5)必须在第一个波段读取后执行否则写入时波段数不对。叠加后的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:], dtypeint16) 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年编码漂移”的亏之后养成的习惯。希望帮到你。本文还有配套的精品资源点击获取