这里插一句,HWSD2 全称是 Harmonized World Soil Database version 2.0,也就是协调世界土壤数据库第二版。不少人第一次打开它都能看见一堆 .mdb 后缀的数据库文件和一张全球 ASCII 栅格图,但真正要用它做研究、做分析、出一张漂亮的全球土壤属性分布图时,就会立刻卡在同一个地方:数据是“分类代码”,不是“属性值”。这篇文章就是专门解决这个问题的,从数据解析、字段连接、栅格化出图,到踩坑排查,完整捋一遍思路和操作。
1. HWSD2 数据构成与前期准备
1.1 数据集里到底装了什么
HWSD2 是 FAO 和 IIASA 联合维护的全球土壤数据产品,比起第一版,它最大的进步是把土壤分层信息做得更细致,同时统一了不同国家土壤分类系统之间的差异。整个数据集大体上由两部分组成:
第一部分是全球栅格文件,通常以 ASCII Grid 格式提供,全图每个栅格像元都对应一个“土壤单元编号”。这里需要注意,栅格的值并不是你想用的“有机碳含量”或“pH值”,而是一个指向数据库的 ID 号。这个 ID 号才是数据集的灵魂,所有的属性值都要靠它去关联。
第二部分是属性数据库,常见的格式有 Microsoft Access 的 .mdb 文件,也有直接提供 CSV 或 SQLite 版本的情况。数据库里每一条记录都对应一个土壤单元,记录里包含了非常多的属性字段,比如有机碳、pH、容重、阳离子交换量、沙粒黏粒含量,甚至还有参考深度和石砾含量。
我刚拿到 HWSD2 的时候犯过一个错误,以为通过 ArcGIS 直接打开栅格属性表,把各个像元的值改成图层符号就行。后来发现那只是“分类显示”,任何一个想用在统计、建模、制图里的真实数值都没法从栅格属性表里直接读出来。所以操作的第一个关键动作,就是要明白:栅格负责定位,数据库负责属性,两者之间靠 ID 连接。
1.2 为什么不能直接拿原始数据出图
很多人会把 HWSD2 的紫红色全球预览图误以为是一张“有机碳分布图”,或者“土壤类型图”。其实那个紫色预览图只是整个原始栅格的默认符号化结果,代表的是图层库里每个土壤单元的索引值在渐变显示。你要想做成“全球土壤有机碳密度图”,至少得走下面这几步:
- 读取原始栅格中每个像元的 MU_GLOBAL(或 ID)值;
- 在属性数据库中筛选出一个目标属性字段,比如 T_OC(顶层有机碳含量);
- 把属性值通过连接赋值给栅格对应的每个像元;
- 重新输出为另一张栅格,这时每个像元值才是你真正需要的属性数值。
如果只停留在“打开栅格、符号化一下”的层面,导出的图片放到论文里一定会被审稿人或同行质疑。因为数据的空间分布、数值范围、缺测处理这些信息完全没有体现出来,人家根本不知道你到底有没有把数据库关联成功。
所以我的经验是,在做任何提取和栅格化之前,先把 HWSD2 自带的数据字典或字段说明表完整过一遍。字段名称虽然带一定规律,比如 T_OC 表示表层有机碳,S_OC 表示下层有机碳,但不同版本里对深度的定义并不一致,有的版本 T 表示 0-30 cm,有的版本改成了 0-20 cm 或其它划分,一定不能靠猜。
2. 属性表关键字段解读
2.1 顶层与底层字段的差异
HWSD2 属性库里几乎每个土壤单元都有两套土壤理化性质,分别对应顶层土壤和底层土壤。顶层和底层的划分深度在不同数据版本里略有不一样,但通常默认是 0-30 cm 和 30-100 cm。字段名称里常以“T_”开头代表顶层,以“S_”开头代表底层。
我在做项目时经常只需要表层数据,比如全球表层土壤有机碳制图。这时如果不去区分 T 和 S,直接把全部字段导出来,生成的栅格在数值上会有非常大的分层跳变,尤其是底层土壤有机碳普遍低于表层,如果混在一起,图面会出现一块一块的“断层感”。
建议提取前先确认字段含义。常见字段大致有这些:
| 字段名 | 含义 | 单位 | 常见取值范围 |
|---|---|---|---|
| T_OC / S_OC | 有机碳含量 | % | 0-30 |
| T_PH_H2O / S_PH_H2O | 土壤 pH(水浸) | 无单位 | 4-9 |
| T_CEC_SOIL / S_CEC_SOIL | 阳离子交换量 | cmol/kg | 0-100 |
| T_BULK_DENSITY / S_BULK_DENSITY | 土壤容重 | kg/dm3 | 0.8-1.8 |
| T_SAND / S_SAND | 砂粒含量 | % | 0-100 |
| T_SILT / S_SILT | 粉粒含量 | % | 0-100 |
| T_CLAY / S_CLAY | 黏粒含量 | % | 0-100 |
| T_GRAVEL / S_GRAVEL | 石砾含量 | % | 0-100 |
| T_REF_DEPTH / S_REF_DEPTH | 参考深度 | cm | 0-300 |
我一般建议先把这些字段在 Excel 里做一轮统计,看看最大值、最小值和缺失值比例。有些栅格像元没有对应的数据库记录,连接之后就会变成 NoData。如果缺失太多,要么考虑用邻域插值填补,要么在制图说明里写明不包含某些地区。
2.2 MU_GLOBAL 与数据库主键的对应关系
HWSD2 原始栅格的值其实在不同版本里有不同的字段名,常见的是 MU_GLOBAL,也有的版本直接叫 ID,或者叫 MU。它本质上是土壤制图单元的全局唯一编号,和属性数据库里的主键是一一对应的。
我在做连接操作时,最常遇到的问题是数据库里主键有重复。因为不同国家提交的数据质量参差不齐,某些土壤单元可能在库里出现了多条记录,有的记录来自旧版,有的来自修正数据。如果直接做一对多连接,栅格属性表或矢量属性表里会突然多出很多重复行,后续栅格化时会出现一个像元叠加多个值的情况,软件通常默认用第一个匹配值,但你根本不知道它是哪一条。
做法是连接之前先对数据库按主键做去重,或者用 GIS 工具进行一对一连接。如果不确定重复情况,可以用 Excel 的删除重复项功能快速筛查。数据库里如果有多条记录对应同一个 MU_GLOBAL,就要看数据版本里是否有“优先级”字段,没有的话最好用最新修订时间的记录。
还有一点容易被忽略:ASCII 栅格里的整型值在导入部分软件时可能发生精度变化,比如值 3200012 会被读成 3200010 或科学计数法显示。这种问题在连接环节特别致命,明明看着是两个相同的数字,软件却说无法匹配。解决方案是用长整型读取原始栅格属性表,并在连接前将两侧字段统一转为长整型格式。
3. 实战操作:从属性连接到全球栅格出图
3.1 方案 A:ArcGIS 里的“栅格转面 + 属性连接”
如果你习惯 ArcGIS 的图形界面操作,这条路最直接:先把全球栅格转成矢量面,然后把属性数据库表连接到矢量面属性表,最后按目标字段把矢量面重新转成栅格。
在 ArcToolbox 里找到“转换工具”下的“从栅格转面”工具,输入 HWSD2 原始栅格。这里有个关键参数:是否简化面。如果勾选了简化面,软件会合并相邻且值相同的面,同时清除一些锯齿边。从出图角度说,全球尺度下勾选简化会让图面更干净,但处理时间更长;从数据精度角度说,不简化可以保留每一个原始像元边界。
面生成后,打开属性表会看到每条记录都有 GRIDCODE 字段,这个值就是栅格像元的值,也就是 MU_GLOBAL。然后右键图层,选择连接和关联中的“连接”,把 GRIDCODE 与数据库主键连接起来。连接时一定要在“连接字段”里精确选择,不要靠软件自动匹配,否则可能出现大量空值。
连接成功以后,可以直接使用“面转栅格”工具,输入要素是已连接的面要素,选择需要输出的字段,比如 T_OC,输出像元大小最好设置为与原始栅格一致,通常是 30 弧秒左右。这样输出的栅格值就是真实土壤属性值,可以用来做后续的裁剪、统计和制图。
这套方案的优点是操作直观,适合不常写代码的人;缺点是全球尺度的数据面数量非常大,转面后文件动辄几个 GB,普通笔记本运行起来风扇会像飞机起飞一样。实测下来,8 GB 内存跑全局转面大概率会卡死,建议先按大洲或区域裁剪后再处理,或者用下面这种更轻量的方法。
3.2 方案 B:Python + GDAL/rasterio 高效批处理
用代码处理 HWSD2 其实没有想象中复杂,它的核心就三步:读栅格 DN 值,建立 DN 到属性的字典映射,替换并输出新栅格。即使你不经常写 Python,照着框架改字段名也能跑通。
以下是我测试过的可直接运行的思路:
import rasterio import pandas as pd import numpy as np # 1. 读取原始栅格的像元值 src_raster = r"HWSD2_global.asc" attrs_csv = r"HWSD2_attributes.csv" with rasterio.open(src_raster) as src: data = src.read(1).astype(np.float64) profile = src.profile # 拿到栅格元数据,之后输出时原样使用 # 2. 读取属性表,这里假设你已把mdb转成了csv df = pd.read_csv(attrs_csv, encoding="utf-8-sig") df = df.drop_duplicates(subset=["MU_GLOBAL"]) # 3. 建立映射字典:MU_GLOBAL -> T_OC mapping = dict(zip(df["MU_GLOBAL"], df["T_OC"])) # 4. 用向量化映射替换原始值 new_data = np.vectorize(lambda x: mapping.get(int(x), np.nan))(data) # 5. 写出新栅格 profile.update(dtype=rasterio.float64, nodata=np.nan) with rasterio.open(r"TOC_global.tif", "w", **profile) as dst: dst.write(new_data.astype(rasterio.float64), 1)这段代码的核心优势在于没有转面过程,直接一步到位。全球尺度的栅格可能有两亿个像元,但本方法的耗时主要卡在读取和写入,一般不会超过 20 分钟,普通电脑都能跑。
需要说明一点:上面代码里“假设你已把 mdb 转成了 csv”,这个步骤很多人会卡住。最简单的转换方式是使用 Microsoft Access 打开后另存,或者用 Python 的 pandas 加 pyodbc 直接读取 Access 数据库表。如果不想折腾,直接用 ArcGIS 的“表转表”工具,把 Access 表导出成 DBF 或 CSV 也行。
我还建议把映射字典的结果做一次回检:随机抽取 20 个像元,对比输出栅格值和数据库源表数值是否一致。这个方法看上去笨,但能够快速发现字段映射错位、数据类型未对齐之类的问题,比什么都强。
3.3 如何在 ArcGIS/ENVI 中快速检查结果
在完成栅格提取之后,不要直接拿出去用,先做一个基本目视检查。我说一下快速的验证流程:
打开结果栅格属性表,看最小值、最大值是否落在合理区间。比如有机碳的合理范围通常在 0-30% 之间;如果看到负值、超过 100 的怪值,说明字段连接有问题,或者数据库里混了其它非土壤属性值。
接下来检查 NoData 的空间分布。很多人在全球图里都会发现不少区域的空值,这是正常的,因为湖泊、冰川、城市和一些没有调查数据的地区本身就没有土壤数据库条记录。但如果 NoData 区域呈现规整的矩形、条纹状,则很可能是坐标系偏移或字段匹配失败造成的。
在 ENVI 里可以使用波段运算工具把 NoData 统一替换为指定的背景值,方便出图。例如想在图上把无数据区域显示为灰色而不是透明,用波段运算写一句(b1 lt 0)*255 + (b1 ge 0)*b1就能实现,具体数值按你自己底图的范围调整。
ArcGIS 的栅格计算器里,判断栅格是否为空值用 IsNull 函数,清理空值则用 SetNull 函数嵌套条件;如果你更习惯工具面板,可以直接检索“设置空值”工具,输入条件表达式与原栅格即可。这些操作本身不难,但新手的误区常常是把“空值”和“数值 0”混为一谈,导致出图时大片区域显示成黑色。
4. 常见问题与排查技巧实录
4.1 栅格转面后出现大量空值怎么办
这是我在处理 HWSD2 时遇到最多的坑。栅格转面后,面要素属性表里的 GRIDCODE 字段时常出现大面积的空值,不是真正的空值,而是 2147483647 这类整型最大值,也可能是 -9999。这些值对应数据集里的背景值、水域或冰盖像元。
如果直接在面属性表里做连接,这些背景值会连接不上数据库,于是变成大片空白面。解决方法是先在 ArcToolbox 里用条件分析或栅格计算器,把特定背景值区域设置为 NoData,再执行栅格转面。表达式类似:
SetNull("HWSD2_raw" == -9999, "HWSD2_raw")这样转出来的面就不会包含那些无效区域,后续连接效率更高,图面也更干净。在实际项目中,我通常会把全球水域和冰川单独处理,另存为掩膜,方便后用。
4.2 属性连接后字段无法正常显示
连接成功后,在面属性表里可以看到目标字段,但制作成栅格后符号化始终是一片黑,或者根本没有分类。这种情况多半是面转栅格时字段类型选错了。比如有机碳字段原本是浮点型,但你转栅格时选择了输出整型,软件会对数值取整,0.23% 变成 0%,结果自然是黑屏。
解决办法是手动指定输出栅格像素类型为浮点型。ArcGIS 的“面转栅格”工具里并没有直接的像素类型参数,但你可以在环境设置里把“输出栅格像素类型”设为 FLOAT。Python 方案则直接在 profile 里写dtype=rasterio.float32。
还有另一层原因:连接表用的是 Excel 文件时,字段名超过 10 个字符可能被截断,导致面属性表找不到你要的字段。比如 T_REF_DEPTH 改名成了 T_REF_DEP,不是原来的字段名,连接时手一滑就选错了。用 CSV 或 DBF 文件做连接更稳妥。
4.3 坐标系不一致导致错位
HWSD2 原始数据通常采用地理坐标系,以弧秒为单位,但部分环境在读取 ASCII Grid 时默认给了 WGS84 之外的定义,或者和另一个数据叠加时自动套用了已有图层的投影,结果右下角、左上角位置出现明显偏移或旋转。
遇到类似情况时,先检查输出栅格的投影信息是否与原始栅格一致。如果只是在出图时与其它图层叠加,可以临时在 ArcGIS 的图层属性里设置投影,但不建议修改原始文件。最严谨的做法是在数据处理阶段就保持统一的坐标系,比如使用 WGS84 地理坐标系,到了制图阶段再根据出图范围动态投影到 Albers 等面积投影。
4.4 全球文件太大,设备卡顿怎么办
8 GB 内存的电脑处理全球尺度的 HWSD2 转面确实吃力。给一个亲测可行的建议:不要一次处理全球,使用全球陆地矢量边界裁剪出目标区域,再在区域范围内做转面和属性连接。制图时如果论文研究区是中国、非洲或某个流域,这样做完全够用,性能还能提升数倍。
Python 方案也可以做分块读取。rasterio 提供了窗口读取功能,可以先按经度分成六个分块,每块处理完写回目标栅格。这样即使数据总量有 2 亿像元,单次内存峰值也不会太高。
还有一个小技巧:如果只是提取某一属性做全球地图,没必要保存中间面文件,直接从编码后的栅格用查表索引输出结果就行。这样硬盘只多一个结果栅格,完全绕开“面”这个中间产物的容量压力。
5. 制图技巧与成果扩展
5.1 用重分类与渐变符号化提升出图效果
属性栅格生成后,默认的灰白色渐变很难体现土壤属性的空间规律。通常我会对栅格做分段重分类,再给每一个区间赋予特定的配色。比如全球表层土壤有机碳,用量级分类区间 0-1%、1-2%、2-4%、4-8%、8-12%、12% 以上,用浅黄到深棕的渐变,图面直观且不夸张。
在 ArcGIS 里可以用“重分类”工具操作。需要注意的一点是,重分类会改变栅格值,如果后续还要做面积统计或建模计算,建议保留原始浮点栅格,另存一个重分类结果用于制图展示。千万不要拿符号化的图去做后续提取,数值早就不是真实含量了。
对于只想做简单的数据预览,可以在符号系统里选择“分类”并手动设置断点,这不会改变数值。如果要发布共享或者投稿,再用重分类结果出图更好。
5.2 从单属性栅格扩展到多属性产品
HWSD2 的价值不仅限于制作某一个属性的栅格。你可以把有机碳、pH、砂粒、黏粒、容重这些关键字段都按照相同的流程批量生成,从而得到一个“全球土壤关键属性栅格集”。
有了这个基础数据后,能做的事就很多:
- 估算区域土壤有机碳储量,单位面积土壤碳密度乘以不同深度的土层厚度,再加上容重校正;
- 做土壤适耕性评价,把 pH、有机质、质地等级代入评价模型;
- 与气象、地形数据叠加,用统计模型制作土壤属性空间预测图;
- 在流域模型里作为下垫面土壤参数输入。
我自己的实操习惯是:把提取脚本写成函数,输入属性字段名,自动生成对应的结果栅格,然后统一输出到一个文件夹。这样以后项目需要哪几个属性,修改一个列表参数,几分钟就能跑完全部结果。建议你也早点打破“手工一个个操作”的思维,代码化以后,重复工作会轻松十倍。
5.3 空间连接时如何避免重复值污染结果
前面提到数据库主键重复的问题,还有一个容易被忽略的切入点:属性数据库里有些记录是“组合土壤单元”,代表一个区域内混合了两种或更多土壤类型。这种记录在 HWSD2 里经常出现,字段里会用“份额”或“比例”表示每种土壤的覆盖面积。
如果这种情况影响到了你提取的准确性,常用的处理思路有两种。第一种是直接忽略,适用于全球制图这种宏观尺度,因为组合单元本身的属性值是按面积加权的综合值。第二种是按比例拆分,用面积权重重新分配属性到子单元,适用于精度要求较高的区域研究。
从操作实现上说,第二种方式比较复杂,需要把原始栅格的每个像元面积计算出来,结合比例字段拆分。一般不建议新手起步就直接做,等基础流程熟练后再尝试也不迟。
6. 个人实操心得与建议
HWSD2 这套数据我前后用过不下二十次,最深刻的体会是:数据集本身并不复杂,真正决定成果质量的是你对属性表字段含义的掌握以及每一步之间是否做了充分验证。很多人卡在“转面”“连接”“转栅格”这些步骤里反复重试,大部分原因不是软件不熟,而是没理解栅格 DN 值和属性数据库之间的指代关系。
如果你准备在自己的项目里用这套流程,建议按下面的次序来推进:
- 先把原始栅格属性表打开,确认 DN 值的范围,比如在 0-30000 之间,判断是否存在负值或奇异的背景值;
- 再打开属性数据库,确认唯一 ID 数量和栅格非背景像元值个数是否基本一致;
- 随后抽几个样本点,从原始栅格读取位置值,再从数据库里反向查属性,两边人工核对;
- 全部无误后再批量生成目标属性栅格。
还有一个细节想提醒:HWSD2 的属性值代表的是“土壤单元的属性”,而不是严格意义上的“每个像元的实测值”。在制图和统计时,建议在成果说明里明确标注这一点,否则容易被人误读为高分辨率实测栅格。数据本身是制图综合的产物,空间分辨率名义上是 30 弧秒,实际有效精度往往低于这个水平,用的时候心里要有数。
最后再分享一个很小的技巧:做全球栅格时,饼图、柱状图这些统计图表放在小比例尺下往往看不清,反而用密度分割加连续色带更能表达空间趋势。如果你还要叠加地形阴影底图,记得把土壤栅格做一点透明度处理,比如设成 70% 透明度,这样既能看到地形起伏背景,又能清楚读出土壤属性的分布规律,最终出图效果会专业很多。