简介：中国1km土地利用数据1980—2015压缩包面向地理学、环境科学、城乡规划与资源管理等领域的研究者和学生，提供覆盖1980至2015年多个时点的1km格网土地利用类型信息，可用于分析近四十年土地变迁、城市扩张、生态保护与农业动态。包内共126个文件，以adf、nit、dat、dir等ArcGIS栅格单元文件为主，并含log、doc、xml等说明与辅助文件，压缩包约21.83MB，文件按年份与数据集组织。已有2904人学习下载。数据以栅格形式记录耕地、林地、草地、建设用地、水域等类型，便于对比不同年份识别显著变化，也可与社会经济数据结合构建土地利用变化模型。对开展区域案例研究、政策评估和可持续分析具有较高参考价值。

1. 1km 分辨率能干什么：先把这套土地利用栅格的适用边界说清楚

一个做县域生态评估的 A 同学，拿到 1980–2015 的中国 1km 土地利用数据之后，直接按像元数乘以 1 平方公里，去算耕地转出的面积。第一次跑出县城建成区扩张 8%，换一个投影再跑变成 3%。他以为是数据坏了，其实这份土地利用数据里的问题几乎都集中在“单位”和“编码”这两件事上。这套资源是 1980 到 2015 年逐期的全国 1 公里栅格，按年份分包，解压之后就是可以直接读进 GDAL 的类别栅格，省掉了自己去拼接、配准、对齐的整段流程。它适合做长时序、大范围的地类结构分析，比如耕地流向、城镇化扩张、林地草地的年际波动，也适合给生态模型当驱动数据。但它不适合地块级权属判断，也分不清城市内部的用地功能。1 公里就是一个像元一平方公里，边界模糊、细小地物会被吞掉，这是分辨率本身决定的，跟数据质量无关。

2. 解压之后先别画图：把目录、编码表和元数据核一遍

拿到一个压缩包，最省时间的做法不是先把它塞进 QGIS 看一眼出图效果，而是花二十分钟把三件事确定下来：文件是怎么组织的、编码是什么意思、栅格自己的元数据靠不靠谱。这三件事没定下来就动手算面积，后面所有的数字都是悬空的。我把这一步叫体检，体检做完，后面写脚本才有依据。

2.1 摊开压缩包：年份目录与文件命名的识别规律

先看清单再解压，不要一口气全铺开。很多全国尺度的栅格包体量大，解压完才发现里面是 ArcInfo Grid 那种一堆.adf的目录结构，单独拷一个文件出去就坏了，得整个目录一起搬。用清单先把结构摸清，能省掉一次返工。

# 先看清单，确认里面是单文件栅格（tif/img）还是 ArcInfo Grid 目录 unzip -l "中国1km土地利用数据1980-2015.zip" | head -40 # 解压到独立目录，别和别的项目数据混在一起 unzip -q -d ./lu_raw "中国1km土地利用数据1980-2015.zip" # 列出所有栅格候选，按路径排序后年份顺序一眼可见 find ./lu_raw -type f \( -iname "*.tif" -o -iname "*.img" \) | sort

Windows 上没有unzip就用 7-Zip 的图形界面做同样的事，重点是先看清单、再解压到独立目录。清单里如果出现.prj、.tfw、.aux.xml这些小文件，说明栅格带坐标和统计附属信息，搬家的时候要跟主文件一起走，漏了.prj有些软件会默认按像素坐标处理，位置直接跑偏。

解压完把年份抽出来，建一张“年份 → 文件路径”的映射表。后面所有脚本都靠这张表循环，不要在每个脚本里重新拼文件名。

import re from pathlib import Path files = sorted(Path("./lu_raw").rglob("*.tif")) year_map = {} for f in files: m = re.search(r"(19|20)\d{2}", f.stem) # 从文件名里抓四位年份 if m: year_map[int(m.group(0))] = f print(sorted(year_map)) # 期望看到 1980、1990、1995、2000、2005、2010、2015 之类的序列

这段代码的逻辑很直白：遍历解压目录下所有.tif，用正则从文件名主干里抓年份，建成字典。如果打印出来的年份少于预期，说明命名规则不是纯年份，可能带了版本后缀或者分幅编号，这时候先把files的前十个路径打出来看规律，别硬套正则。年份数量对不上是最容易被忽略的第一个坑，少一期数据做时序分析，结论能差出好几年。

