简介：这份资源为2018年全国土地利用30米分辨率遥感数据，面向GIS、遥感、城乡规划、生态环保等方向的研究人员与学生，用于土地覆盖分类、时空变化分析与制图实践。数据以30米栅格像元刻画耕地、林地、草地、建设用地、水域等地类，遵循国家土地分类标准，可直接用于空间统计与趋势研判。压缩包共10个文件，约811.49MB，包含tif栅格主数据、dbf属性表、tfw坐标信息、ovr金字塔、xml元数据、pdf说明文档、xlsx分类标准及jpg色标参考，兼顾数据读取、分类对照与可视化需求。目前已有9295人学习下载，说明其在教学与科研中具有较高参考价值。读者可借此掌握遥感土地利用数据的组织方式、分类体系与GIS处理流程，为城市扩张、耕地保护、生态评估等课题提供基础数据支撑。
1. 遥感全国土地利用30m数据：从一张图到一套可复现的落地链路
手里有一份全国范围、30m 分辨率的土地利用栅格，第一反应往往不是“真好看”，而是“这玩意儿怎么用”。全国 30m 土地利用数据（业内常说的 CLCD 系列就是典型代表）本质是一张覆盖国土的类别栅格，每个像元存一个类别码，常见的是耕地、林地、草地、水域、建设用地、未利用地六大类，也有按一级/二级分类体系细分的版本。它解决的核心问题是：在不需要自己跑分类模型的前提下，拿到一份时间序列一致、空间无缝、类别体系统一的底图，用来做变化检测、生态遥感指数计算、统计报表或者给遥感图像标注做先验。适合谁？做国土空间规划、生态评估、农业遥感、碳汇估算的从业者，以及想拿现成标签训练遥感随机森林、SegFormer 这类模型的学生和工程师。这一章先把“它是什么、边界在哪”讲清楚，后面几章再动手。
2. 拿到数据先别急着算：30m 土地利用的坐标、分类与质量核验
2.1 为什么 30m 分辨率决定了你能做什么、不能做什么
30m 这个数字不是随便定的，它对应 Landsat 系列多光谱影像的空间分辨率，也是国内土地利用产品最主流的一档。一个 30m×30m 的像元，实际覆盖 900 平方米，约等于 1.35 亩。这个尺度意味着：你能可靠识别成片农田、连片林地、大中型水体、城市建成区，但识别不了一条乡道、一栋独立小楼、一条田埂。很多新手拿 30m 数据去做地块级农田识别，结果边界糊成一团，这就是尺度错配。
选型上要分清两类需求。第一类是统计与趋势分析，比如算某县十年间建设用地扩张速率、林地转耕地的面积，30m 完全够用，全国覆盖、逐年更新、类别一致，比你自己拼影像分类省几个月。第二类是精细制图，比如高分遥感影像农田地块智能识别，那需要亚米级或米级影像，30m 只能当先验掩膜，用来缩小搜索范围，不能当最终边界。
还有一个容易被忽略的点：30m 产品的类别精度在不同区域差异很大。平原农业区耕地提取通常很稳，山区林地与草地、灌丛的混分就明显增多，城乡结合部的建设用地和裸地也容易互相串。所以拿到数据后，第一件事不是跑统计，而是做区域性的质量核验。
2.2 坐标系统与投影：别让统计面积悄悄偏了
全国 30m 土地利用栅格常见的坐标系是地理坐标（如 CGCS2000 或 WGS84），单位是度，像元大小约 0.000269 度。问题来了：地理坐标下每个像元的实际地面面积随纬度变化，直接按像元计数乘 900 平方米算面积，在高纬度会偏。正确做法是先投影到等面积投影再统计。
下面这段用 Python 做投影转换和面积统计，GDAL 和 rasterio 都行，这里用 rasterio 演示：
import rasterio from rasterio.warp import calculate_default_transform, reproject, Resampling import numpy as np src_path = "CLCD_2020_national.tif" dst_path = "CLCD_2020_albers.tif" # 目标：Albers 等面积投影，适合全国尺度面积统计 dst_crs = "+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" with rasterio.open(src_path) as src: transform, width, height = calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds) kwargs = src.meta.copy() kwargs.update({ "crs": dst_crs, "transform": transform, "width": width, "height": height, "compress": "lzw" # 全国数据量大，压缩能省一半以上空间 }) with rasterio.open(dst_path, "w", **kwargs) 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=dst_crs, resampling=Resampling.nearest # 类别栅格必须用最近邻，不能双线性 )逻辑说明：calculate_default_transform自动算出目标投影下的范围和像元数；Resampling.nearest是关键，类别码是离散整数，用双线性或立方插值会造出 3.5 这种不存在的类别。参数上，lat_1、lat_2是 Albers 标准纬线，全国常用 25 和 47，lon_0取 105 居中。投影后像元约 30m，此时按像元计数乘 900 平方米才准。
2.3 分类体系对齐：一级类和二级类别混用
不同来源的土地利用数据分类体系不一样。常见的一级类六类，二级类可能到二十多类。做统计前必须确认你手里这份的类别码定义，别拿一套码表去套另一套数据。下面这张对照表是我常用的核对方式：
| 类别码 | 一级类 | 常见二级类 | 统计注意 |
|---|---|---|---|
| 1 | 耕地 | 水田、旱地 | 山区旱地与草地易混 |
| 2 | 林地 | 有林地、灌木、疏林 | 与灌丛边界模糊 |
| 3 | 草地 | 高覆盖、中低覆盖 | 与未利用地交叉 |
| 4 | 水域 | 河渠、湖泊、水库 | 季节性水体易漏 |
| 5 | 建设用地 | 城镇、农村、工矿 | 城乡结合部偏多 |
| 6 | 未利用地 | 裸地、沙地、盐碱 | 与建设用地互串 |
核验方法很直接：裁一块你熟悉的区域，把栅格类别和卫星底图叠着看，重点看边界地带。发现系统性偏移，要么是坐标系没对齐，要么是类别码理解错了。
3. 用 Python 把全国 30m 土地利用跑成统计表：分块、掩膜与三大类汇总
3.1 全国数据不能一次性读进内存：分块读取的正确姿势
全国 30m 栅格动辄几十 GB，直接read()整幅进内存，机器内存不够就崩。血泪经验是必须分块（window）处理。rasterio 的block_windows或者手动切 window 都行。下面按行政边界做分区统计的骨架：
import rasterio import numpy as np import geopandas as gpd from rasterio.mask import mask landuse_path = "CLCD_2020_albers.tif" boundary_path = "county_boundary.shp" # 县级行政边界，需与栅格同投影 gdf = gpd.read_file(boundary_path) with rasterio.open(landuse_path) as src: for idx, row in gdf.iterrows(): geom = [row.geometry.__geo_interface__] try: out_image, out_transform = mask(src, geom, crop=True, nodata=0) except ValueError: continue # 边界与栅格无交集 data = out_image[0] # 统计每个类别像元数 classes, counts = np.unique(data[data > 0], return_counts=True) area = {int(c): int(n) * 900 / 1e6 for c, n in zip(classes, counts)} # 平方公里 print(row.get("NAME", idx), area)逻辑说明：mask按几何裁剪，crop=True只返回边界外接矩形范围，省内存；nodata=0把边界外置零，统计时用data > 0过滤。参数上，900是投影后像元面积（平方米），除以1e6转平方公里。注意边界 shp 必须和栅格同投影，否则 mask 会报错或裁出空图。
3.2 土地利用三大类分类统计工具：耕地、生态、建设用地的归并逻辑
热搜里常出现“土地利用三大类分类统计工具”，本质是把六类归并成三大功能类：生产类（耕地）、生态类（林地+草地+水域）、生活类（建设用地），未利用地单列或并入生态。归并不是简单相加，要按你的分析目标定。做生态遥感指数，水域和林草地要分开权重；做建设用地扩张，未利用地转建设用地的部分要单独追踪。
# 六类 -> 三大类映射 mapping = { 1: "production", # 耕地 2: "ecological", # 林地 3: "ecological", # 草地 4: "ecological", # 水域 5: "living", # 建设用地 6: "other" # 未利用地 } def summarize_three(area_dict): result = {"production": 0, "ecological": 0, "living": 0, "other": 0} for code, km2 in area_dict.items(): result[mapping.get(code, "other")] += km2 return result逻辑说明：映射表是核心，改一个类别归属，整张统计表就变。参数上，mapping用字典便于替换；如果数据是二级类，先把二级码归到一级码再套这张表。这一步做完，你就能输出每个行政单元的三大类面积和占比，直接进报表。
3.3 变化检测：两期栅格相减得到转移矩阵
土地利用最有价值的用法是变化检测。两期同投影、同类别码的栅格，逐像元比较就能得到转移矩阵。注意必须保证两期数据分类体系和投影完全一致，否则矩阵全是噪声。
import numpy as np import rasterio with rasterio.open("CLCD_2010_albers.tif") as s1, \ rasterio.open("CLCD_2020_albers.tif") as s2: a1 = s1.read(1) a2 = s2.read(1) valid = (a1 > 0) & (a2 > 0) # 组合成 from*10+to 的编码，统计转移 combo = a1[valid].astype(np.int32) * 10 + a2[valid].astype(np.int32) codes, counts = np.unique(combo, return_counts=True) for c, n in zip(codes, counts): frm, to = divmod(int(c), 10) if frm != to: print(f"{frm} -> {to}: {n * 900 / 1e6:.2f} km2")逻辑说明：a1*10+a2把“从哪类到哪类”编码成一个整数，divmod拆回来。参数上，乘 10 是因为类别码是个位数，若类别码超过 9 要改成乘 100。valid掩膜排除两期任一为 nodata 的像元，否则边界会造出假变化。这一步跑完，耕地转建设、林地转耕地这些关键流向一目了然。
4. 避坑与排查：30m 土地利用落地时最容易翻车的五件事
4.1 现象：统计面积和官方公布对不上，差百分之十几
原因：多半是投影没转，直接在地理坐标下按像元计数。纬度越高，像元实际面积越小，全国平均下来偏差可观。另一个原因是边界裁剪时用了外接矩形而非精确掩膜，把边界外像元算进来了。
解决：先转 Albers 等面积投影再统计；裁剪用rasterio.mask精确到几何，别用 bounding box。核验时拿一个你熟悉的县，和统计年鉴对一下，偏差应在几个百分点内。
4.2 现象：类别栅格重采样后出现一堆没见过的类别码
原因：用了双线性或立方插值。类别码是离散整数，插值会算出 2.7、4.3 这种值，取整后变成乱七八糟的码。
解决：类别栅格的一切重采样、投影转换，resampling必须设nearest。这条没有例外，遥感数字图像处理里类别图和连续量图的处理逻辑是两套。
4.3 现象：变化检测结果里出现大量“未利用地转建设用地”的假变化
原因：两期数据分类体系或版本不一致，或者其中一期在城乡结合部把裸地标成了建设用地。也可能是 nodata 处理不一致，边界像元被当成真实变化。
解决：先确认两期数据同源同版本；用valid掩膜排除 nodata；对城乡结合部做抽样目视核验，必要时用高分影像辅助判断。假变化集中的区域，往往是分类精度最弱的地方。
4.4 现象：全国数据跑统计跑到一半内存爆掉
原因：一次性read()整幅，或者用 GeoPandas 把全国边界和栅格做叠加时没分块。
解决：按 window 或按行政单元逐个裁剪处理，处理完立即释放；用crop=True减少返回数据量；全国任务建议拆成省或县并行跑，单机也能扛。
4.5 现象：拿 30m 数据训练分割模型，精度上不去
原因：把 30m 类别栅格当成了像素级标签去训练 SegFormer 这类模型，但 30m 标签本身边界模糊，和米级影像的像素对不齐，模型学到的是噪声边界。
解决：30m 数据更适合当粗标签或先验掩膜，训练时降采样或做标签松弛；精细分割要用高分影像加人工标注，30m 只用来约束大区域类别。遥感图像标注里，标签尺度和影像尺度匹配是第一条原则。
5. 进阶：把 30m 土地利用接进生态遥感指数与模型训练链路
走到这一步，数据能统计、能检测变化了，接下来是把它变成更大分析链路的一环。最常见的两个方向：生态遥感指数计算，和给遥感随机森林、SegFormer 提供先验。
生态遥感指数（如植被覆盖度、遥感生态指数 RSEI）通常需要土地利用做掩膜或分区。比如算植被覆盖度时，用土地利用把建设用地和水域剔掉，只对林草耕统计，结果更干净。做法是把 30m 类别栅格重采样到指数栅格的分辨率（同样用 nearest），生成布尔掩膜：
import rasterio import numpy as np from rasterio.warp import reproject, Resampling # 把土地利用重采样到 NDVI 栅格网格，生成植被区掩膜 with rasterio.open("landuse.tif") as lu, rasterio.open("ndvi.tif") as nd: lu_resampled = np.empty(nd.shape, dtype=np.uint8) reproject( source=rasterio.band(lu, 1), destination=lu_resampled, src_transform=lu.transform, src_crs=lu.crs, dst_transform=nd.transform, dst_crs=nd.crs, resampling=Resampling.nearest ) veg_mask = np.isin(lu_resampled, [1, 2, 3]) # 耕林草 ndvi = nd.read(1) ndvi_veg = np.where(veg_mask, ndvi, np.nan) print("植被区 NDVI 均值:", np.nanmean(ndvi_veg))逻辑说明：reproject把类别栅格对齐到 NDVI 网格，np.isin生成植被布尔掩膜，np.where把非植被区置 NaN，nanmean忽略 NaN 求均值。参数上，类别列表[1,2,3]按你的分类体系改；如果 NDVI 有缩放因子（如乘了 10000），记得先还原。
给模型训练用时，思路反过来：把 30m 类别作为弱标签，训练一个粗分割模型，再用它去预标注高分影像，人工只做修正。这就是局部聚焦算法辅助标记的高分遥感影像农田地块智能识别的常见套路——30m 给大区域先验，人工聚焦在边界疑难处，标注效率能提不少。但记住，30m 标签不能直接当高分影像的像素级真值，中间必须有人工核验环节，否则误差会被模型放大。
我自己踩过最深的一个坑，是早期图省事，直接在地理坐标下统计全国耕地面积，结果和年鉴差了近一成，排查了两天才发现是投影问题。从那以后，任何栅格统计前先问一句：投影转了吗，重采样用的 nearest 吗。这两个习惯帮我省了无数后悔药。希望帮到你。
本文还有配套的精品资源，点击获取
