简介：这份资源是中国1km土地利用遥感监测数据集，覆盖1980至2015年间共45个年份，面向地理学、环境科学、城乡规划及经济学等领域的研究者与政策分析人员，用于探究近四十年土地利用变迁、城市扩张、生态保护与农业动态等议题。压缩包共126个文件，以adf、nit、dat等ArcInfo栅格与属性格式为主，辅以log日志、dir目录、001分块及doc说明文档，整体约21.83MB，按年份分文件组织，便于逐年对比分析。目前已有2896人学习下载，具备一定参考热度。数据以1km×1km格网记录耕地、林地、草地、建设用地、水域等类型，可用于识别显著变化区域、评估生态政策成效、研究耕地分布调整，并与其他社会经济数据融合构建土地利用变化模型，为城市规划、灾害评估与气候变化适应研究提供科学依据。
1. 从一份 1km 土地利用数据说起：1980-2015 年三十六年变化怎么用起来
做国土空间规划、生态评估或者碳汇核算的同行，大概率都遇到过同一个尴尬：想拉一份长时序的土地利用数据，要么年份断档，要么分辨率粗到县级都勉强，要么格式乱得没法直接进 GIS。这份「中国1km土地利用数据1980-2015.zip」解决的正是这个痛点——它把 1980 到 2015 年共 36 年的土地利用栅格按 1km 网格铺开，覆盖全国，年份连续，是少见的能直接拿来做长时序对比的基础底图。土地利用数据本身不稀奇，稀奇的是「连续 36 年 + 全国 1km」这个组合，它让你能在同一套分类体系下看耕地、林地、草地、水域、建设用地和未利用地的此消彼长。适合谁？做 LUCC 变化检测的、算生态系统服务价值的、给水文模型准备下垫面输入的，以及需要一份可靠底图来交叉验证自己遥感分类结果的人。下面我按「拿到手怎么用、参数怎么设、哪里会翻车」的顺序拆一遍。
2. 数据拿到手先做什么：解压、目录结构与坐标系统核对
2.1 解压后的目录长什么样，先别急着往 GIS 里拖
压缩包解开后，常见做法是按年份分文件夹，每个文件夹里放当年的一幅全国栅格，命名里带年份和分辨率标识。我一般先不打开 GIS，而是用命令行把目录树和文件大小扫一遍，确认年份有没有缺、单幅文件是不是完整。这一步花不了两分钟，但能避免后面跑到一半发现某年文件是 0 字节的翻车。
# 列出所有年份目录，确认 1980-2015 是否连续 ls -d */ | sort # 统计每个年份目录下的文件数量和总大小 for d in */; do echo -n "$d : " find "$d" -type f | wc -l du -sh "$d" done逻辑说明：第一段确认年份目录连续性，36 年应该是 36 个目录（或按十年分组）；第二段逐个目录统计文件数和体积，正常情况每年应该只有一幅全国栅格，如果某年出现多幅，要么是分幅存储，要么是混进了辅助文件。参数上注意du -sh给的是磁盘占用，和实际像素数据量有差异，只用来做横向对比。
2.2 坐标系统和 NoData 值必须先确认，否则面积统计全是错的
土地利用数据最常见的坑就是坐标系。国内这类产品多数用 Albers 等积投影或经纬度投影，两者算出来的面积能差出百分之几。用gdalinfo看一眼元数据，重点盯四个字段：投影、行列数、像素大小、NoData 值。
# 查看单幅栅格的元数据 gdalinfo 2015/landuse_2015.tif | head -40 # 批量导出所有年份的投影和像素大小，方便对比 for f in */*.tif; do echo "=== $f ===" gdalinfo "$f" | grep -E "Pixel Size|PROJCRS|NoData" done逻辑说明：gdalinfo输出的Pixel Size如果是(0.008333, 0.008333)这类小数，说明是经纬度投影，1km 只是赤道附近的近似；如果是(1000, 1000)且单位是米，才是真正的等积投影。NoData 值常见是 0 或 -9999，但有些产品把 0 当作「未利用地」的有效值，这时候如果你按 0 去裁边，会把真实的未利用地整片抹掉。参数上，head -40只是防止输出太长，实际核对时要把PROJCRS整段看完。
提示：如果发现各年份投影不一致，别急着统一重投影，先确认原始产品是不是本来就分了两套投影体系，盲目重采样会引入额外误差。
2.3 分类体系对照表，别拿自己的编码去套
土地利用数据一般用一级类或二级类编码，不同产品编码规则不同。拿到数据后第一件事是找到配套的分类说明，把编码和地类名对起来。下面这张表是我按常见一级类整理的对照框架，具体编码以你手里那份的说明为准。
| 编码 | 一级地类 | 常见二级类举例 | 面积统计注意 |
|---|---|---|---|
| 1 | 耕地 | 水田、旱地 | 注意区分水田旱地，影响蒸散估算 |
| 2 | 林地 | 有林地、灌木林 | 疏林地容易被误判为草地 |
| 3 | 草地 | 高覆盖、中覆盖、低覆盖 | 覆盖度阈值影响碳储量 |
| 4 | 水域 | 河渠、湖泊、水库 | 年际波动大，需单独核对 |
| 5 | 建设用地 | 城镇、农村、工矿 | 增长最快，变化检测重点 |
| 6 | 未利用地 | 沙地、裸地、盐碱地 | 与 NoData 值易混淆 |
逻辑说明：这张表的作用是让你在写重分类脚本前先把映射关系定死。参数上，一级类做长时序对比最稳，二级类适合做特定地类的精细分析，但二级类年际一致性往往不如一级类，跨年份比较时建议先聚合到一级类。
3. 用 Python 做长时序变化统计：从单幅读取到 36 年趋势
3.1 用 rasterio 批量读取并统计各地类面积
单幅读取谁都会，难的是 36 年批量跑还不爆内存。我一般用 rasterio 分块读取，配合 numpy 做像素计数，最后换算成面积。下面这段脚本是核心骨架。
import rasterio import numpy as np import pandas as pd from pathlib import Path # 地类编码到名称的映射，按你手里数据的说明改 CLASS_MAP = {1: "耕地", 2: "林地", 3: "草地", 4: "水域", 5: "建设用地", 6: "未利用地"} def count_by_class(tif_path, nodata=-9999): with rasterio.open(tif_path) as src: arr = src.read(1) # 把 NoData 屏蔽掉，避免混进统计 valid = arr != nodata # 像素面积：投影单位是米时直接用，经纬度需另算 pixel_area = abs(src.transform.a * src.transform.e) result = {} for code, name in CLASS_MAP.items(): cnt = np.sum((arr == code) & valid) result[name] = cnt * pixel_area / 1e6 # 转平方公里 return result rows = [] for tif in sorted(Path(".").glob("*/*.tif")): year = int(tif.parent.name) stats = count_by_class(tif) stats["year"] = year rows.append(stats) df = pd.DataFrame(rows).set_index("year").sort_index() df.to_csv("landuse_area_1980_2015.csv") print(df.head())逻辑说明：src.transform.a和src.transform.e分别是 x 和 y 方向的像素分辨率，相乘得到单个像素的面积。如果投影是经纬度，这个值只是近似，需要按纬度做余弦校正，否则高纬度地区面积会偏大。参数上nodata要和你gdalinfo里看到的一致，1e6是把平方米转成平方公里。这段脚本跑 36 年数据，内存占用取决于单幅大小，1km 全国栅格一般单幅在几百 MB 量级，分块读取更稳，但上面这种一次性读取在 16G 内存机器上通常也扛得住。
3.2 变化检测：转移矩阵怎么算才不出错
面积统计只是第一步，真正有价值的是「哪类变成了哪类」。转移矩阵用 numpy 的bincount或者 pandas 的crosstab都能做，但两个年份的栅格必须严格对齐——行列数、投影、范围完全一致。
import rasterio import numpy as np import pandas as pd def transition_matrix(tif_a, tif_b, nodata=-9999): with rasterio.open(tif_a) as sa, rasterio.open(tif_b) as sb: a = sa.read(1) b = sb.read(1) # 对齐检查，不对齐直接报错，别硬算 assert a.shape == b.shape, "两期栅格行列数不一致" mask = (a != nodata) & (b != nodata) # 把两期编码组合成一个整数，便于统计 combo = a[mask].astype(np.int32) * 100 + b[mask] codes, counts = np.unique(combo, return_counts=True) records = [] for c, n in zip(codes, counts): from_code, to_code = divmod(c, 100) records.append({"from": from_code, "to": to_code, "pixels": n}) return pd.DataFrame(records) tm = transition_matrix("1990/landuse_1990.tif", "2015/landuse_2015.tif") tm_pivot = tm.pivot(index="from", columns="to", values="pixels").fillna(0) print(tm_pivot)逻辑说明：combo = a * 100 + b是把两期编码拼成一个唯一整数，前提是编码小于 100，如果编码超过两位需要调整倍数。divmod再拆回来。参数上assert那行是关键防线，行列数不一致时直接中断，比算出一堆错值强。转移矩阵的对角线是未变化部分，非对角线才是真正的转移量，做生态评估时重点看耕地转建设用地、林地转草地这几条线。
3.3 把结果落到地图上：出图前的符号系统
统计表出来之后，多数人还要出一张变化图。这里有个容易被忽略的点：分类配色要固定，不能每年换一套，否则读者没法对比。我一般把配色写成一个字典，出图时统一调用。
import matplotlib.pyplot as plt from matplotlib.colors import ListedColormap import rasterio COLORS = ["#FFFF00", "#00B050", "#A6D96A", "#00B0F0", "#FF0000", "#C0C0C0"] cmap = ListedColormap(COLORS) with rasterio.open("2015/landuse_2015.tif") as src: arr = src.read(1) arr = np.where(arr == -9999, np.nan, arr) plt.figure(figsize=(12, 8)) plt.imshow(arr, cmap=cmap, vmin=1, vmax=6) plt.colorbar(ticks=[1,2,3,4,5,6], label="地类编码") plt.title("2015 年土地利用") plt.savefig("landuse_2015.png", dpi=200, bbox_inches="tight")逻辑说明：ListedColormap把六个地类固定成六种颜色，vmin和vmax锁死范围，避免 matplotlib 自动归一化导致颜色漂移。参数上dpi=200是出图精度，bbox_inches="tight"去掉多余白边。注意 NoData 要先转成np.nan，否则会被当成有效值参与配色。
4. 避坑与排查：这份数据最容易翻车的五个地方
4.1 现象：面积统计出来逐年跳变，某年突然翻倍
原因：多半是某年的 NoData 值和其他年份不一致，或者该年文件本身是分幅存储，你只读了一幅。解决：先用gdalinfo把 36 年的 NoData 值全部导出对比，再把每年的文件数核对一遍，确认没有漏读分幅。
4.2 现象：转移矩阵对角线数值异常小，看起来像全变了
原因：两期栅格的投影或范围没对齐，像素错位导致同一位置被算成不同地类。解决：在transition_matrix里加assert检查 shape，再用gdalinfo核对两期的PROJCRS和Origin，不一致就先重投影对齐。
4.3 现象：建设用地面积明显偏小，和统计年鉴对不上
原因：1km 分辨率下，农村宅基地和零星工矿用地被混合像元稀释，分类时容易被归到耕地。解决：这是数据本身的尺度限制，不是 bug。做建设用地分析时要么接受这个偏差，要么用更高分辨率数据做局部校正，别硬调分类阈值去凑。
4.4 现象：经纬度投影下算出的面积比实际大
原因：经纬度投影的像素面积随纬度变化，直接用transform.a * transform.e在高纬度会偏大。解决：按行做余弦校正，或者先重投影到 Albers 等积投影再统计。我一般选后者，一步到位。
4.5 现象：出图时颜色和分类对不上，水域变成了红色
原因：matplotlib 的ListedColormap按vmin到vmax均匀分配颜色，如果实际编码不是连续的 1 到 6，颜色就会错位。解决：出图前先np.unique看实际编码，确保COLORS列表长度和编码数量一致，必要时用BoundaryNorm显式指定边界。
5. 进阶技巧：把 36 年数据压成一条变化轨迹
前面都是按年份独立处理，真正能出洞察的是把 36 年串成每个像素的变化轨迹。思路是把所有年份堆叠成一个三维数组，然后沿时间轴做编码压缩。下面这段是我常用的做法，把每个像素 36 年的地类序列转成一个字符串，再统计高频轨迹。
import rasterio import numpy as np import pandas as pd from pathlib import Path tifs = sorted(Path(".").glob("*/*.tif")) stack = [] for tif in tifs: with rasterio.open(tif) as src: stack.append(src.read(1)) stack = np.stack(stack, axis=0) # 形状: (36, rows, cols) # 把每个像素的时间序列转成字符串轨迹 rows, cols = stack.shape[1], stack.shape[2] flat = stack.reshape(stack.shape[0], -1) traj = ["".join(map(str, flat[:, i])) for i in range(flat.shape[1])] # 统计最高频的 10 条轨迹 from collections import Counter top = Counter(traj).most_common(10) for t, c in top: print(f"{t} : {c} 像素")逻辑说明：np.stack把 36 幅栅格堆成三维数组，reshape后每一列就是一个像素的完整时间序列。转成字符串是为了方便用Counter统计。参数上注意内存，36 幅 1km 全国栅格堆叠后可能到几个 GB，内存紧张就分块处理。这个方法的输出能直接告诉你「哪些像素 36 年一直是耕地」「哪些像素经历了林地到建设用地的转变」，比逐年统计更有决策价值。
注意：轨迹字符串的长度等于年份数，如果中间有缺年份，轨迹会错位，所以第一步的年份连续性检查不能省。
从那以后我每次拿到长时序栅格数据，都强制先跑一遍年份连续性、投影一致性和 NoData 值三项核对，再动任何分析脚本。这三步花十分钟，能省掉后面几小时的返工。希望帮到你。
本文还有配套的精品资源，点击获取