2.2 分类编码表：一级类看趋势，二级类看结构

这一类全国土地利用栅格，主流用的是两级分类：一级类 6 个，二级类二十多个，用一个整数编码表示。一级类够做趋势和总量，二级类才用来看内部结构，比如耕地里水田和旱地的比例、林地里是有林地还是灌木林。

这张表要当成“常见口径”来用，不能当成“这份数据一定如此”来用。不同分发批次里，二级类可能做过合并，也可能某些年份只给了一级类编码。我一般会先跑一遍值域统计，看实际出现了哪些整数，再和这张表对照，出现的编码是不是都在表里、有没有超出 67 的值。有超出就说明这一期用了别的口径，得单独处理。

提示：先跑值域统计再决定用哪套编码解释数据，比对着文档猜快得多，也准得多。

把二级编码折算回一级，最省事的写法是整除：

def to_level1(code: int) -> int: """两位二级编码折算成一级编码：11->1，43->4，67->6""" return code // 10 if code >= 10 else code

注意这个函数只对两位编码成立。如果某期数据里混进了三位编码，整除会算出 10 以上的值，这时候要改成查表。这也是为什么我坚持先看值域：函数写得再漂亮，口径不对一样算错。

2.3 用 gdalinfo 和 rasterio 给数据做一次体检

体检的核心是三件事：坐标系是什么、像元多大、NoData 是怎么声明的。这三项决定了后面面积能不能直接乘。

gdalinfo -stats "LU_2015.tif"

-stats会顺带算出最小最大值和直方图，输出里的Coordinate System is、Pixel Size、NoData Value三个字段要重点看。如果Pixel Size显示的是类似(0.00833, -0.00833)这种小数，说明这份数据是经纬度网格，不是米制投影——这一点非常关键，后面第 3 章会展开。

再用 rasterio 批量体检一遍，把关键属性整理成能直接比较的结构：

import rasterio import numpy as np with rasterio.open("LU_2015.tif") as src: print("CRS :", src.crs) print("宽高 :", src.width, src.height) print("像元尺寸 :", src.res) # (x, y)，单位跟随 CRS print("NoData :", src.nodata) print("数据类型 :", src.dtypes[0]) arr = src.read(1) vals, counts = np.unique(arr, return_counts=True) for v, c in zip(vals, counts): print(f"编码 {int(v):>3} 像元 {int(c):>12,}")

src.res返回的是像元在x、y方向上的尺寸，单位跟 CRS 走：米制投影就是米，经纬度就是度。src.nodata如果打印出None，说明栅格没有声明无效值，这是最常见的隐患之一，海洋、境外区域可能是 0 或者 255 这种占位数字，会被当成真实地类统计进去。np.unique那一段是看值域分布，顺便能发现数据里到底有没有二级编码：如果只出现 1 到 6 这六个数，那这份就是一级类版本，别按二级类去拆。

体检做完，应该能回答三个问题：这是等积投影还是经纬度网格、像元实际物理面积是多少、无效值怎么处理。三个问题都有答案了，再往下算面积才不会翻车。

3. 面积算不对的根源：投影、像元面积与 categorical 重采样

面积统计是这套数据被用得最多的场景，也是错得最多的场景。错法通常有三种：拿经纬度网格直接当等面积网格用、按像元计数乘一个固定系数、重采样的时候用了插值算法。这三种错误各有各的隐蔽性，跑出来的数字看着都很像那么回事。

3.1 Albers 等积投影与经纬度网格：同一块地的两个答案

全国尺度的土地利用栅格，常见做法是采用 Albers 等积圆锥投影，一套典型参数是双标准纬线 25°N 和 47°N、中央经线 105°E。选等积投影的理由很实际：土地利用统计只关心面积，不关心形状和距离是否保真，等积投影能保证任意区域的面积比例不变，做汇总时才不会被投影变形带偏。如果拿到的是发布方给的 Albers 版本，像元是 1000 米 × 1000 米，一个像元就是一平方公里，面积直接乘就行。

