简介中国30米分辨率土壤类型栅格数据面向GIS、土壤与生态环境研究者可用于土壤资源调查、农业区划及生态评估等场景。数据分为东三省与全国不含东三省、港澳台两套Soil_Raster栅格并配套土壤类型代码表便于属性查阅与制图分析。压缩包共15个文件以tif栅格、tfw坐标配准、xml元数据、ovr金字塔及dbf属性表为主另有xlsx类型对照表整体249.54MB。已有130人学习下载。拿到后可直接在ArcGIS或QGIS中加载结合土壤代码表快速提取、重分类与专题制图节省数据预处理时间是一份实用且可直接入库的基础地理数据。1. 一张30m的土壤栅格能省掉你半个野外季做区域生态评估或农业区划时最尴尬的不是缺数据而是手里的数据尺度对不上需求。国家级的土壤图常见是1:100万或1km分辨率放在省级尺度还能看个轮廓一旦切到县或流域河道两侧的土壤变化、山前平原与低山丘陵的过渡带全都糊成一片。而这份30m分辨率的中国土壤类型数据像素大小相当于把一个县切成上百万个30m方格每个格子都带一个土壤类型编码配合附带的土壤类型代码表你可以直接按代码筛选、重分类、叠加分析不需要再从纸质土壤图或扫描PDF里一点点数字化。适合谁用它做土地利用变化的遥感分类、跑SWAT或InVEST模型、算区域土壤有机碳储量的人都能把它作为基础底图。但拿到手后不能直接当黑匣子用分类体系是发生分类还是系统分类、坐标系是WGS84还是Albers、代码表和栅格像素值是否完全对应这几个问题不先核对清楚后边的每个分析步骤都会跟着翻车。下面按我处理这类栅格数据的习惯顺序从数据结构讲到实际落地再讲踩过的坑。2. 30m土壤类型数据的内部组成栅格、属性表和代码表之间的关系2.1 先分清三样东西tif栅格、属性表和你要用的代码表这类型数据在交付时通常是一个标准GeoTIFF文件加一份Excel或CSV代码表。栅格文件里每个像素的取值是一个整数这个整数本身没有意义只有对照代码表才能翻译成「棕壤」「红壤」「水稻土」这样的类型名称。听起来很简单但实际使用时有一个容易被忽略的层次问题代码表里往往不止一层对应关系。常见做法是代码表至少包含四列代码值、土壤类型中文名、英文名或拼音缩写以及所属土纲或土类。以中国土壤发生分类为例代码可能编排到土类级别也可能到亚类级别。如果是亚类级别一个土类下边会挂着多个亚类这时你如果按土类做重分类就需要先在代码表里做一次分组汇总不能直接在栅格上操作。我习惯先把代码表导入到DataFrame里检查「高级分类单元」这一列的去重数量以此判断这份数据到底细到什么程度。import pandas as pd # codebook_path换成实际路径编码常见为gbk或utf-8 df pd.read_excel(/path/to/soil_codebook.xlsx, engineopenpyxl) # 看一眼列名 print(df.columns.tolist()) # 检查代码是否唯一 print(代码唯一性:, df[code].is_unique) # 如果存在父级分类列统计父类数量 if parent_group in df.columns: print(父类数量:, df[parent_group].nunique()) print(df.groupby(parent_group)[code].count())第一段代码的功能是确认代码表里的代码值是否有重复第二段是看高级分类单元的数量。如果代码不唯一说明栅格里同一个像素值可能对应多个含义需要立即联系数据提供方确认如果父类数量明显偏少而代码行数很多说明这份表细分到了亚类级别。这两点决定了后续所有筛选和统计的粒度。2.2 先判断分类体系发生分类和系统分类不能混用中国土壤分类体系有两套并行一套是传统的发生分类强调成土过程和地带性常见名称如暗棕壤、黑钙土、栗钙土、红壤、黄壤等另一套是系统分类参考国际土壤分类体系用诊断层和诊断特性来定义名称形式是「淋溶土」「雏形土」「人为土」等。30m土壤栅格数据在使用前最要紧的一件事就是确认代码表里的类型名称属于哪一套体系。有一个现象值得注意数据本身分辨率做到30m编制年代越近的产品越倾向于使用系统分类。原因在于系统分类是用定量指标划分的便于在计算机里做空间推演而旧的发生分类更多依赖剖面形态描述数字化之后容易出现边界争议。但也正因为如此很多人拿到的数据可能是「新瓶装旧酒」——栅格的原始信息来自第二次土壤普查的发生分类图后处理时把类型名称翻译成了系统分类代码表里两套名字并存。import rasterio import numpy as np with rasterio.open(/path/to/soil_type_30m.tif) as src: print(CRS:, src.crs) print(像元大小:, src.res) print(数据范围:, src.bounds) # 读取全图或按窗口读取先看像素值的取值范围 data src.read(1) vals np.unique(data) print(像素值数量:, vals.size) print(取值分布(前50个):, vals[:50])这段代码是拿到栅格后我必跑的第一步。CRS决定你后面要不要做重投影像元大小如果是0.000277度相当于经纬度坐标那么在地图上的实际地面宽度在不同纬度不一样如果像元大小是30或30.xx说明原数据已经做过投影。像素值的取值范围则和代码表直接对照如果最大值超过代码表里的行数大概率是栅格内混入了背景值或NoData需要单独掩膜。2.3 属性表RAT不是必须但有时能救急很多GeoTIFF会在内部嵌一个栅格属性表Raster Attribute Table记录像素值和对应类型名称的映射。这套机制在ArcGIS里打开图例时能直接看到中文名但在QGIS里不一定显示在GDAL读取时需要用专门的接口读取。我处理过一份数据代码表文件是损坏的反而是栅格内嵌的RAT帮了大忙。from osgeo import gdal ds gdal.Open(/path/to/soil_type_30m.tif) rat ds.GetRasterBand(1).GetDefaultRAT() if rat is not None: row_count rat.GetRowCount() print(RAT行数:, row_count) col_count rat.GetColumnCount() for i in range(min(row_count, 5)): # 第0列通常是像素值其余列可能是土壤名称等字段 row {} for j in range(col_count): row[rat.GetNameOfCol(j)] rat.GetValueAsString(i, j) print(row)有一点需要注意RAT里的第0列才是像素值本身其他列是描述字段。名字可能叫ClassName也可能叫SoilName或CATEGORY不同数据源的命名完全不同。读取RAT的目的不是替代代码表而是做交叉验证——如果代码表和RAT对同一像素值的中文名不一致说明数据在分发前经过一次重分类两类名称之间可能存在转换误差。3. 数据落地实操Python和桌面GIS两种路径处理30m土壤栅格3.1 Python最小可跑流程裁剪、重投影和代码表关联能做到什么程度拿到数据后最常用的三个操作按行政区裁剪、统一坐标系、像素级提取统计。我用rasterio加geopandas组合来完成原因是这两个库都可以在内存里处理不需要额外起ArcGIS环境。下面给出一套完整的最小流程从读取行政边界到产出裁剪后的栅格和一张按土壤类型统计的面积表。import geopandas as gpd import rasterio from rasterio.mask import mask from rasterio.warp import calculate_default_transform, reproject, Resampling import numpy as np import pandas as pd # 1. 读行政区边界 boundary gpd.read_file(/path/to/boundary.gpkg) # 确保栅格和矢量的坐标系一致不一致时先统一到栅格的CRS if boundary.crs ! rasterio.open(/path/to/soil_type_30m.tif).crs: boundary boundary.to_crs(rasterio.open(/path/to/soil_type_30m.tif).crs) # 2. 按边界裁剪 with rasterio.open(/path/to/soil_type_30m.tif) as src: out_image, out_transform mask(src, boundary.geometry, cropTrue, nodata0) out_meta src.meta.copy() out_meta.update({ driver: GTiff, height: out_image.shape[1], width: out_image.shape[2], transform: out_transform, nodata: 0 }) with rasterio.open(/path/to/soil_type_clip.tif, w, **out_meta) as dst: dst.write(out_image) # 3. 裁剪结果统计像素值计数 面积计算 with rasterio.open(/path/to/soil_type_clip.tif) as clipped: data clipped.read(1) nodata clipped.nodata pixel_area clipped.res[0] * clipped.res[1] # 平面坐标系下才成立 valid_mask data ! nodata values, counts np.unique(data[valid_mask], return_countsTrue) stat_df pd.DataFrame({code: values, pixel_count: counts}) # 30m分辨率栅格每个像素面积约为0.09公顷平面投影下 stat_df[area_ha] stat_df[pixel_count] * pixel_area / 10000这套代码里有三个参数值得注意。第一是mask函数里的nodata参数如果不显式指定栅格原文件的nodata值可能不是0裁剪后背景值会和实际土壤类型混淆统计出的每个类型面积都会偏大。第二是计算面积时用的pixel_area只在投影坐标系下才准确如果栅格是经纬度坐标需要用像素宽度乘纬度余弦做修正或者干脆先重投影到Albers等面积投影再算。第三是输出文件的nodata统一设成0这个值不能和土壤类型编码表里的任何合法编码冲突否则后续筛选时容易误伤。如果栅格是经纬度坐标系重投影这步不能省。常见目标投影是Albers等积圆锥投影Krasovsky_1940_Albers参数上中央经线、双标准纬线要和目标区域匹配。但重投影会带来一个隐患重采样会改变像素值。对土壤类型栅格这种类别数据必须使用最近邻法。# 4. 重投影到Albers等积投影示例参数按数据范围微调 dst_crs projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm no_defs with rasterio.open(/path/to/soil_type_clip.tif) as src: transform, width, height calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds, resolution30 ) kwargs src.meta.copy() kwargs.update({ crs: dst_crs, transform: transform, width: width, height: height, nodata: 0 }) with rasterio.open(/path/to/soil_type_albers.tif, w, **kwargs) as dst: for i in range(1, src.count 1): reproject( sourcerasterio.band(src, i), destinationrasterio.band(dst, i), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingResampling.nearest )重投影里最常见的问题不是选错坐标系而是选了双线性或三次卷积插值。土壤类型是类别变量0到1之间没有过渡意义插值会凭空产生不存在的编码值。我见过有人用双线性法重采样后栅格里出现一批新的整数代码表里根本查不到整份数据因此报废。Resampling.nearest保证输出值永远来自输入值集合。3.2 桌面GIS处理路径ArcGIS和QGIS里的关键参数位如果你不写Python桌面GIS也能完成同样工作但有几个参数位容易看漏。ArcGIS里推荐流程是先用Project Raster工具重投影Environment里的Resampling Technique选择NEAREST然后用Clip工具按面裁剪Environment里的Snapping Raster选原始土壤栅格避免输出范围出现错位最后用Reclassify或Join Field把代码表关联到栅格属性表。QGIS里流程更简单Raster菜单下的WarpReproject选「Nearest Neighbour」然后Clipper裁剪最后用「Add Vector Fields in Attribute Table」把代码表CSV以连接方式关联进来。但QGIS在栅格和属性表关联上默认只支持数值字段如果代码表里的土壤名称是中文关联后能显示但导出会被截断需要检查输出编码。桌面路径的优点是每一步能看到中间结果缺点是无法留痕参数记录靠截图或手写笔记。如果团队协作建议至少把关键参数存成QGIS模型或ArcGIS Toolbox脚本方便复现。3.3 分辨率基准30m栅格和矢量图斑边界在什么情况下会冲突30m分辨率的特殊之处在于它介于Landsat遥感解译和传统土壤普查图之间。第二次土壤普查的原始图件很多是1:50万甚至1:100万比例尺数字化后转成30m栅格图斑边界会呈现明显的锯齿状。这并不代表数据精度达到30m只是存储单元变小了。你在做叠加分析时如果另一个数据是1:100万比例尺的矢量边界两者叠加出的细碎多边形都是假象。我一般会在数据分析前做一个降尺度测试把30m栅格聚合成250m或1km看土壤类型的面积比例变化是否超过10%。如果聚合后比例变化很小说明原数据在空间上比较平滑如果变化剧烈说明数据里有大量孤立像素这些位置存在分类噪声。翻车现场通常是没做这个测试就直接进入模型参数提取结果每个图斑里都带了一堆异常点。4. 避坑指南拿到30m土壤数据后容易翻车的五个典型场景4.1 现象裁剪后边界处多出一圈带状异常像素原因行政边界矢量的坐标系和栅格不一致或者矢量本身带有非常细碎的边界点。ArcGIS裁剪时如果没开Snapping Raster输出范围会和栅格像元边界错开半个像元导致边界处出现重采样产生的像素。解决裁剪前统一两套数据的坐标系ArcGIS里在Environment栏指定Snapping RasterPython流程用rasterio.mask并传入boundary.geometry输出transform会自动对齐边界但保险起见再用reproject对齐到原始栅格的网格。4.2 现象代码表和栅格像素值能对上但土壤类型面积比例和文献差异巨大原因大概率是分类体系判断错误。例如代码表里写的是「简育湿润雏形土」这类系统分类名称你却按发生分类的「棕壤」去找参考值自然对不上。文献里公开的统计数据往往基于发生分类面积汇总。解决先确定代码表的体系归属再看文献中统计口径是否一致。如果代码表混用了两套体系发生分类土类名、系统分类土纲名并存每一行都要单独甄别不能按列名粗分。我的习惯是在代码表里加一列「体系标签」逐行标注后再汇总。4.3 现象重采样后栅格里出现代码表中不存在的新值原因重投影或重采样时用了双线性、三次卷积等插值法。土壤类型是类别变量插值产生的中间值没有实际意义但ArcGIS默认重采样方法往往不是最近邻QGIS里Warp对话框的默认选项也容易被人忽略。解决把Project Raster和Warp的Resampling Type强制改为NEAREST。如果错误已经发生唯一补救办法是重新下载原始数据再跑一遍不要尝试用栅格计算器修数值。4.4 现象面积统计结果比官方数据多出百分之十几原因栅格中把NoData背景值编码为0而代码表里0可能没有对应类型统计时忘了排除。更有一种隐蔽场景离岛、南海诸岛在图幅范围外被填充成某个有效编码但官方统计不含这些岛屿。解决先读取栅格stats确认nodata值再用代码表全表比对栅格内的全部取值最后在统计脚本里显式排除nodata。南海诸岛的问题建议用国家发布的标准行政区边界裁剪后再统计。4.5 现象同一像素位置在ArcGIS和QGIS里显示的类型名称不一样原因软件对栅格属性表的读取方式不同。ArcGIS优先读取内部RATQGIS则显示的是栅格原始数值或颜色映射表。如果你在QGIS里看到的是数字自行对照代码表时把数字看成了另一列编码就会产生错觉。解决把代码表和RAT都导出来做三向核对——栅格像素值、RAT名称、代码表名称。三处一致才能动手分析。5. 进阶自检把30m土壤类型数据真正用进模型前的三道验证工序模型对输入数据质量的要求比制图展示高一个量级。制图时一个像素的误差看不出来模型参数错一格碳储量或径流系数可能差出一个数量级。所以我额外做三道验证。第一道是空间一致性抽查。选一个县或小流域把30m土壤栅格叠加到1:5万或1:10万地形图上检查土壤类型边界是否和地貌界线吻合。具体做法是用县内典型剖面点的土壤类型记录做验证点逐个用rasterio提取像素值比对代码表名称。我遇到过的问题是山地暗棕壤与河谷草甸土的边界偏移超过两个像素这类偏差通常源于原始图件数字化时的配准误差无法通过后处理修正只能在模型参数分配时保留它的不确定性。第二道是聚合稳定性测试。把30m栅格分别聚合到90m、250m、1km计算每个类型的面积比例画一条面积随分辨率变化的曲线。如果曲线在某个尺度上出现断裂说明该尺度下类型图斑破碎度异常模型的网格尺寸最好不要小于这个尺度。对于30m数据土壤图斑的实际最小面积通常远大于一个像素如果聚合到250m后面积比例变化超过5%一般建议模型网格取250m而不是30m。from rasterio.enums import Resampling def aggregate_ratio(input_tif, scale_factor): with rasterio.open(input_tif) as src: # scale_factor3表示聚合成90m30m*390m data src.read( 1, out_shape(1, src.height // scale_factor, src.width // scale_factor), resamplingResampling.mode ) vals, counts np.unique(data, return_countsTrue) total counts.sum() ratios {int(v): c / total for v, c in zip(vals, counts)} return ratios for factor in [1, 3, 8, 33]: # 30m, 90m, 240m, 990m ratios aggregate_ratio(/path/to/soil_type_albers.tif, factor) print(f聚合尺度{factor}:, ratios)注意这段代码用Resampling.mode进行聚合对类别栅格而言mode众数比majority多数更稳因为mode对并列情况有确定处理逻辑。计算出的面积比例可以导成表后续写进模型文档。第三道是代码表语义终查。这里的重点是确认代码表里每一行的中文类型名在国家标准分类中有明确对应不能有「其他」「未分类」这类占位项。占位项在制图时无碍但在模型里如果土壤属性数据库查不到这类土种的物理参数整个像元的模拟结果就是空的。我会在代码表末尾加一列「模型属性匹配状态」把有明确土壤理化参数的类型标为OK查不到的类型标为MANUAL后续再为MANUAL类型指定平均参数或最邻近土壤的属性。这三道工序做完这份30m数据才算真正从「能看」变成「能用」。我早期做碳储量估算时跳过第二道聚合测试直接用30m像元统计结果一个县域的土壤有机碳总量比同行的估算高出40%排查了很久发现是大量孤立小图斑中的有机质异常值没有被尺度效应平滑掉。从那以后每次拿到新的土壤栅格数据都会先跑一遍这三道检验希望这套流程也能帮你在项目里少走这一步弯路。本文还有配套的精品资源点击获取