简介这份资源为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 projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm 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( sourcerasterio.band(src, 1), destinationrasterio.band(dst, 1), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingResampling.nearest # 类别栅格必须用最近邻不能双线性 )逻辑说明calculate_default_transform自动算出目标投影下的范围和像元数Resampling.nearest是关键类别码是离散整数用双线性或立方插值会造出 3.5 这种不存在的类别。参数上lat_1、lat_2是 Albers 标准纬线全国常用 25 和 47lon_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, cropTrue, nodata0) except ValueError: continue # 边界与栅格无交集 data out_image[0] # 统计每个类别像元数 classes, counts np.unique(data[data 0], return_countsTrue) area {int(c): int(n) * 900 / 1e6 for c, n in zip(classes, counts)} # 平方公里 print(row.get(NAME, idx), area)逻辑说明mask按几何裁剪cropTrue只返回边界外接矩形范围省内存nodata0把边界外置零统计时用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*10to 的编码统计转移 combo a1[valid].astype(np.int32) * 10 a2[valid].astype(np.int32) codes, counts np.unique(combo, return_countsTrue) 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*10a2把“从哪类到哪类”编码成一个整数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 或按行政单元逐个裁剪处理处理完立即释放用cropTrue减少返回数据量全国任务建议拆成省或县并行跑单机也能扛。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, dtypenp.uint8) reproject( sourcerasterio.band(lu, 1), destinationlu_resampled, src_transformlu.transform, src_crslu.crs, dst_transformnd.transform, dst_crsnd.crs, resamplingResampling.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把非植被区置 NaNnanmean忽略 NaN 求均值。参数上类别列表[1,2,3]按你的分类体系改如果 NDVI 有缩放因子如乘了 10000记得先还原。给模型训练用时思路反过来把 30m 类别作为弱标签训练一个粗分割模型再用它去预标注高分影像人工只做修正。这就是局部聚焦算法辅助标记的高分遥感影像农田地块智能识别的常见套路——30m 给大区域先验人工聚焦在边界疑难处标注效率能提不少。但记住30m 标签不能直接当高分影像的像素级真值中间必须有人工核验环节否则误差会被模型放大。我自己踩过最深的一个坑是早期图省事直接在地理坐标下统计全国耕地面积结果和年鉴差了近一成排查了两天才发现是投影问题。从那以后任何栅格统计前先问一句投影转了吗重采样用的 nearest 吗。这两个习惯帮我省了无数后悔药。希望帮到你。本文还有配套的精品资源点击获取