从去年开始我一直在做全国尺度的农业环境模型最头疼的从来都不是算法本身而是“数据荒”。你要一张全国范围的土壤属性分布图传统土壤图拿出来基本都是第二次土壤普查时期的外业调查成果图斑粗糙属性表里翻来翻去就是质地、pH、有机质那几项往模型里一塞结果自己心里都没底。后面花了大量时间筛选数据源最后锁定了一套中国高分辨率国家土壤信息网格基本属性数据集2010–2018年1km栅格、TIFF格式、16项精细属性覆盖密度、质地、厚度等关键物理指标。今天就把这套数据的使用心得、技术细节和我踩过的坑一次性说清楚。这套数据到底解决了什么问题简单说它把几十年积累的土壤调查点数据、环境协变量和机器学习预测方法揉在一起生成了一张张全国连续分布的土壤属性栅格图空间分辨率1km时间范围是2010–2018年属性维度做到了16项。不管是做土壤碳库估算、作物适种区划分还是跑水文模型、生态模型这套数据完全可以当主力底图不用再自己拿点位数据辛辛苦苦做空间插值。适合的人群很明确GIS方向的研究生、搞土壤调查和制图的专业人员、做环境模拟和农业区划的工程师以及所有在项目里被“土壤数据从哪来”困扰的朋友。1. 数据集为什么会以“网格”和“栅格”的面貌出现1.1 传统土壤图 vs 土壤信息网格过去我们在项目里用到的土壤数据多数是矢量多边形也就是把一块区域按土壤类型或者土种边界切成一堆面每个面挂一条属性记录。这种数据的好处是接近传统调查成果坏处也相当明显边界往往很生硬同一个多边形内的土壤属性被当成完全均匀的实际调查采样点也稀少图斑之间经常出现莫名其妙的跳变。做模型的人拿到这种面数据要么只能给整个多边形赋同一个值要么还要自己做“面转点再插值”的二次加工工作量巨大。土壤信息网格则完全是另一套思路。它把整个研究区域划分成规则的格网每个格网单元像元都独立给出一套属性预测值不再有“同一个多边形内完全相同”这种强制假设。这样得到的空间分布更平滑也更符合连续型土壤属性的实际变化规律。再加上预测模型会引入高程、气候、植被、母岩等协变量网格之间的差异背后是有环境逻辑支撑的不是单纯机械插值。1.2 1km分辨率意味着什么标题里写的“高分辨率”不是自吹在土壤属性制图这个领域1km确实已经属于比较精细的全国尺度产品。你可能觉得1km很粗毕竟现在遥感影像动不动就是10米、30米但土壤属性这个东西和地表覆盖不一样它的空间变异尺度本身就比较大而且没有任何传感器能直接从天上观测到土壤容重、质地或者有机碳含量。要在全国范围内得到连续的1km栅格背后需要大量样点、多源协变量和成熟的预测模型工程量和计算量都相当可观。1km栅格在实际应用中到底对应多大范围如果数据采用经纬度坐标1km近似为0.008333度也就是1/120度赤道附近一个像元大概是1km乘1km如果采用等积投影坐标像元尺寸会直接定义成1000m乘1000m。对全国尺度分析来说这个粒度足够捕捉气候带、地形带和大尺度土壤类型带之间的差异但你要是做某个县、某个乡镇的精确耕作区划1km就显得太糙了建议再用局部的高精度采样数据来修正。1.3 时间范围2010–2018年意味着什么这个时间范围不是随口标的它代表的是数据源中环境协变量和观测站点数据所覆盖的主要时段。土壤属性的变化相对缓慢但也不是一成不变尤其是有机碳、养分含量这些受人类活动影响明显的指标几十年间完全可能发生显著变化。这套数据用2010–2018年的协变量和调查数据做预测比直接用上世纪80年代第二次土壤普查结果要新得多也更能反映当下的土壤状态。不过要提醒一点这套数据本质上是“空间预测产品”不是某个时间点的实测全量普查结果。2010–2018这个范围更准确的理解是“基准期”你把它用在2015年前后的项目里是最合适的如果拿去代表2023年或者2024年的土壤状况理论上会有一定偏差但多数情况下不会影响趋势判断。2. 16项精细属性到底有哪些每一项怎么用2.1 物理属性密度、质地、厚度标题明确点出来的三项里“密度”一般指土壤容重bulk density单位通常是g/cm³数值范围多数在1.0到1.6之间。容重是计算土壤碳储量、水分储量最关键的参数之一模型里经常需要它来把质量分数换算成单位面积的储量没有这个数你后面所有计算都得靠猜。我一般用这套数据里的容重栅格做基础输入实测效果比按土壤类型赋经验值靠谱得多。“质地”相关属性通常不是单一字段而是拆分成了沙粒、粉粒、黏粒三组百分含量有时还会附带基于美国农部制或国际制的土壤质地分类结果。质地直接决定土壤的持水特性、导水特性和抗侵蚀能力水文模型和作物模型里几乎都是核心参数。拿到这三项后你还能用土壤质地三角图自行判断土类也可以用来计算土壤有效水含量。“厚度”指的是土壤剖面厚度常见单位是cm这个值对计算土体中的总碳储量、总养分库容特别重要做区域碳汇评估的时候厚度和容重两个数据必须配合使用。2.2 化学属性有机碳、pH、CEC、养分指标物理属性之外这套数据里还会有一批化学属性。土壤有机碳SOC含量通常以g/kg为单位是土壤质量的核心指标也是碳循环模型最重要的输入参数之一。pH值直接决定土壤酸碱状况做农作物适种区划、重金属有效性评估时离不开它。阳离子交换量CEC反映土壤保肥能力单位多为cmol/kg做土壤肥力评价和污染物迁移模拟时需要。养分指标一般包括全氮、有效磷、速效钾等虽然单位和提取方法在不同来源数据里会有差异但做宏观农业区划和施肥策略初判完全够用。这些化学属性在空间上的分布规律比物理属性更复杂受人类耕作、施肥、土地利用历史的影响很大。使用时要额外留意数值的合理性比如pH在4到9之外、有机碳含量突然高到几十g/kg很可能是有问题的像元需要结合局部实际情况判断。2.3 水文/衍生属性与属性单位坑如果数据版本内容再丰富一点16项里通常还会包含一些水文衍生指标比如田间持水量、萎蔫系数、饱和导水率、砾石含量等。这些指标多数是从质地、有机碳等基础属性通过传递函数pedotransfer function推导出来的虽然有一定不确定性但作为区域尺度模型输入已经很有价值。用这套数据最容易出的事故就是单位。同样是土壤有机碳有的数据用g/kg有的用%有的用kg C/m²表示碳密度差一个数量级是常事。我在实际项目里就遇到过同事直接把g/kg当成%代入模型结果整个区域的碳储量凭空大了10倍。建议拿到数据的第一个动作就是把每个属性文件的单位、深度档位、无数据值全部记录下来做成一张元数据表贴在项目文档里后面不管谁接手都不会出大错。2.4 属性数据从栅格到表格的提取栅格数据本身是“一张图对应一个属性”但很多时候我们需要的是“按行政区或者样地点提取属性值”这时候就要从栅格里做值提取。QGIS里用Point Sampling Tool插件或者工具箱里的“Sample raster values”可以按点提取ArcGIS里用Extract Multi Values To Points一次能导出多个栅格到点要素属性表。更普适的做法是直接用Python的rasterio库配合geopandas批量读取示例代码如下import rasterio import geopandas as gpd import pandas as pd points gpd.read_file(sample_points.shp) rst_path bulk_density_1km.tif with rasterio.open(rst_path) as src: coords [(x, y) for x, y in zip(points.geometry.x, points.geometry.y)] values [v[0] for v in src.sample(coords)] points[bulk_density] values points.to_file(sample_points_with_bd.shp)这段代码的本质是把点的X、Y坐标丢给栅格数据集的sample方法按最邻近或双线性方式取到像元值。注意坐标参考必须和栅格一致不一致的话要先做坐标转换。我自己的习惯是先把点数据统一投影到栅格的坐标系再做提取省得犯空间错位这种低级错误。3. TIFF格式背后的技术细节3.1 GeoTIFF结构与时序这套数据用的是TIFF格式准确说是GeoTIFF也就是在标准TIFF基础上嵌入了地理空间元数据。打开文件时软件能自动读取坐标系、地理范围、像元尺寸无需另配世界文件。如果你拿到一个扩展名是.tif的文件但软件里显示没有投影信息十有八九是元数据被剥离了只剩普通TIFF。GeoTIFF既可以存单波段灰度图也可以存多波段文件。这套数据如果是每个属性一个独立TIFF那就是16个单波段文件如果打包成一个多波段TIFF读取时你会在软件里看到16个波段。两种方式各有好处单波段文件方便单独加载和算数运算多波段文件方便整体管理和批量处理。在使用前务必确认每一层的属性名和单位别让波段顺序弄混。3.2 坐标系、分辨率和像元尺寸换算1km分辨率在不同坐标系下的定义方式不同。如果是经纬度坐标如WGS84或CGCS2000分辨率写的是0.008333333度相当于赤道附近约1km但高纬度地区经度方向上实际距离会缩短所以严格来说全球等经纬度格网只有赤道附近才真正是“正方形”。如果是投影坐标如Albers等积投影、Lambert投影分辨率会直接写1000单位是米。做全国尺度分析时我通常建议先投影到Albers等积投影或兰伯特等角投影再开始重采样和面积统计否则用经纬度算面积会出现较大误差。如果你在ArcGIS里用Project Raster工具做投影转换注意选择合适的重采样方法连续属性用双线性或三次卷积分类属性用最邻近法这一点下面还会细讲。3.3 文件大小与存储优化1km格网覆盖中国陆域约960万平方公里大约对应960万个有效像元。一个Float32类型的单波段TIFF未压缩时约为38MB960万×4字节。16项属性全部展开就是约600MB如果分深度档位存储总量会成倍增加。实际下载时看到的文件可能比这个小因为很多产品做了压缩如LZW或Deflate大幅降低存储占用。不建议因为文件大就转成ESRI Grid或其他私有格式GeoTIFF是跨平台通用性最好的栅格格式QGIS、ArcGIS、R、Python都能无障碍读取。如果实在觉得读取慢可以先用GDAL将数据转成COGCloud Optimized GeoTIFF这样在QGIS里能实现按需读取打开速度会快很多。3.4 QGIS/ArcGIS 加载步骤加载TIFF属于GIS入门操作但几个细节值得注意。QGIS里直接拖拽.tif文件到图层面板或者用“Layer - Add Layer - Add Raster Layer”选择文件。ArcGIS Pro里用“Add Data”拖进去就能看到。加载后第一件事是查看图层属性里的“Source”信息确认坐标系、像元大小、范围和像元类型是否符合预期。如果数据没有自动显示先检查是不是没有数据值导致显示范围全是异常大值再设置透明度或拉伸方式让正常值显示出来。4. 实操环境从加载到制图的完整流程4.1 裁剪研究区并保持栅格对齐大多数人拿到全国数据后第一步就是裁剪到自己的研究区。QGIS里用“Clip Raster by Extent”工具箱搜“clip raster”选择一个矢量边界要素作为裁剪范围输出TIFF。ArcGIS里对应“Clip”工具栅格裁剪时勾选“Use Input Features for Clipping Geometry”可以精确按边界裁剪避免矩形多余区域影响后期统计。但裁剪只是第一步更关键的问题是“网格原点、范围、分辨率逐格对应”。不同栅格图层如果来源不同像元网格往往错位直接做波段计算或叠加分析就会出现“同一个区域使用的像元位置对不齐”的问题。解决这个问题最稳妥的办法是统一重采样到同一个网格模板具体操作是在QGIS用“Raster - Align Raster”或者在GDAL里用gdalwarp先确定一个基准图层然后把其他所有图层按基准图层的分辨率、范围和网格原点重采样。gdalwarp -t_srs EPSG:4326 -tr 0.008333333 0.008333333 -r bilinear -tap -cutline study_area.shp -crop_to_cutline input_soc.tif output_soc_clip.tif这里面有几个参数值得解释。-tr指定像元大小-tap的作用是让输出网格原点与目标网格对齐-r bilinear表示双线性重采样适合连续型变量。如果做的是土壤类型这类分类变量应改成-r near最邻近法否则类别值会被平滑成无意义的中间值。ArcGIS里则可以用“Raster Calculator”配合“Snap Raster”环境设置在环境参数里指定参考栅格所有输出都会自动对齐。4.2 栅格重分类与栅格转面栅格重分类在QGIS里的操作是“Raster Layers - Raster Calculator”或者SAGA/Grass提供的“Reclassify”工具。我常用的是QGIS的“Reclassify by Table”工具可以设置阈值区间映射成新值。举个实际例子如果我要把土壤容重数据从连续值转成“紧实/疏松”两个等级就可以设置1.0–1.3为低密级1.3–1.6为高密级。重分类之后很多人会想把栅格结果转成面要素比如把“土壤质地类型”栅格转成带属性表的多边形。QGIS里“Raster - Conversion - PolygonizeRaster to Vector”就可以做。但这里有个大坑如果对连续属性直接做栅格转面每个不同像元值都会被划成一个独立多边形结果就是研究区里密密麻麻全是碎斑根本没法用。正确做法是先重分类成少量几个类别再转面最后用“Dissolve融合”按类别字段合并。ArcGIS里的“Raster to Polygon”工具也遵循同样的逻辑转面后记得把“在主键字段中融合”选项打开。4.3 矢量图与栅格图对比使用很多人纠结到底是该用矢量图还是栅格图这个问题我在实际工作中总结的经验是做统计分析和建模用栅格做制图和拓扑分析用矢量。栅格的优势是每个像元都有值数学运算和模型模拟非常方便矢量的优势是边界精确、属性表结构清晰、能直接做拓扑查询和空间分析。在QGIS中你可以把栅格作为底图叠加矢量边界查看空间分布再用“Zonal statistics”分区统计工具计算每个行政区内土壤属性的平均值、最大值、最小值。这个工作流比单纯用栅格转矢量更高效因为不需要处理成千上万个微型多边形只需要一个行政区矢量文件和一个属性栅格就够了。4.4 投影栅格报错000151 的排查ArcGIS里用“Project Raster”工具时经常遇到报错代码000151这个问题在QQ群和知乎被反复问过。根据我自己的排查经验绝大多数情况是因为源栅格没有官方定义的投影信息或者投影信息定义错误。ArcGIS无法判断源坐标系自然无法完成转换。解决步骤很简单第一步打开栅格属性查看坐标系如果显示未知或空值先用“Define Projection”工具手工指定正确的坐标系第二步检查栅格实际空间范围是否和定义的坐标系大范围匹配比如明明是经纬度数据你定义成投影坐标范围会差到离谱第三步再执行Project Raster目标坐标系选清楚。一个小技巧如果源栅格范围在0到180、0到60左右多半是经纬度如果范围是几百万、几千万米的数量级那就是投影坐标。5. 我在使用中踩过的坑与排查技巧5.1 深度档位千万别搞混土壤属性数据经常按深度分成多个层次比如0–5cm、5–15cm、15–30cm、30–60cm、60–100cm。这套数据如果按层发布每个深度档位都是一个独立栅格文件。最容易被忽略的是你在报告里写“0–30cm有机碳储量”结果把0–5cm和5–15cm的数据直接相加压根没考虑中间缺了15–30cm。这种错误在自动批处理脚本里特别容易发生因为你可能循环读文件时读错了层。我的做法是下载数据后第一时间建立文件命名规则表把每个文件的深度档位写清楚然后写脚本检查文件名里的深度标识和实际读取到的属性值范围是否一致。宁可慢一步不能错一步。5.2 无数据值处理不当栅格数据里很多区域没有有效的土壤观测或预测值比如水域、永久积雪区、某些极端地形区数据生产方会用无数据值填充常见的有-9999、0、-3.4028235e38等。如果你直接拿这些值参与统计或做栅格计算结果会被严重污染。用Python读取时要显式指定nodata值并把它排除到计算外在QGIS里则可以通过图层属性设置“No Data Value”ArcGIS也在环境设置里处理。我遇到过最无语的情况是用栅格计算器做乘法没有先设置无数据值结果无数据区域变成了0下游面积统计时这些区域被当成“土壤属性为0的裸地”整个出来的体量少了十分之一。所以在一开始就把无数据值处理好这是数据预处理的重中之重。5.3 坐标系混乱导致的对不齐问题使用多来源数据时坐标系混乱几乎是必然的。有的数据是WGS84经纬度有的是CGCS2000有的是Albers还有的是UTM分带。表面上看加载到软件里可能因为动态投影显示位置差不多但一旦你做像元级别的计算就会发现两幅图“错位”了。解决的办法前面提过用gdalwarp统一坐标系和分辨率同时检查坐标范围是否在合理区间比如中国区域的经纬度范围大致是73°E到135°E、18°N到54°N如果你的影像范围出现负数或偏离这个范围坐标系八成有问题。5.4 数据质量检查的小习惯栅格数据再专业也不是“绝对真值”它本质上是模型预测结果。拿到数据后我建议先做几项快速检查第一用分区统计对比几个典型行政区的平均值和你对当地土壤状况的经验认识比如东北黑土区有机碳含量应该明显高于西北荒漠区第二检查极端值比例是否异常比如有没有负值的容重、有没有大于2.0的容重第三随机抽取若干点位和已有文献或野外实测数据做直观对比。这样虽然不能完全验证数据精度但至少能把重大错误挡在门外。6. 这套数据后续还能怎么扩展使用最后分享一个我最近在用的新思路。这套数据虽然是1km分辨率但完全可以和更高分辨率的遥感数据结合实现一定程度的尺度降尺度。比如用NDVI、DEM、土地利用数据做回归模型把1km土壤属性降尺度到250m或者100m当然这个过程要注意不确定性传递不能盲目追求高分辨率而忽略预测误差。这个方向很多课题组都在探索如果你想深入研究土壤属性空间预测这套数据其实是最好的训练和验证样本库。从我个人的使用经验来说这套数据最大的价值不是那16个栅格文件本身而是它提供了一种“用可重复的空间预测方式做土壤属性制图”的范式。你不需要再到处翻几十年前的土壤普查报告不需要对着缺漏百出的属性表发愁只需要把数据下载好、处理好、嵌入到自己的工作流里就能腾出精力去解决真正的科学问题。希望这篇拆解能帮你少走点弯路让数据在你手里发挥出应有的价值。