麻烦的是经纬度版本。经纬度网格里，一个像元覆盖的东西向距离随纬度收缩，同样“1 度 × 1 度”的方格，在 30°N 和 45°N 覆盖的实际面积能差两成左右。粗算像元物理面积的公式是这样：

import math def cell_area_km2(px_x_deg: float, px_y_deg: float, lat_center: float) -> float: """经纬度网格单个像元的近似面积（km²）""" dx = abs(px_x_deg) * 111.320 * math.cos(math.radians(lat_center)) # 经向实际距离 dy = abs(px_y_deg) * 111.132 # 纬向实际距离 return dx * dy

px_x_deg、px_y_deg是像元的经纬度尺寸，lat_center是像元中心纬度。这个近似在单个像元上误差不大，但做全国汇总时，不同纬度带要用不同系数，逐像元算一遍才准，写起来不难但容易漏。我一般会直接绕开这个问题：先在等积投影下统计，实在要用经纬度版本，就先重投影过去再算面积。

注意：判断依据是src.crs.is_geographic。返回True就是经纬度网格，这时候把像元数乘以一个固定的平方公里数，一定会偏。

3.2 用 rasterio 批量统计分类面积

确定了投影和像元面积，统计就是纯计数。这段代码我用了很多次，改的主要是输入输出路径：

import numpy as np import rasterio def class_area_km2(tif_path: str, pixel_area: float = None) -> dict: """返回 {编码: 面积km²}。pixel_area 为 None 时按米制像元自动推算。""" with rasterio.open(tif_path) as src: arr = src.read(1) nodata = src.nodata # 米制投影下，像元面积 = res_x * res_y，除以 1e6 换成 km² if pixel_area is None: pixel_area = src.res[0] * src.res[1] / 1e6 if nodata is not None: arr = arr[arr != nodata] # 剔除无效值，别让占位数字混进统计 vals, counts = np.unique(arr, return_counts=True) return {int(v): float(c * pixel_area) for v, c in zip(vals, counts)}

逐行说明一下。src.res[0] * src.res[1]得到像元面积，单位是平方米，除以1e6换成平方公里；如果是 1000 米分辨率，结果正好是 1.0。arr != nodata这一步看起来简单，但前提是nodata有值，所以前面体检时发现nodata为None的，要在这之前手动指定一个占位值。最后返回字典，键是编码、值是面积，后面聚合一级类、跨年份比较都基于这个结构。

拿到二级类面积之后聚合到一级类：

def aggregate_level1(area_dict: dict) -> dict: out = {} for code, area in area_dict.items(): k = code // 10 if code >= 10 else code out[k] = out.get(k, 0.0) + area return out

这套写法在图幅边界上会有个细节：跨图幅的像元被重复计数或者漏计。全国数据如果是一整幅，不用管；如果是分幅拼接出来的，统计前要先做一次镶嵌，别对着单幅挨个加。

3.3 重采样到别的分辨率：categorical 只用最近邻

有人为了“让图更好看”或者“和别的 500 米数据对齐”，会把这份 1 公里数据重采样到 500 米。技术上没问题，但必须用最近邻。双线性会在两个编码之间插值，比如 1（耕地）和 5（建设用地）中间插出 3，而 3 在编码表里是草地——凭空多出一个不存在的草地像元，面积统计立刻失真。

