尧图网页设计 
简介本资源为中国土地利用现状遥感监测数据合集面向地理信息、遥感解译、国土空间规划及生态环境研究方向的科研人员与高校师生可用于长时序土地利用变化分析、区域生态评估与制图建模。数据覆盖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」这个压缩包的人十有八九不是想欣赏数据而是想尽快把某块区域的土地利用类型统计出来。这类数据通常来自全国尺度的遥感解译成果按年份、按省份或按流域分幅存放解压后常见的是栅格影像GeoTIFF 或 GRID 格式加一份分类编码表。它解决的核心问题是不用自己从原始卫星影像开始做分类直接拿到已经解译好的土地利用类型图做面积统计、变化检测、景观格局分析。适合做国土空间规划、生态评估、耕地保护研究以及需要快速出图出数的工程场景。但很多人卡在第一步——解压后一堆文件夹不知道从哪下手或者用 ArcGIS 打开发现坐标系对不上、分类值看不懂。这篇就按我实际处理这类数据的顺序把解压、校验、裁剪、统计、出图整条链路讲清楚参数怎么设、坑在哪都落到具体操作上。2. 解压后先别急着打开目录结构与数据校验2.1 典型目录长什么样这类全国尺度的土地利用遥感监测数据解压后一般不会是一个大文件而是按行政区划或图幅分目录。常见结构是根目录下按省份拼音或行政区代码建文件夹每个省份文件夹里再按年份或期次放栅格文件。栅格文件命名通常包含区域代码、年份、土地利用分类标识。分类体系多数采用二级分类一级类如耕地、林地、草地、水域、建设用地、未利用地二级类在一级类下再细分比如耕地分水田、旱地。分类编码表可能是单独的 Excel 或 DBF 文件也可能直接写在栅格属性表里。先做一件事把目录树打印出来确认文件数量和命名规律。不要用鼠标一个个点直接命令行。# 查看解压后的目录结构只看两层避免输出爆炸 find ./中国土地利用现状遥感监测数据 -maxdepth 2 -type d | sort # 统计栅格文件数量确认是否完整 find ./中国土地利用现状遥感监测数据 -name *.tif -o -name *.img | wc -l # 查看单个栅格文件的基本信息需要 GDAL 环境 gdalinfo ./中国土地利用现状遥感监测数据/某省/某市_2020.tiffind的-maxdepth 2控制目录深度避免全国数据几万个文件夹把终端刷爆。gdalinfo输出里重点看几个字段Size is看行列数Coordinate System is看坐标系Band 1下的Type看数据类型NoData Value看无效值。如果Coordinate System is显示的是地理坐标系经纬度后面做面积统计就要注意了——经纬度下直接算面积会随纬度变化必须投影后再算。2.2 校验三件事坐标系、无效值、分类值范围坐标系不统一是这类数据最常见的翻车点。全国数据可能混用 CGCS2000 地理坐标和 Albers 等面积投影分省数据可能各省自己投影。校验方法批量跑gdalinfo把坐标系信息抓出来对比。# 批量提取所有栅格的坐标系和行列数输出到文本方便比对 for f in $(find ./中国土地利用现状遥感监测数据 -name *.tif); do echo $f gdalinfo $f | grep -E Coordinate System|Size is|NoData Value done crs_check.txt跑完打开crs_check.txt如果发现有的文件是GEOGCS[China Geodetic Coordinate System 2000]有的是PROJCS[Albers_Conic_Equal_Area]那就必须统一。统一的原则做面积统计用等面积投影做位置叠加用地理坐标或统一投影。我一般全部转成 Albers 等面积投影中央经线取 105°E双标准纬线取 25°N 和 47°N这是全国尺度常用的参数。无效值也要看。土地利用分类栅格通常用 0 或 255 表示无效或背景如果NoData Value没设统计时会把 0 当成一个类别算进去面积就错了。分类值范围用gdalinfo -hist看直方图或者用 Python 读一下唯一值。import rasterio import numpy as np # 读取一个栅格查看唯一值和各自像元数 with rasterio.open(./中国土地利用现状遥感监测数据/某省/某市_2020.tif) as src: data src.read(1) unique, counts np.unique(data, return_countsTrue) for u, c in zip(unique, counts): print(f分类值 {u}: {c} 个像元)这段代码用rasterio读第一个波段np.unique返回唯一值和计数。如果发现分类值里有 0 且数量巨大基本就是背景值后续统计要排除。如果分类值出现 100 以上但编码表里没有说明数据可能被重采样过或编码体系不同需要找原始说明确认。3. 裁剪到目标区域按行政区或按自定义范围3.1 用行政区矢量裁剪全国数据动辄几十 GB直接统计全国没必要先裁到目标区域。最常见需求是按省、市、县裁剪。需要准备一份行政区划矢量shp 或 GeoJSON坐标系要和栅格一致不一致先转。# 用 gdalwarp 按矢量裁剪-cutline 指定矢量-crop_to_cutline 裁到边界 gdalwarp -cutline ./boundary/某市.shp \ -crop_to_cutline \ -dstnodata 0 \ -tr 30 30 \ -t_srs EPSG:4526 \ ./中国土地利用现状遥感监测数据/某省/某市_2020.tif \ ./output/某市_2020_clip.tif参数逐个说-cutline是裁剪边界矢量-crop_to_cutline让输出范围紧贴矢量外接矩形-dstnodata 0把裁剪外的区域设为 0方便后续排除-tr 30 30指定输出像元大小 30 米如果原始数据是 30 米就保持一致不要随意重采样-t_srs EPSG:4526指定输出投影这里用的是 CGCS2000 的 3 度带投影具体带号按目标区域经度算。如果矢量坐标系和栅格不一致gdalwarp会自动做投影转换但前提是矢量有正确的.prj文件。裁剪后一定要再跑一次gdalinfo确认输出范围、像元大小、NoData 值都对。我见过有人裁完发现像元大小变成 0.0001 度就是因为没指定-trGDAL 按地理坐标默认输出了。3.2 按自定义范围裁剪有时候目标区域不是行政区比如流域、保护区、项目红线。这时候用-te指定地理范围或者用-cutline传自定义矢量。-te的参数顺序是xmin ymin xmax ymax单位跟-t_srs一致。# 按自定义矩形范围裁剪范围用投影坐标 gdalwarp -te 500000 3200000 560000 3260000 \ -t_srs EPSG:4526 \ -tr 30 30 \ -dstnodata 0 \ ./中国土地利用现状遥感监测数据/某省/某市_2020.tif \ ./output/custom_extent.tif-te的四个值必须落在原数据范围内否则输出可能是空图。不确定范围就先gdalinfo看原数据的Upper Left和Lower Right坐标。自定义范围裁剪适合做固定窗口的对比分析比如同一区域多年份变化检测范围固定了才好逐像元比对。4. 面积统计从像元计数到平方公里4.1 统计逻辑与分类编码映射面积统计的本质是每个分类值有多少个像元乘以单个像元面积得到总面积。单个像元面积 像元宽 × 像元高30 米分辨率就是 900 平方米。但前提是数据已经投影到等面积坐标系否则像元面积不恒定。分类编码映射是另一个关键。栅格里的值可能是 1、2、3但你需要知道 1 代表耕地还是林地。编码表通常在数据说明文档里或者栅格的属性表RAT里。如果都没有只能根据分类体系反推但风险很大。我一般先把唯一值和编码表对一遍确认没有遗漏。import rasterio import numpy as np import pandas as pd # 分类编码映射按实际数据说明修改 code_map { 1: 耕地, 2: 林地, 3: 草地, 4: 水域, 5: 建设用地, 6: 未利用地 } # 像元面积30 米分辨率下为 900 平方米 pixel_area 30 * 30 with rasterio.open(./output/某市_2020_clip.tif) as src: data src.read(1) nodata src.nodata # 排除 NoData if nodata is not None: data np.where(data nodata, 0, data) unique, counts np.unique(data, return_countsTrue) records [] for u, c in zip(unique, counts): if u 0: continue # 跳过背景值 area_km2 c * pixel_area / 1e6 # 平方米转平方公里 records.append({ 分类值: u, 地类: code_map.get(u, 未知), 像元数: c, 面积_平方公里: round(area_km2, 2) }) df pd.DataFrame(records) df.to_csv(./output/某市_2020_面积统计.csv, indexFalse, encodingutf-8-sig) print(df)code_map必须按实际编码表改不能照抄。pixel_area如果数据不是 30 米按实际分辨率算。nodata从src.nodata读如果没设就手动指定。输出 CSV 用utf-8-sig编码Excel 打开不乱码。统计结果里如果出现「未知」地类说明编码表不全需要回去查数据说明。4.2 用 GDAL 命令行快速统计不想写 Python 的话gdalinfo -hist能出直方图但不够直观。更直接的是用gdal_sieve或gdallocationinfo不过批量统计还是 Python 方便。如果只是快速看一个文件可以用gdalinfo -stats然后看Band 1的统计信息但那个不分类别。# 快速查看栅格直方图每个分类值的像元数 gdalinfo -hist -json ./output/某市_2020_clip.tif | \ python -c import sys,json; djson.load(sys.stdin); print(d[bands][0][histogram])这条命令把gdalinfo的 JSON 输出管道给 Python直接打印直方图数组。数组下标就是分类值值是像元数。适合快速验证但不适合出正式报表。5. 避坑与排查坐标系、NoData、分类值错位5.1 面积算出来偏大或偏小现象统计出的耕地面积比官方公布数据大 20% 以上。原因栅格还是地理坐标系像元面积按固定值算但经纬度下像元实际面积随纬度变化低纬度地区偏小高纬度地区偏大。解决统计前先gdalwarp -t_srs转到等面积投影Albers 或 Lambert 都可以参数按区域选。5.2 裁剪后边缘出现大量 0 值现象裁剪后的图边缘一圈 0统计时如果没排除 NoData0 被当成一个地类。原因gdalwarp裁剪时默认用 0 填充边界外区域如果原数据 NoData 不是 0就冲突了。解决裁剪时显式指定-dstnodata值和原数据 NoData 一致或者统计时手动排除 0。5.3 分类值对不上编码表现象栅格里出现 11、12、21 这样的值但编码表只有 1 到 6。原因数据可能用了二级分类编码比如 11 代表水田、12 代表旱地编码表只给了一级类。解决找二级分类编码表或者按一级类合并。合并逻辑11、12 归为耕地21、22 归为林地以此类推。5.4 多个年份数据叠加时范围不一致现象做变化检测时两年数据裁剪范围差几个像元逐像元比对报错。原因不同年份数据可能用了不同的图幅分幅或者裁剪时-te范围没对齐。解决以其中一年为基准用gdalwarp把其他年份重采样到相同的范围和分辨率-te和-tr都显式指定。5.5 大文件处理内存爆掉现象全国数据直接读进内存做统计Python 进程被 kill。原因全国 30 米栅格几十 GB一次性read()内存扛不住。解决分块读取用rasterio的block_windows或者gdal的Translate分块处理。或者先按省裁剪再逐省统计。# 分块读取统计避免内存爆掉 with rasterio.open(./中国土地利用现状遥感监测数据/全国_2020.tif) as src: counts {} for ji, window in src.block_windows(1): block src.read(1, windowwindow) unique, c np.unique(block, return_countsTrue) for u, cnt in zip(unique, c): counts[u] counts.get(u, 0) cnt print(counts)block_windows按栅格内部块遍历每次只读一块内存占用可控。统计结果累加到字典里最后统一算面积。6. 进阶变化检测与出图技巧变化检测是这类数据的高频进阶用法。核心逻辑两年分类图逐像元比对生成转移矩阵。转移矩阵的行是前一期地类列是后一期地类值是从行地类转到列地类的像元数或面积。这个矩阵能直接看出耕地转建设用地、林地转耕地等关键变化。import rasterio import numpy as np import pandas as pd with rasterio.open(./output/某市_2010_clip.tif) as src1: data1 src1.read(1) with rasterio.open(./output/某市_2020_clip.tif) as src2: data2 src2.read(1) # 确保两期数据形状一致 assert data1.shape data2.shape, 两期数据行列数不一致先对齐 # 生成转移矩阵 mask (data1 0) (data2 0) transition np.zeros((7, 7), dtypenp.int64) # 假设分类值 1-6 for i in range(data1.shape[0]): for j in range(data1.shape[1]): if mask[i, j]: transition[data1[i, j], data2[i, j]] 1 df pd.DataFrame(transition, index[f2010_{i} for i in range(7)], columns[f2020_{i} for i in range(7)]) df.to_csv(./output/转移矩阵.csv, encodingutf-8-sig) print(df)逐像元循环在 Python 里很慢全国数据不要这么写。实际工程中我用numpy的向量化操作或者pandas的crosstab。这里为了逻辑清晰用了循环小区域可以跑大区域换成pd.crosstab(data1.flatten(), data2.flatten())。出图方面土地利用图的关键是配色。分类图不要用连续色带要用定性配色每个地类一个固定颜色。QGIS 里可以加载栅格后手动设色或者用gdal_translate转成调色板 PNG。出图前记得加图例、比例尺、指北针这些在 QGIS 打印布局里做。我一般会把统计出的面积表直接贴在图的角落一张图把空间分布和数量结构都说清楚。最后说个血泪经验这类数据拿到手先花半小时把坐标系、NoData、分类编码这三件事确认清楚再动手裁剪统计。我见过太多人跳过校验直接跑结果面积算错、图对不上回头返工的时间够把校验做十遍。数据说明文档如果找不到就去翻压缩包里的readme或者说明.txt实在没有就按唯一值和分类体系反推但一定要在报告里注明推断依据。希望帮到你。本文还有配套的精品资源点击获取 中国土地利用遥感监测数据实战：从解压到面积统计与变化检测
简介本资源为中国土地利用现状遥感监测数据合集面向地理信息、遥感解译、国土空间规划及生态环境研究方向的科研人员与高校师生可用于长时序土地利用变化分析、区域生态评估与制图建模。数据覆盖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」这个压缩包的人十有八九不是想欣赏数据而是想尽快把某块区域的土地利用类型统计出来。这类数据通常来自全国尺度的遥感解译成果按年份、按省份或按流域分幅存放解压后常见的是栅格影像GeoTIFF 或 GRID 格式加一份分类编码表。它解决的核心问题是不用自己从原始卫星影像开始做分类直接拿到已经解译好的土地利用类型图做面积统计、变化检测、景观格局分析。适合做国土空间规划、生态评估、耕地保护研究以及需要快速出图出数的工程场景。但很多人卡在第一步——解压后一堆文件夹不知道从哪下手或者用 ArcGIS 打开发现坐标系对不上、分类值看不懂。这篇就按我实际处理这类数据的顺序把解压、校验、裁剪、统计、出图整条链路讲清楚参数怎么设、坑在哪都落到具体操作上。2. 解压后先别急着打开目录结构与数据校验2.1 典型目录长什么样这类全国尺度的土地利用遥感监测数据解压后一般不会是一个大文件而是按行政区划或图幅分目录。常见结构是根目录下按省份拼音或行政区代码建文件夹每个省份文件夹里再按年份或期次放栅格文件。栅格文件命名通常包含区域代码、年份、土地利用分类标识。分类体系多数采用二级分类一级类如耕地、林地、草地、水域、建设用地、未利用地二级类在一级类下再细分比如耕地分水田、旱地。分类编码表可能是单独的 Excel 或 DBF 文件也可能直接写在栅格属性表里。先做一件事把目录树打印出来确认文件数量和命名规律。不要用鼠标一个个点直接命令行。# 查看解压后的目录结构只看两层避免输出爆炸 find ./中国土地利用现状遥感监测数据 -maxdepth 2 -type d | sort # 统计栅格文件数量确认是否完整 find ./中国土地利用现状遥感监测数据 -name *.tif -o -name *.img | wc -l # 查看单个栅格文件的基本信息需要 GDAL 环境 gdalinfo ./中国土地利用现状遥感监测数据/某省/某市_2020.tiffind的-maxdepth 2控制目录深度避免全国数据几万个文件夹把终端刷爆。gdalinfo输出里重点看几个字段Size is看行列数Coordinate System is看坐标系Band 1下的Type看数据类型NoData Value看无效值。如果Coordinate System is显示的是地理坐标系经纬度后面做面积统计就要注意了——经纬度下直接算面积会随纬度变化必须投影后再算。2.2 校验三件事坐标系、无效值、分类值范围坐标系不统一是这类数据最常见的翻车点。全国数据可能混用 CGCS2000 地理坐标和 Albers 等面积投影分省数据可能各省自己投影。校验方法批量跑gdalinfo把坐标系信息抓出来对比。# 批量提取所有栅格的坐标系和行列数输出到文本方便比对 for f in $(find ./中国土地利用现状遥感监测数据 -name *.tif); do echo $f gdalinfo $f | grep -E Coordinate System|Size is|NoData Value done crs_check.txt跑完打开crs_check.txt如果发现有的文件是GEOGCS[China Geodetic Coordinate System 2000]有的是PROJCS[Albers_Conic_Equal_Area]那就必须统一。统一的原则做面积统计用等面积投影做位置叠加用地理坐标或统一投影。我一般全部转成 Albers 等面积投影中央经线取 105°E双标准纬线取 25°N 和 47°N这是全国尺度常用的参数。无效值也要看。土地利用分类栅格通常用 0 或 255 表示无效或背景如果NoData Value没设统计时会把 0 当成一个类别算进去面积就错了。分类值范围用gdalinfo -hist看直方图或者用 Python 读一下唯一值。import rasterio import numpy as np # 读取一个栅格查看唯一值和各自像元数 with rasterio.open(./中国土地利用现状遥感监测数据/某省/某市_2020.tif) as src: data src.read(1) unique, counts np.unique(data, return_countsTrue) for u, c in zip(unique, counts): print(f分类值 {u}: {c} 个像元)这段代码用rasterio读第一个波段np.unique返回唯一值和计数。如果发现分类值里有 0 且数量巨大基本就是背景值后续统计要排除。如果分类值出现 100 以上但编码表里没有说明数据可能被重采样过或编码体系不同需要找原始说明确认。3. 裁剪到目标区域按行政区或按自定义范围3.1 用行政区矢量裁剪全国数据动辄几十 GB直接统计全国没必要先裁到目标区域。最常见需求是按省、市、县裁剪。需要准备一份行政区划矢量shp 或 GeoJSON坐标系要和栅格一致不一致先转。# 用 gdalwarp 按矢量裁剪-cutline 指定矢量-crop_to_cutline 裁到边界 gdalwarp -cutline ./boundary/某市.shp \ -crop_to_cutline \ -dstnodata 0 \ -tr 30 30 \ -t_srs EPSG:4526 \ ./中国土地利用现状遥感监测数据/某省/某市_2020.tif \ ./output/某市_2020_clip.tif参数逐个说-cutline是裁剪边界矢量-crop_to_cutline让输出范围紧贴矢量外接矩形-dstnodata 0把裁剪外的区域设为 0方便后续排除-tr 30 30指定输出像元大小 30 米如果原始数据是 30 米就保持一致不要随意重采样-t_srs EPSG:4526指定输出投影这里用的是 CGCS2000 的 3 度带投影具体带号按目标区域经度算。如果矢量坐标系和栅格不一致gdalwarp会自动做投影转换但前提是矢量有正确的.prj文件。裁剪后一定要再跑一次gdalinfo确认输出范围、像元大小、NoData 值都对。我见过有人裁完发现像元大小变成 0.0001 度就是因为没指定-trGDAL 按地理坐标默认输出了。3.2 按自定义范围裁剪有时候目标区域不是行政区比如流域、保护区、项目红线。这时候用-te指定地理范围或者用-cutline传自定义矢量。-te的参数顺序是xmin ymin xmax ymax单位跟-t_srs一致。# 按自定义矩形范围裁剪范围用投影坐标 gdalwarp -te 500000 3200000 560000 3260000 \ -t_srs EPSG:4526 \ -tr 30 30 \ -dstnodata 0 \ ./中国土地利用现状遥感监测数据/某省/某市_2020.tif \ ./output/custom_extent.tif-te的四个值必须落在原数据范围内否则输出可能是空图。不确定范围就先gdalinfo看原数据的Upper Left和Lower Right坐标。自定义范围裁剪适合做固定窗口的对比分析比如同一区域多年份变化检测范围固定了才好逐像元比对。4. 面积统计从像元计数到平方公里4.1 统计逻辑与分类编码映射面积统计的本质是每个分类值有多少个像元乘以单个像元面积得到总面积。单个像元面积 像元宽 × 像元高30 米分辨率就是 900 平方米。但前提是数据已经投影到等面积坐标系否则像元面积不恒定。分类编码映射是另一个关键。栅格里的值可能是 1、2、3但你需要知道 1 代表耕地还是林地。编码表通常在数据说明文档里或者栅格的属性表RAT里。如果都没有只能根据分类体系反推但风险很大。我一般先把唯一值和编码表对一遍确认没有遗漏。import rasterio import numpy as np import pandas as pd # 分类编码映射按实际数据说明修改 code_map { 1: 耕地, 2: 林地, 3: 草地, 4: 水域, 5: 建设用地, 6: 未利用地 } # 像元面积30 米分辨率下为 900 平方米 pixel_area 30 * 30 with rasterio.open(./output/某市_2020_clip.tif) as src: data src.read(1) nodata src.nodata # 排除 NoData if nodata is not None: data np.where(data nodata, 0, data) unique, counts np.unique(data, return_countsTrue) records [] for u, c in zip(unique, counts): if u 0: continue # 跳过背景值 area_km2 c * pixel_area / 1e6 # 平方米转平方公里 records.append({ 分类值: u, 地类: code_map.get(u, 未知), 像元数: c, 面积_平方公里: round(area_km2, 2) }) df pd.DataFrame(records) df.to_csv(./output/某市_2020_面积统计.csv, indexFalse, encodingutf-8-sig) print(df)code_map必须按实际编码表改不能照抄。pixel_area如果数据不是 30 米按实际分辨率算。nodata从src.nodata读如果没设就手动指定。输出 CSV 用utf-8-sig编码Excel 打开不乱码。统计结果里如果出现「未知」地类说明编码表不全需要回去查数据说明。4.2 用 GDAL 命令行快速统计不想写 Python 的话gdalinfo -hist能出直方图但不够直观。更直接的是用gdal_sieve或gdallocationinfo不过批量统计还是 Python 方便。如果只是快速看一个文件可以用gdalinfo -stats然后看Band 1的统计信息但那个不分类别。# 快速查看栅格直方图每个分类值的像元数 gdalinfo -hist -json ./output/某市_2020_clip.tif | \ python -c import sys,json; djson.load(sys.stdin); print(d[bands][0][histogram])这条命令把gdalinfo的 JSON 输出管道给 Python直接打印直方图数组。数组下标就是分类值值是像元数。适合快速验证但不适合出正式报表。5. 避坑与排查坐标系、NoData、分类值错位5.1 面积算出来偏大或偏小现象统计出的耕地面积比官方公布数据大 20% 以上。原因栅格还是地理坐标系像元面积按固定值算但经纬度下像元实际面积随纬度变化低纬度地区偏小高纬度地区偏大。解决统计前先gdalwarp -t_srs转到等面积投影Albers 或 Lambert 都可以参数按区域选。5.2 裁剪后边缘出现大量 0 值现象裁剪后的图边缘一圈 0统计时如果没排除 NoData0 被当成一个地类。原因gdalwarp裁剪时默认用 0 填充边界外区域如果原数据 NoData 不是 0就冲突了。解决裁剪时显式指定-dstnodata值和原数据 NoData 一致或者统计时手动排除 0。5.3 分类值对不上编码表现象栅格里出现 11、12、21 这样的值但编码表只有 1 到 6。原因数据可能用了二级分类编码比如 11 代表水田、12 代表旱地编码表只给了一级类。解决找二级分类编码表或者按一级类合并。合并逻辑11、12 归为耕地21、22 归为林地以此类推。5.4 多个年份数据叠加时范围不一致现象做变化检测时两年数据裁剪范围差几个像元逐像元比对报错。原因不同年份数据可能用了不同的图幅分幅或者裁剪时-te范围没对齐。解决以其中一年为基准用gdalwarp把其他年份重采样到相同的范围和分辨率-te和-tr都显式指定。5.5 大文件处理内存爆掉现象全国数据直接读进内存做统计Python 进程被 kill。原因全国 30 米栅格几十 GB一次性read()内存扛不住。解决分块读取用rasterio的block_windows或者gdal的Translate分块处理。或者先按省裁剪再逐省统计。# 分块读取统计避免内存爆掉 with rasterio.open(./中国土地利用现状遥感监测数据/全国_2020.tif) as src: counts {} for ji, window in src.block_windows(1): block src.read(1, windowwindow) unique, c np.unique(block, return_countsTrue) for u, cnt in zip(unique, c): counts[u] counts.get(u, 0) cnt print(counts)block_windows按栅格内部块遍历每次只读一块内存占用可控。统计结果累加到字典里最后统一算面积。6. 进阶变化检测与出图技巧变化检测是这类数据的高频进阶用法。核心逻辑两年分类图逐像元比对生成转移矩阵。转移矩阵的行是前一期地类列是后一期地类值是从行地类转到列地类的像元数或面积。这个矩阵能直接看出耕地转建设用地、林地转耕地等关键变化。import rasterio import numpy as np import pandas as pd with rasterio.open(./output/某市_2010_clip.tif) as src1: data1 src1.read(1) with rasterio.open(./output/某市_2020_clip.tif) as src2: data2 src2.read(1) # 确保两期数据形状一致 assert data1.shape data2.shape, 两期数据行列数不一致先对齐 # 生成转移矩阵 mask (data1 0) (data2 0) transition np.zeros((7, 7), dtypenp.int64) # 假设分类值 1-6 for i in range(data1.shape[0]): for j in range(data1.shape[1]): if mask[i, j]: transition[data1[i, j], data2[i, j]] 1 df pd.DataFrame(transition, index[f2010_{i} for i in range(7)], columns[f2020_{i} for i in range(7)]) df.to_csv(./output/转移矩阵.csv, encodingutf-8-sig) print(df)逐像元循环在 Python 里很慢全国数据不要这么写。实际工程中我用numpy的向量化操作或者pandas的crosstab。这里为了逻辑清晰用了循环小区域可以跑大区域换成pd.crosstab(data1.flatten(), data2.flatten())。出图方面土地利用图的关键是配色。分类图不要用连续色带要用定性配色每个地类一个固定颜色。QGIS 里可以加载栅格后手动设色或者用gdal_translate转成调色板 PNG。出图前记得加图例、比例尺、指北针这些在 QGIS 打印布局里做。我一般会把统计出的面积表直接贴在图的角落一张图把空间分布和数量结构都说清楚。最后说个血泪经验这类数据拿到手先花半小时把坐标系、NoData、分类编码这三件事确认清楚再动手裁剪统计。我见过太多人跳过校验直接跑结果面积算错、图对不上回头返工的时间够把校验做十遍。数据说明文档如果找不到就去翻压缩包里的readme或者说明.txt实在没有就按唯一值和分类体系反推但一定要在报告里注明推断依据。希望帮到你。本文还有配套的精品资源点击获取