简介本资源为中国土地利用现状遥感监测数据合集面向地理信息、遥感解译、国土空间规划及生态环境研究领域的科研人员与高年级学生可支撑长时序土地利用变化分析、区域对比与制图建模等场景。数据覆盖1980、1990、1995、2000、2005、2010、2015及2020年来源为资源环境科学与数据中心其中2020年数据基于Landsat 8影像在2015年基础上人工目视解译生成全国范围已完成。土地利用类型包含耕地、林地、草地、水域、居民地和未利用土地6个一级类型及25个二级类型分类体系完整便于直接开展统计与空间分析。压缩包共163个文件以adf、dat、nit等ArcInfo栅格与属性文件为主辅以rar分卷、log日志、doc说明及xml元数据整体约26.25MB目录结构清晰便于按年份与区域检索。目前已有4497人学习下载适合需要长时序土地利用基础数据的研究者参考使用。1. 拿到「中国土地利用现状遥感监测数据.rar」先别急着解压它到底能干什么如果你手里正躺着一个叫「中国土地利用现状遥感监测数据.rar」的压缩包第一反应大概率是双击解压然后对着一堆命名规整的栅格文件发懵。这份数据在国内土地覆被、生态评估、耕地变化、城市扩张研究里出现频率极高很多论文、报告、国土空间规划的基础底图就是它。它解决的核心问题是在统一分类体系下把全国每一块地表覆盖类型落到像元级别让你能按年份、按区域统计耕地、林地、草地、水域、建设用地和未利用地的面积与空间分布。适合谁做遥感应用、地理信息、生态遥感、土地资源管理、双碳核算的从业者以及需要一份权威底图来验证自己分类算法的人。但先别急着跑统计这份数据的坑从解压那一刻就开始了。2. 先搞懂分类体系和年份口径不然统计出来的耕地面积全是错的2.1 这套数据到底分了几类一级类和二级类怎么对应中国土地利用现状遥感监测数据最常见的版本采用两级分类体系。一级类通常是 6 个耕地、林地、草地、水域、建设用地、未利用地。二级类在一级类下细分比如耕地再分水田、旱地林地分有林地、灌木林、疏林地等。不同年份、不同生产批次二级类的数量和命名会有细微差别这是第一个容易翻车的地方。我一般拿到数据后第一件事不是打开 ArcGIS而是先找到随数据附带的分类说明文件。如果压缩包里没有就去翻数据生产单位的公开文档。分类对照关系必须落到一张表里否则后面做面积统计时你会把「有林地」和「灌木林」混在一起或者把「水库水面」当成「滩涂」。一级类常见二级类编码示例统计时注意耕地水田、旱地11、12水田和旱地灌溉条件不同碳核算要分开林地有林地、灌木林、疏林地21、22、23疏林地容易和草地混淆草地高覆盖、中覆盖、低覆盖31、32、33覆盖度阈值不同年份可能调整水域河渠、湖泊、水库、滩涂41、42、43、44滩涂季节性变化大建设用地城镇、农村、工矿、交通51、52、53、54农村建设用地口径变化最频繁未利用地沙地、戈壁、盐碱地、沼泽61、62、63、64沼泽有时被归到水域这张表不是让你背而是让你在写统计脚本时有个对照依据。编码体系一旦搞错后面所有面积数字都是空中楼阁。2.2 年份口径和分辨率为什么同一块地两年面积对不上这套数据通常按年份发布常见的有 1980 年代末、1990 年代、2000 年、2005 年、2010 年、2015 年、2020 年等节点。每个年份的生产方式、影像源、分类标准可能有调整。比如早期版本依赖 Landsat TM后期加入国产高分影像分类精度和最小图斑面积都会变。分辨率方面常见的是 1km 栅格和 30m 栅格两种。1km 数据适合全国尺度趋势分析30m 数据适合省域或流域尺度。如果你拿 1km 数据去算一个县的耕地变化边界混合像元会严重干扰结果。我一般会先确认三件事投影坐标系是什么、像元大小是多少、NoData 值填的是什么。这三项不确认后面做掩膜和重采样必出问题。提示不同年份的数据不要直接相减做变化检测先统一投影、重采样到同一网格再对齐分类编码。3. 用 Python 把 rar 解压、读取、裁切到研究区一套能复现的最小流程3.1 解压和目录结构梳理别用鼠标一个个点拿到 .rar 文件Windows 上很多人习惯用 WinRAR 手动解压。但如果你的研究区涉及多个年份、多个分幅手动操作迟早出错。我一般用命令行工具批量处理。Linux 和 macOS 下用 unrarWindows 下可以用 7z 命令行。# 查看压缩包内文件列表先不急着解压 unrar l 中国土地利用现状遥感监测数据.rar # 解压到指定目录保留目录结构 unrar x 中国土地利用现状遥感监测数据.rar ./landuse_data/ # 如果压缩包分卷确保所有分卷在同一目录 # unrar x 中国土地利用现状遥感监测数据.part1.rar ./landuse_data/解压后常见目录结构是按年份分文件夹每个年份下再按省份或分幅存放 GeoTIFF 或 IMG 文件。先别急着写代码用find或dir把文件清单导出来确认命名规律。# 导出文件清单方便后续批量处理 find ./landuse_data -name *.tif -o -name *.img file_list.txt wc -l file_list.txt逻辑说明unrar l只列出内容不写磁盘适合先确认压缩包完整性。unrar x保留目录结构避免文件散落一地。find同时匹配 tif 和 img是因为不同批次数据格式可能不统一。参数上如果你的 rar 有密码加-p参数但这类公开数据一般没有。3.2 用 rasterio 读取并裁切到研究区边界读取栅格我首选 rasterio比 GDAL 原生绑定更友好和 numpy 配合也顺。下面这段代码完成三件事读取土地利用栅格、读取研究区矢量边界、按边界裁切并统计各类型面积。import rasterio from rasterio.mask import mask import geopandas as gpd import numpy as np from collections import Counter # 1. 读取研究区边界 study_area gpd.read_file(./boundary/study_area.shp) # 确保坐标系一致这里假设边界是 WGS84 study_area study_area.to_crs(EPSG:4326) # 2. 打开土地利用栅格 with rasterio.open(./landuse_data/2020/landuse_2020.tif) as src: # 打印元数据确认投影、像元大小、NoData print(src.crs, src.res, src.nodata) # 3. 按边界裁切 geometries study_area.geometry.values out_image, out_transform mask(src, geometries, cropTrue) out_meta src.meta.copy() out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) # 4. 统计各类型像元数 data out_image[0] # 排除 NoData valid data[data ! src.nodata] counts Counter(valid.flatten()) for code, cnt in sorted(counts.items()): print(f编码 {int(code)}: {cnt} 个像元) # 5. 像元数转面积假设像元 30m pixel_area 30 * 30 # 平方米 for code, cnt in sorted(counts.items()): area_km2 cnt * pixel_area / 1e6 print(f编码 {int(code)}: {area_km2:.2f} 平方公里)逻辑说明mask函数按矢量边界裁切栅格cropTrue会缩小输出范围减少内存占用。src.nodata必须读取否则 NoData 会被当成有效值参与统计。像元面积换算时如果数据是地理坐标系经纬度30m 不是固定值需要投影到等面积坐标系再算。参数上study_area.to_crs要和栅格 CRS 一致不一致会裁切出空数组。注意如果研究区跨多个分幅先做镶嵌再裁切不要逐幅裁切后合并否则边界像元统计会重复。3.3 批量处理多年份数据并导出面积表单年份跑通后把流程封装成函数批量处理所有年份。下面这段代码遍历年份文件夹输出一张 CSV 面积表。import os import pandas as pd def stats_by_year(tif_path, boundary, nodataNone): with rasterio.open(tif_path) as src: if nodata is None: nodata src.nodata out_image, _ mask(src, boundary.geometry.values, cropTrue) data out_image[0] valid data[data ! nodata] counts Counter(valid.flatten()) pixel_area src.res[0] * src.res[1] records [] for code, cnt in counts.items(): records.append({ code: int(code), pixel_count: cnt, area_km2: cnt * pixel_area / 1e6 }) return pd.DataFrame(records) years [1990, 2000, 2010, 2020] all_stats [] for y in years: tif f./landuse_data/{y}/landuse_{y}.tif if os.path.exists(tif): df stats_by_year(tif, study_area) df[year] y all_stats.append(df) result pd.concat(all_stats, ignore_indexTrue) result.to_csv(./output/landuse_area_by_year.csv, indexFalse) print(result.head())逻辑说明把单年份逻辑封装成函数方便复用。src.res返回像元宽高相乘得像元面积。如果数据是地理坐标系res单位是度不能直接乘必须先投影。参数上nodata允许外部传入因为有些文件元数据里 NoData 标记不规范需要手动指定比如 0 或 255。4. 投影、重采样和边界对齐三个最容易让面积统计翻车的操作4.1 地理坐标系直接算面积为什么你的平方公里数偏小很多人拿到栅格看元数据里写着 WGS84就直接用像元大小乘个数算面积。这是血泪经验里最常见的一条。经纬度坐标系下一个像元在地表的东西方向距离随纬度变化赤道附近约 30m到了北纬 40 度只剩约 23m。你直接按 30m 算高纬度地区面积会偏大低纬度地区偏小全国尺度上误差能到百分之十几。正确做法是先投影到等面积坐标系。国内常用 Albers 等面积投影中央经线一般取 105°E双标准纬线取 25°N 和 47°N。用 rasterio 或 gdalwarp 都可以。# 用 gdalwarp 投影到 Albers 等面积投影 gdalwarp -t_srs projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm no_defs \ -tr 30 30 -r near \ -srcnodata 255 -dstnodata 255 \ ./landuse_data/2020/landuse_2020.tif \ ./landuse_data/2020/landuse_2020_albers.tif逻辑说明-t_srs指定目标投影-tr 30 30设置输出像元 30m-r near表示最近邻重采样分类数据必须用最近邻不能用双线性否则会出现不存在的类别编码。-srcnodata和-dstnodata保持一致避免 NoData 被重采样成有效值。4.2 重采样时分类数据只能用最近邻双线性会造出假类别分类栅格和连续栅格的重采样逻辑完全不同。NDVI 这种连续值可以用双线性或立方卷积但土地利用编码是离散整数双线性会把 11 和 12 平均成 11.5这个编码在分类体系里根本不存在。我见过有人重采样后统计出编码 15、25 这种奇怪值就是插值惹的祸。最近邻重采样的缺点是边界像元会偏移但至少类别是真实的。如果你需要更精细的边界应该用更高分辨率数据而不是靠插值。重采样方法适用数据类型对分类数据的影响最近邻分类、整数编码保持原始编码边界略偏移双线性连续值、NDVI产生不存在的编码值立方卷积连续值、影像同上且更平滑众数分类数据保持编码但计算量大4.3 多期数据网格对齐像元不重合导致变化检测全错做变化检测时两期数据的网格必须完全对齐。如果一期是 30m另一期是 1km或者投影参数差一点像元边界就错开了。你直接相减会得到大量伪变化。正确做法是以某一期为基准网格把其他期重采样到同一网格。from rasterio.warp import reproject, Resampling # 以 2020 年为基准网格 with rasterio.open(./landuse_data/2020/landuse_2020_albers.tif) as ref: ref_meta ref.meta.copy() ref_data ref.read(1) # 把 2010 年数据重采样到基准网格 with rasterio.open(./landuse_data/2010/landuse_2010_albers.tif) as src: dst np.empty_like(ref_data) reproject( sourcerasterio.band(src, 1), destinationdst, src_transformsrc.transform, src_crssrc.crs, dst_transformref_meta[transform], dst_crsref_meta[crs], resamplingResampling.nearest ) # 保存对齐后的数据 ref_meta.update({dtype: uint8, nodata: 255}) with rasterio.open(./landuse_data/2010/landuse_2010_aligned.tif, w, **ref_meta) as dst_file: dst_file.write(dst, 1)逻辑说明reproject把源数据映射到目标网格Resampling.nearest保证分类编码不变。dst_transform和dst_crs来自基准数据这样两期数据每个像元位置一一对应。参数上dst数组形状必须和基准数据一致否则写入会报错。5. 避坑与排查土地利用遥感监测数据最常见的 5 个翻车现场5.1 现象统计出的建设用地面积逐年下降原因不同年份的建设用地编码口径变了。早期版本农村建设用地可能归在耕地附属后期单独编码。或者你用的分类对照表只对上了一级类二级类编码错位。解决逐年核对分类说明建立年份-编码映射表。不要用一套编码表跑所有年份。如果找不到说明用频率分布对比看哪些编码在某年突然出现或消失。5.2 现象裁切后研究区边界外出现大量有效值原因矢量边界和栅格坐标系不一致或者边界有自相交、空洞。mask函数默认按几何范围裁切如果几何不闭合会裁出多余区域。解决先study_area study_area.buffer(0)修复几何再to_crs对齐坐标系。裁切后用边界再做一次精确掩膜把外部像元设为 NoData。5.3 现象面积统计结果比官方公布数据大很多原因像元面积算错最常见的是地理坐标系直接乘或者像元大小读的是经纬度度数而不是米。另一个原因是 NoData 没排除被当成未利用地统计了。解决确认投影单位投影到等面积坐标系后再算。检查src.nodata如果元数据里是 None手动指定常见值如 0、255、-9999逐个测试。5.4 现象变化检测结果里出现大量 11 变 11 的伪变化原因两期数据网格没对齐或者分类编码版本不同但编码值相同。比如 2010 年的 11 是水田2020 年的 11 可能变成了旱地。解决先做网格对齐再核对两期分类体系。如果编码含义变了先重编码到统一体系再做变化检测。变化检测前用np.where排除 NoData 和相同值。5.5 现象大文件读取时内存溢出原因全国 30m 数据单年份可能几十 GB直接src.read()全读进内存会爆。解决用窗口读取分块处理。rasterio 支持block_windows按块统计后累加。或者先用gdalwarp裁切到研究区再读取小文件。# 分块统计示例 with rasterio.open(tif_path) as src: counts Counter() for ji, window in src.block_windows(1): block src.read(1, windowwindow) valid block[block ! src.nodata] counts.update(valid.flatten())逻辑说明block_windows按栅格内部块遍历每次只读一块内存占用可控。counts.update累加各块统计结果。参数上块大小由数据创建时决定一般 256x256 或 512x512。6. 从面积统计到空间格局一个我常用的验证技巧面积统计只是第一步真正让这份数据发挥价值的是空间格局分析。我一般会做一个简单但有效的交叉验证把统计出的各类型面积和同区域、同年份的公开统计年鉴数据对比。不是要求完全一致而是看量级和趋势是否合理。如果耕地面积差了一个数量级那肯定是投影或编码出了问题。进阶用法上我习惯把土地利用数据和 DEM、NDVI、夜间灯光数据叠加做地形梯度分析或城市扩张分析。比如用 DEM 提取坡度统计不同坡度带上的耕地分布能快速识别陡坡开垦区域。下面这段代码演示如何按坡度带统计耕地像元。import rasterio import numpy as np from rasterio.mask import mask # 读取土地利用和 DEM with rasterio.open(./landuse_albers.tif) as lu_src: lu_data, lu_transform mask(lu_src, study_area.geometry.values, cropTrue) lu lu_data[0] with rasterio.open(./dem_albers.tif) as dem_src: dem_data, _ mask(dem_src, study_area.geometry.values, cropTrue) dem dem_data[0] # 计算坡度这里简化用梯度近似 dy, dx np.gradient(dem) slope np.degrees(np.arctan(np.sqrt(dx**2 dy**2))) # 定义坡度带 slope_bins [0, 2, 6, 15, 25, 90] slope_labels [0-2, 2-6, 6-15, 15-25, 25] # 统计耕地编码 11、12在各坡度带的分布 farmland_mask np.isin(lu, [11, 12]) for i in range(len(slope_bins)-1): mask_bin (slope slope_bins[i]) (slope slope_bins[i1]) count np.sum(farmland_mask mask_bin) print(f坡度 {slope_labels[i]} 度: {count} 个耕地像元)逻辑说明np.gradient计算 DEM 在 x 和 y 方向的梯度arctan转成角度得到坡度。np.isin筛选耕地编码。参数上坡度带划分依据不同研究目的调整做水土流失评估时常用 15 度和 25 度作为阈值。这个技巧的好处是你不需要复杂的模型就能从土地利用数据里挖出有生态意义的信息。我做了这么多年最深的教训就是别一上来就追求花哨的算法先把投影、编码、NoData 这三件事做对比什么都强。希望帮到你。本文还有配套的精品资源点击获取
中国土地利用遥感监测数据实战：从解压到面积统计的避坑指南
简介本资源为中国土地利用现状遥感监测数据合集面向地理信息、遥感解译、国土空间规划及生态环境研究领域的科研人员与高年级学生可支撑长时序土地利用变化分析、区域对比与制图建模等场景。数据覆盖1980、1990、1995、2000、2005、2010、2015及2020年来源为资源环境科学与数据中心其中2020年数据基于Landsat 8影像在2015年基础上人工目视解译生成全国范围已完成。土地利用类型包含耕地、林地、草地、水域、居民地和未利用土地6个一级类型及25个二级类型分类体系完整便于直接开展统计与空间分析。压缩包共163个文件以adf、dat、nit等ArcInfo栅格与属性文件为主辅以rar分卷、log日志、doc说明及xml元数据整体约26.25MB目录结构清晰便于按年份与区域检索。目前已有4497人学习下载适合需要长时序土地利用基础数据的研究者参考使用。1. 拿到「中国土地利用现状遥感监测数据.rar」先别急着解压它到底能干什么如果你手里正躺着一个叫「中国土地利用现状遥感监测数据.rar」的压缩包第一反应大概率是双击解压然后对着一堆命名规整的栅格文件发懵。这份数据在国内土地覆被、生态评估、耕地变化、城市扩张研究里出现频率极高很多论文、报告、国土空间规划的基础底图就是它。它解决的核心问题是在统一分类体系下把全国每一块地表覆盖类型落到像元级别让你能按年份、按区域统计耕地、林地、草地、水域、建设用地和未利用地的面积与空间分布。适合谁做遥感应用、地理信息、生态遥感、土地资源管理、双碳核算的从业者以及需要一份权威底图来验证自己分类算法的人。但先别急着跑统计这份数据的坑从解压那一刻就开始了。2. 先搞懂分类体系和年份口径不然统计出来的耕地面积全是错的2.1 这套数据到底分了几类一级类和二级类怎么对应中国土地利用现状遥感监测数据最常见的版本采用两级分类体系。一级类通常是 6 个耕地、林地、草地、水域、建设用地、未利用地。二级类在一级类下细分比如耕地再分水田、旱地林地分有林地、灌木林、疏林地等。不同年份、不同生产批次二级类的数量和命名会有细微差别这是第一个容易翻车的地方。我一般拿到数据后第一件事不是打开 ArcGIS而是先找到随数据附带的分类说明文件。如果压缩包里没有就去翻数据生产单位的公开文档。分类对照关系必须落到一张表里否则后面做面积统计时你会把「有林地」和「灌木林」混在一起或者把「水库水面」当成「滩涂」。一级类常见二级类编码示例统计时注意耕地水田、旱地11、12水田和旱地灌溉条件不同碳核算要分开林地有林地、灌木林、疏林地21、22、23疏林地容易和草地混淆草地高覆盖、中覆盖、低覆盖31、32、33覆盖度阈值不同年份可能调整水域河渠、湖泊、水库、滩涂41、42、43、44滩涂季节性变化大建设用地城镇、农村、工矿、交通51、52、53、54农村建设用地口径变化最频繁未利用地沙地、戈壁、盐碱地、沼泽61、62、63、64沼泽有时被归到水域这张表不是让你背而是让你在写统计脚本时有个对照依据。编码体系一旦搞错后面所有面积数字都是空中楼阁。2.2 年份口径和分辨率为什么同一块地两年面积对不上这套数据通常按年份发布常见的有 1980 年代末、1990 年代、2000 年、2005 年、2010 年、2015 年、2020 年等节点。每个年份的生产方式、影像源、分类标准可能有调整。比如早期版本依赖 Landsat TM后期加入国产高分影像分类精度和最小图斑面积都会变。分辨率方面常见的是 1km 栅格和 30m 栅格两种。1km 数据适合全国尺度趋势分析30m 数据适合省域或流域尺度。如果你拿 1km 数据去算一个县的耕地变化边界混合像元会严重干扰结果。我一般会先确认三件事投影坐标系是什么、像元大小是多少、NoData 值填的是什么。这三项不确认后面做掩膜和重采样必出问题。提示不同年份的数据不要直接相减做变化检测先统一投影、重采样到同一网格再对齐分类编码。3. 用 Python 把 rar 解压、读取、裁切到研究区一套能复现的最小流程3.1 解压和目录结构梳理别用鼠标一个个点拿到 .rar 文件Windows 上很多人习惯用 WinRAR 手动解压。但如果你的研究区涉及多个年份、多个分幅手动操作迟早出错。我一般用命令行工具批量处理。Linux 和 macOS 下用 unrarWindows 下可以用 7z 命令行。# 查看压缩包内文件列表先不急着解压 unrar l 中国土地利用现状遥感监测数据.rar # 解压到指定目录保留目录结构 unrar x 中国土地利用现状遥感监测数据.rar ./landuse_data/ # 如果压缩包分卷确保所有分卷在同一目录 # unrar x 中国土地利用现状遥感监测数据.part1.rar ./landuse_data/解压后常见目录结构是按年份分文件夹每个年份下再按省份或分幅存放 GeoTIFF 或 IMG 文件。先别急着写代码用find或dir把文件清单导出来确认命名规律。# 导出文件清单方便后续批量处理 find ./landuse_data -name *.tif -o -name *.img file_list.txt wc -l file_list.txt逻辑说明unrar l只列出内容不写磁盘适合先确认压缩包完整性。unrar x保留目录结构避免文件散落一地。find同时匹配 tif 和 img是因为不同批次数据格式可能不统一。参数上如果你的 rar 有密码加-p参数但这类公开数据一般没有。3.2 用 rasterio 读取并裁切到研究区边界读取栅格我首选 rasterio比 GDAL 原生绑定更友好和 numpy 配合也顺。下面这段代码完成三件事读取土地利用栅格、读取研究区矢量边界、按边界裁切并统计各类型面积。import rasterio from rasterio.mask import mask import geopandas as gpd import numpy as np from collections import Counter # 1. 读取研究区边界 study_area gpd.read_file(./boundary/study_area.shp) # 确保坐标系一致这里假设边界是 WGS84 study_area study_area.to_crs(EPSG:4326) # 2. 打开土地利用栅格 with rasterio.open(./landuse_data/2020/landuse_2020.tif) as src: # 打印元数据确认投影、像元大小、NoData print(src.crs, src.res, src.nodata) # 3. 按边界裁切 geometries study_area.geometry.values out_image, out_transform mask(src, geometries, cropTrue) out_meta src.meta.copy() out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) # 4. 统计各类型像元数 data out_image[0] # 排除 NoData valid data[data ! src.nodata] counts Counter(valid.flatten()) for code, cnt in sorted(counts.items()): print(f编码 {int(code)}: {cnt} 个像元) # 5. 像元数转面积假设像元 30m pixel_area 30 * 30 # 平方米 for code, cnt in sorted(counts.items()): area_km2 cnt * pixel_area / 1e6 print(f编码 {int(code)}: {area_km2:.2f} 平方公里)逻辑说明mask函数按矢量边界裁切栅格cropTrue会缩小输出范围减少内存占用。src.nodata必须读取否则 NoData 会被当成有效值参与统计。像元面积换算时如果数据是地理坐标系经纬度30m 不是固定值需要投影到等面积坐标系再算。参数上study_area.to_crs要和栅格 CRS 一致不一致会裁切出空数组。注意如果研究区跨多个分幅先做镶嵌再裁切不要逐幅裁切后合并否则边界像元统计会重复。3.3 批量处理多年份数据并导出面积表单年份跑通后把流程封装成函数批量处理所有年份。下面这段代码遍历年份文件夹输出一张 CSV 面积表。import os import pandas as pd def stats_by_year(tif_path, boundary, nodataNone): with rasterio.open(tif_path) as src: if nodata is None: nodata src.nodata out_image, _ mask(src, boundary.geometry.values, cropTrue) data out_image[0] valid data[data ! nodata] counts Counter(valid.flatten()) pixel_area src.res[0] * src.res[1] records [] for code, cnt in counts.items(): records.append({ code: int(code), pixel_count: cnt, area_km2: cnt * pixel_area / 1e6 }) return pd.DataFrame(records) years [1990, 2000, 2010, 2020] all_stats [] for y in years: tif f./landuse_data/{y}/landuse_{y}.tif if os.path.exists(tif): df stats_by_year(tif, study_area) df[year] y all_stats.append(df) result pd.concat(all_stats, ignore_indexTrue) result.to_csv(./output/landuse_area_by_year.csv, indexFalse) print(result.head())逻辑说明把单年份逻辑封装成函数方便复用。src.res返回像元宽高相乘得像元面积。如果数据是地理坐标系res单位是度不能直接乘必须先投影。参数上nodata允许外部传入因为有些文件元数据里 NoData 标记不规范需要手动指定比如 0 或 255。4. 投影、重采样和边界对齐三个最容易让面积统计翻车的操作4.1 地理坐标系直接算面积为什么你的平方公里数偏小很多人拿到栅格看元数据里写着 WGS84就直接用像元大小乘个数算面积。这是血泪经验里最常见的一条。经纬度坐标系下一个像元在地表的东西方向距离随纬度变化赤道附近约 30m到了北纬 40 度只剩约 23m。你直接按 30m 算高纬度地区面积会偏大低纬度地区偏小全国尺度上误差能到百分之十几。正确做法是先投影到等面积坐标系。国内常用 Albers 等面积投影中央经线一般取 105°E双标准纬线取 25°N 和 47°N。用 rasterio 或 gdalwarp 都可以。# 用 gdalwarp 投影到 Albers 等面积投影 gdalwarp -t_srs projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm no_defs \ -tr 30 30 -r near \ -srcnodata 255 -dstnodata 255 \ ./landuse_data/2020/landuse_2020.tif \ ./landuse_data/2020/landuse_2020_albers.tif逻辑说明-t_srs指定目标投影-tr 30 30设置输出像元 30m-r near表示最近邻重采样分类数据必须用最近邻不能用双线性否则会出现不存在的类别编码。-srcnodata和-dstnodata保持一致避免 NoData 被重采样成有效值。4.2 重采样时分类数据只能用最近邻双线性会造出假类别分类栅格和连续栅格的重采样逻辑完全不同。NDVI 这种连续值可以用双线性或立方卷积但土地利用编码是离散整数双线性会把 11 和 12 平均成 11.5这个编码在分类体系里根本不存在。我见过有人重采样后统计出编码 15、25 这种奇怪值就是插值惹的祸。最近邻重采样的缺点是边界像元会偏移但至少类别是真实的。如果你需要更精细的边界应该用更高分辨率数据而不是靠插值。重采样方法适用数据类型对分类数据的影响最近邻分类、整数编码保持原始编码边界略偏移双线性连续值、NDVI产生不存在的编码值立方卷积连续值、影像同上且更平滑众数分类数据保持编码但计算量大4.3 多期数据网格对齐像元不重合导致变化检测全错做变化检测时两期数据的网格必须完全对齐。如果一期是 30m另一期是 1km或者投影参数差一点像元边界就错开了。你直接相减会得到大量伪变化。正确做法是以某一期为基准网格把其他期重采样到同一网格。from rasterio.warp import reproject, Resampling # 以 2020 年为基准网格 with rasterio.open(./landuse_data/2020/landuse_2020_albers.tif) as ref: ref_meta ref.meta.copy() ref_data ref.read(1) # 把 2010 年数据重采样到基准网格 with rasterio.open(./landuse_data/2010/landuse_2010_albers.tif) as src: dst np.empty_like(ref_data) reproject( sourcerasterio.band(src, 1), destinationdst, src_transformsrc.transform, src_crssrc.crs, dst_transformref_meta[transform], dst_crsref_meta[crs], resamplingResampling.nearest ) # 保存对齐后的数据 ref_meta.update({dtype: uint8, nodata: 255}) with rasterio.open(./landuse_data/2010/landuse_2010_aligned.tif, w, **ref_meta) as dst_file: dst_file.write(dst, 1)逻辑说明reproject把源数据映射到目标网格Resampling.nearest保证分类编码不变。dst_transform和dst_crs来自基准数据这样两期数据每个像元位置一一对应。参数上dst数组形状必须和基准数据一致否则写入会报错。5. 避坑与排查土地利用遥感监测数据最常见的 5 个翻车现场5.1 现象统计出的建设用地面积逐年下降原因不同年份的建设用地编码口径变了。早期版本农村建设用地可能归在耕地附属后期单独编码。或者你用的分类对照表只对上了一级类二级类编码错位。解决逐年核对分类说明建立年份-编码映射表。不要用一套编码表跑所有年份。如果找不到说明用频率分布对比看哪些编码在某年突然出现或消失。5.2 现象裁切后研究区边界外出现大量有效值原因矢量边界和栅格坐标系不一致或者边界有自相交、空洞。mask函数默认按几何范围裁切如果几何不闭合会裁出多余区域。解决先study_area study_area.buffer(0)修复几何再to_crs对齐坐标系。裁切后用边界再做一次精确掩膜把外部像元设为 NoData。5.3 现象面积统计结果比官方公布数据大很多原因像元面积算错最常见的是地理坐标系直接乘或者像元大小读的是经纬度度数而不是米。另一个原因是 NoData 没排除被当成未利用地统计了。解决确认投影单位投影到等面积坐标系后再算。检查src.nodata如果元数据里是 None手动指定常见值如 0、255、-9999逐个测试。5.4 现象变化检测结果里出现大量 11 变 11 的伪变化原因两期数据网格没对齐或者分类编码版本不同但编码值相同。比如 2010 年的 11 是水田2020 年的 11 可能变成了旱地。解决先做网格对齐再核对两期分类体系。如果编码含义变了先重编码到统一体系再做变化检测。变化检测前用np.where排除 NoData 和相同值。5.5 现象大文件读取时内存溢出原因全国 30m 数据单年份可能几十 GB直接src.read()全读进内存会爆。解决用窗口读取分块处理。rasterio 支持block_windows按块统计后累加。或者先用gdalwarp裁切到研究区再读取小文件。# 分块统计示例 with rasterio.open(tif_path) as src: counts Counter() for ji, window in src.block_windows(1): block src.read(1, windowwindow) valid block[block ! src.nodata] counts.update(valid.flatten())逻辑说明block_windows按栅格内部块遍历每次只读一块内存占用可控。counts.update累加各块统计结果。参数上块大小由数据创建时决定一般 256x256 或 512x512。6. 从面积统计到空间格局一个我常用的验证技巧面积统计只是第一步真正让这份数据发挥价值的是空间格局分析。我一般会做一个简单但有效的交叉验证把统计出的各类型面积和同区域、同年份的公开统计年鉴数据对比。不是要求完全一致而是看量级和趋势是否合理。如果耕地面积差了一个数量级那肯定是投影或编码出了问题。进阶用法上我习惯把土地利用数据和 DEM、NDVI、夜间灯光数据叠加做地形梯度分析或城市扩张分析。比如用 DEM 提取坡度统计不同坡度带上的耕地分布能快速识别陡坡开垦区域。下面这段代码演示如何按坡度带统计耕地像元。import rasterio import numpy as np from rasterio.mask import mask # 读取土地利用和 DEM with rasterio.open(./landuse_albers.tif) as lu_src: lu_data, lu_transform mask(lu_src, study_area.geometry.values, cropTrue) lu lu_data[0] with rasterio.open(./dem_albers.tif) as dem_src: dem_data, _ mask(dem_src, study_area.geometry.values, cropTrue) dem dem_data[0] # 计算坡度这里简化用梯度近似 dy, dx np.gradient(dem) slope np.degrees(np.arctan(np.sqrt(dx**2 dy**2))) # 定义坡度带 slope_bins [0, 2, 6, 15, 25, 90] slope_labels [0-2, 2-6, 6-15, 15-25, 25] # 统计耕地编码 11、12在各坡度带的分布 farmland_mask np.isin(lu, [11, 12]) for i in range(len(slope_bins)-1): mask_bin (slope slope_bins[i]) (slope slope_bins[i1]) count np.sum(farmland_mask mask_bin) print(f坡度 {slope_labels[i]} 度: {count} 个耕地像元)逻辑说明np.gradient计算 DEM 在 x 和 y 方向的梯度arctan转成角度得到坡度。np.isin筛选耕地编码。参数上坡度带划分依据不同研究目的调整做水土流失评估时常用 15 度和 25 度作为阈值。这个技巧的好处是你不需要复杂的模型就能从土地利用数据里挖出有生态意义的信息。我做了这么多年最深的教训就是别一上来就追求花哨的算法先把投影、编码、NoData 这三件事做对比什么都强。希望帮到你。本文还有配套的精品资源点击获取