import rasterio from rasterio.warp import calculate_default_transform, reproject, Resampling src_path, dst_path = "LU_2015.tif", "LU_2015_500m.tif" with rasterio.open(src_path) as src: transform, width, height = calculate_default_transform( src.crs, src.crs, src.width, src.height, *src.bounds, resolution=500) meta = src.meta.copy() meta.update(transform=transform, width=width, height=height) with rasterio.open(dst_path, "w", **meta) as dst: reproject( source=rasterio.band(src, 1), destination=rasterio.band(dst, 1), src_transform=src.transform, src_crs=src.crs, dst_transform=transform, dst_crs=src.crs, resampling=Resampling.nearest, # 类别数据只能最近邻 )

如果要顺便把投影换成 Albers，把上面的dst_crs从src.crs改成目标 CRS 就行，calculate_default_transform会自动算出新的尺寸和变换参数。

这里有个语义上的坑要说清楚：重采样到 500 米不代表精度变成 500 米。原始 1 公里像元里本来就没有的边界信息，插值插不出来，只是把一张粗网格画成了细网格，图好看了，信息量没增加。做面积统计的时候，重采样还会引入边界像元的类别替换误差，全国范围能差出千分之几，跨年份比较时这个误差如果方向一致，趋势会被轻微带偏。我的习惯是统计永远在原始分辨率上做，重采样只用于出图和叠置显示。

4. 避坑排查：五个把结果带偏的现场

这套数据用久了，翻车的场景会重复出现。下面五个是我自己和身边人踩过最多次的，现象看起来各不相同，根子上大多是“没对齐”或者“没声明清楚”。每一条按现象、原因、解决写，遇到类似的可以先照着对一遍。

4.1 两期数据地图上几乎没差，统计出来差了一大截

现象：某省 1980 和 2015 两期叠在一起看，地类边界几乎没有位移，算出来的耕地转出却是两位数百分比。

原因：两期的网格原点或者尺寸对不上。有些年份的裁剪窗口差了半个像元，地图上肉眼看不出来，数组下标一对齐就成了错位比较，每个像元都在和自己的邻居比。

解决：做比较之前先断言两期网格一致，不一致就用重投影对齐到同一套 transform 上。

import rasterio def assert_same_grid(a_path: str, b_path: str, tol: float = 1e-6) -> None: with rasterio.open(a_path) as a, rasterio.open(b_path) as b: assert a.crs == b.crs, "两期坐标系不一致" for i in range(6): # transform 六个分量逐个比 assert abs(a.transform[i] - b.transform[i]) < tol, f"transform[{i}] 不一致" assert (a.width, a.height) == (b.width, b.height), "行列数不一致"

tol是米制下的容差，1e-6 米对栅格来说可以认为完全相等；如果两期本来就是不同投影，这个断言会先报 CRS 不一致，那就先做重投影再比。

4.2 某县面积逐年跳一下，别的县却平滑

现象：全省趋势平滑，单独看某个县，某一类面积每年跳 5% 到 10%，像有人在改数据。

原因：裁切时没有统一 NoData。按行政区裁剪之后，边界外的区域可能被填成 0，也可能被填成 255，还有可能保留原值。填 0 的话会被当成一个真实存在的“编码 0 地类”算进去，面积自然跳。

解决：裁切时显式指定输出 NoData，统计时用掩膜读取。

with rasterio.open("county_2015.tif") as src: data = src.read(1, masked=True) # masked=True 会把 nodata 自动掩掉 valid_ratio = 1 - data.mask.mean() # 有效像元占比，低于 0.8 就要查裁切参数 print("有效像元占比:", round(float(valid_ratio), 3))

masked=True依赖栅格自己的 nodata 声明，声明缺失时这个掩膜是空的，所以裁切那一步必须把 NoData 写进输出文件，不能只靠统计端补救。

4.3 草地面积是别人软件算的两倍

现象：和同行对数，同一范围同一年的草地面积，自己算出来差不多是人家的两倍。

原因：一级类和二级类混着算。数据里是二级编码，用arr == 3去筛草地一个像元都筛不到；改用(arr >= 30) & (arr <= 33)才是对的。反过来，如果数据是一级类版本，用arr // 10聚合会把 1 到 6 全部映射成 0，统计直接归零。

解决：先确认数据是哪个层级，然后用统一口径聚合，别在同一份统计里混用两种写法。判断方法在 2.3 的值域统计那一步已经做了，值域里最大到 67 就是二级，最大到 6 就是一级。

4.4 转格式或者建金字塔之后，冒出 13、27 这种不存在的编码

现象：把栅格转成别的格式、或者在 GIS 软件里建完金字塔再看直方图，出现了一堆编码表里没有的值。

原因：概览图和格式转换默认用了平均或者双线性重采样。类别栅格的概览层必须是最近邻，否则每一层金字塔都是一次插值污染。

解决：命令行建金字塔时显式指定算法，转格式时也一样。

# 显式用最近邻建金字塔，层级 2/4/8/16 gdaladdo -r nearest LU_2015.tif 2 4 8 16 # 转 GeoPackage 时同样指定，别用默认 gdal_translate -of GPKG -r nearest LU_2015.tif LU_2015.gpkg

已经建错的，先删掉.ovr附属文件再重建，不删的话旧的污染概览层还在。

4.5 NoData 没声明，全国汇总面积明显偏大

现象：把全国所有地类面积加起来，比国土面积大出一截，或者某一年莫名其妙多出几万平方公里。

原因：栅格没有 nodata 属性，海洋、境外、无数据区被填成了 0 或者某个占位数字，统计时当成真实像元计入了。

解决：先看值域，找出那个不该出现的编码，再在脚本里把它当无效值剔除。如果占位值是 0，而编码表里本来就没有 0，那arr = arr[arr > 0]就够了。如果有多个占位值，就显式列出白名单，只保留编码表里存在的值：

VALID_CODES = {11, 12, 21, 22, 23, 24, 31, 32, 33, 41, 42, 43, 44, 45, 46, 51, 52, 53, 61, 62, 63, 64, 65, 66, 67} arr = arr[np.isin(arr, list(VALID_CODES))]

白名单写法比黑名单稳，因为占位值可能不止一个。代价是每次数据口径变化都要改白名单，所以我把它放在配置区，不散落在业务逻辑里。

5. 从两期栅格到转移矩阵：变化检测的完整链路

面积统计回答的是“有多少”，转移矩阵回答的是“从哪来到哪去”。后者才是土地利用变化分析真正有价值的部分：耕地减少了多少不重要，重要的是减少的这部分去了建设用地还是转成了林地，两者的政策含义完全不同。这一章的链路是：先对齐网格，再用交集掩膜，最后做交叉表。

5.1 掩膜先对齐，再谈变化

两期数据直接做逐像元比较，前提是它们在同一套网格上。第 4 章那个断言只能告诉你一不一致，不能帮你对齐。实际做法是把后一期重采样到前一期的网格上，最近邻，然后取两期都有效的交集。

import numpy as np import rasterio from rasterio.warp import reproject, Resampling def load_on_grid(path: str, ref_path: str): """把 path 对齐到 ref_path 的网格，返回 (参考期数组, 对齐后数组)""" with rasterio.open(ref_path) as ref: h, w = ref.height, ref.width transform, crs = ref.transform, ref.crs ref_arr = ref.read(1) with rasterio.open(path) as src: dst = np.zeros((h, w), dtype=src.dtypes[0]) reproject( source=rasterio.band(src, 1), destination=dst, src_transform=src.transform, src_crs=src.crs, dst_transform=transform, dst_crs=crs, resampling=Resampling.nearest, src_nodata=src.nodata, dst_nodata=0, # 无效区统一填 0，方便后面掩掉 ) return ref_arr, dst

函数的第二个返回值是已经对齐到参考网格的数组，dst_nodata=0把重采样产生的空白统一成 0，后面用> 0就能一次性掩掉两期的无效区。如果两期本来就是同一套网格，这一步等于做了一次恒等重采样，开销可以忽略，但保留了代码路径的一致性，我不喜欢为“应该一致”的情况写两条分支。

5.2 用 crosstab 生成转移矩阵

有了对齐的两期数组，交叉表一行代码就能出。但行列标签必须统一，不然某些编码只在第一期出现、不在第二期出现时，对角线会错位，算出来的“未变化”面积是假的。

import numpy as np import pandas as pd def transfer_matrix(a: np.ndarray, b: np.ndarray, valid: np.ndarray) -> pd.DataFrame: """a、b 已对齐的两期数组，valid 为两期都有效的布尔掩膜""" va, vb = a[valid].ravel(), b[valid].ravel() codes = np.union1d(np.unique(va), np.unique(vb)) # 两期编码的并集 tab = pd.crosstab(pd.Series(va, name="t1"), pd.Series(vb, name="t2")) # 用并集重排行列，缺失的组合补 0，保证行列标签顺序一致 return tab.reindex(index=codes, columns=codes, fill_value=0)

np.union1d拿到的是两期所有编码的并集，reindex把交叉表强行对齐到这个并集上，缺的组合填 0。这一步看起来多余，但没有它，某些年份某类地类完全消失时，行列长度不同，np.trace取到的对角线就不是真正的“未变化”面积。调用的时候掩膜这样构造：

y1980, y2015 = load_on_grid("LU_2015.tif", "LU_1980.tif") valid = (y1980 > 0) & (y2015 > 0) # 两期都有效的像元才参与比较 tab = transfer_matrix(y1980, y2015, valid)

5.3 把矩阵读成人话：净变化和总动量要分开看

矩阵出来之后，最容易犯的错是只看净变化。净变化等于行和减列和，正数表示这类地增加了，看起来直观，但它会把方向相反的变化互相抵消。

total = int(tab.values.sum()) unchanged = int(np.trace(tab.values)) print(f"有效像元 {total:,}，未变化 {unchanged:,}，变化占比 {1 - unchanged/total:.2%}") net = tab.sum(axis=1) - tab.sum(axis=0) # 每类净变化，正数为增加 print(net.sort_values())

假设某地一年内 100 平方公里耕地转成了林地，同时 90 平方公里林地转成了耕地，净变化只剩“耕地减少 10 平方公里”。报上去的数字看着很温和，实际上发生了 190 平方公里的地类互换。所以做区域评估时，我一般同时报三个数：净变化、总转出、总转入。总转出等于行和减对角线，总转入等于列和减对角线，这三个数放一起，读者才知道这块地到底“动”了多少。

另外，转移矩阵最好在一级类和二级类两个层级各做一次。一级类矩阵适合对外汇报，结论干净；二级类矩阵能看出耕地内部水田旱地的互换，这种变化在一级类里是看不见的。

6. 收尾：一个脚本把全部年份体检完，顺便做三个守恒检查

一份数据七个年份、三十来个分幅，手工一个个打开看元数据，看完第一个小时人就麻了，而且一定会漏。我现在拿到新一期的土地利用数据，会先跑一个批量体检，把每年每幅的关键属性摊成一张表，再用几条断言卡住明显的异常。表里的字段就是前面几章反复提到的那些：年份、坐标系、行列数、像元面积、有效像元数、总面积。

import pandas as pd def audit(year_map: dict) -> pd.DataFrame: rows = [] for year, path in sorted(year_map.items()): with rasterio.open(path) as src: arr = src.read(1) nodata = src.nodata if src.nodata is not None else 0 valid = arr != nodata px_km2 = src.res[0] * src.res[1] / 1e6 rows.append({ "year": year, "crs": src.crs.to_string() if src.crs else "None", "w": src.width, "h": src.height, "pixel_km2": round(px_km2, 4), "valid_px": int(valid.sum()), "total_km2": round(valid.sum() * px_km2, 1), }) return pd.DataFrame(rows) df = audit(year_map) print(df.to_string(index=False))

三个检查我每次都跑，成本几乎为零，但救过我好几次。第一，所有年份的crs必须只有一个取值，出现两个就说明混进了不同投影的版本，后面统计全废。第二，total_km2的年际波动要在一个合理区间内，全国范围用最大最小比小于 1.05 做粗筛，超过就去看是不是某一年边界外被算进来了。第三，pixel_km2在所有年份必须一致，出现 0.98 或者 1.02 这种，说明有年份被重采样过了，原始分辨率信息已经不可逆地丢了。

assert df["crs"].nunique() == 1, "存在多套坐标系" assert df["total_km2"].max() / df["total_km2"].min() < 1.05, "总面积年际波动异常" assert df["pixel_km2"].nunique() == 1, "像元尺寸不一致，可能被重采样过"

断言失败的时候不要急着改阈值，先去看失败那一年的原始文件，绝大多数情况是分发批次不同、参数被打乱，重新取一份对应版本的就行。阈值本身是经验值，全国尺度上海洋和境外区域基本不变，1.05 的宽松度足够覆盖真实的地类变化。

血泪经验是，最开始我图省事，把面积统计和变化检测写在一个脚本里，结果中间某一期投影不对，一路算到最后才报错，排查花了整整一个下午。从那以后我每次拿到新一期的土地利用数据，都强制先跑一遍这个体检表和三条断言，通过了再进分析流程。看着啰嗦，其实是给自己留了后悔药。希望帮到你。

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