简介这是一份面向水文模型研究、农业规划、环境评估与城乡规划等场景的土壤专题数据包旨在为科研与实践者提供全国尺度的土壤本底资料。内容覆盖土壤类型图、质地比例、容重、渗透率、含水量、pH、有机质及氮磷钾养分等关键参数可用于径流模拟、地下水位分析、农田施肥推荐、土壤退化监测与基础设施承载力评估其中红壤、黄壤、黑土、棕壤等类型分布信息可辅助识别区域差异地理坐标与海拔字段便于和DEM、土地利用数据叠加分析。压缩包共37个文件以adf栅格数据为主辅以dat属性数据、nit索引文件、xml元数据与doc说明文档整体体积约9.25MB目录结构保持ArcGIS原生组织方式可减少格式转换成本。数据包内还整理有土壤类型代码表、侵蚀数据和示意插图便于快速核对字段含义与空间分布特征。目前已有7540人学习下载适合需要建立区域水文或环境模型、但缺乏标准化空间土壤参数的研究者、规划师和工程师使用。1. 中国土壤数据集.rar一个压缩包里的全国土壤家底做环境模型、耕地质量评价或水文模拟的人几乎都遇到过同一个场景好不容易拿到一份中国土壤数据集的.rar压缩包满怀期待地解压结果蹦出来几十个TXT、几张TIF、一堆乱码文件名元数据文档还不全。这时候才发现真正的门槛不是数据本身而是怎么把它变成能进模型、能制图、能统计的标准化数据。这类压缩包装的是全国范围的土壤类型分布、理化属性、剖面观测等整理结果价值在于帮你省掉逐篇翻文献、逐县抄统计资料的时间直接拿到一份可入库、可制图、可做模型输入的数据底稿。适合环境建模、农业区划、土壤修复评估的从业者也适合刚入门遥感与GIS的研究生。要提醒的是拿到手第一件事不是解压而是先搞清楚这里面装的坐标系、土壤分类体系和缺测值标记这三样才是真正的坑。2. 解压之前先看懂数据组织从 .rar 到可用的土壤属性表2.1 先看文件清单.rar 里通常装了什么中国土壤数据集这类压缩包虽然来源不同、整理年份不同但内部结构大概率逃不出四类内容土壤类型分布图、土壤理化属性表、剖面观测记录、元数据说明文档。土壤类型图常见的是矢量面文件shp或栅格文件tif理化属性表多为TXT或CSV字段名可能混着中文、拼音、英文缩写剖面记录则是一层一行包含深度、质地、pH、有机质等。先别急着全部解压我一般会先列出压缩包内容确认文件数量和目录层级避免后面出现磁盘空间不足或解压到一半才发现套了好几层压缩包的尴尬。# 先测试压缩包是否完整不完整时不要继续 rar t 中国土壤数据集.rar # 查看压缩包里的文件清单确认嵌套结构和文件大小 rar l 中国土壤数据集.rar # 如果系统没有 rar用 unrar 效果相同 unrar l 中国土壤数据集.rar # 确认没问题后解压x 表示保留压缩包内的目录结构 mkdir -p /data/soil unrar x -o 中国土壤数据集.rar /data/soil/参数说明t是test只校验不输出文件能提前发现压缩包损坏l是list列出全部文件路径这一步能让你看到是否存在“套娃压缩包”x表示按原路径解压-o表示覆盖已有文件避免重复解压时中断。如果你在Linux下解压后发现文件名是一堆乱码不要怪数据集多半是压缩包在Windows下用GBK编码打包而你当前的locale是UTF-8。我一般改用unar代替unrar它能自动检测编码并转成UTF-8文件名比手动convmv转码省事得多。解压完成后建议先建一个data/目录把栅格、矢量、表格、文档分开放并在同级写一个README.txt记录你后续做过的所有重投影、字段重命名和缺测值处理动作。这个习惯会在你三个月后回看数据时救你一命因为这类数据集几乎都没有完善的版本记录你自己不记就没人帮你记。2.2 统一坐标与字段命名预处理的第一步解压出来后很快会碰到第二个问题字段命名极其不统一。同一份数据里“有机质”可能写成“有机质”“OM”“omm”pH可能写成“PH”“ph”全氮可能是“TN”“N_total”“全N”。做数据整合时这些差异会让pandas按列名合并时直接报错或悄悄合并出空值。我的做法是先把所有表格读进来统一做一遍字段重命名转成一套自己项目里的标准命名再往下走。import pandas as pd # 常见编码Windows 下导出的表格多为 GBK/GB2312先试 gbk df pd.read_csv(soil_property.txt, encodinggbk, sep\t) # 字段名统一把中文、拼音、英文缩写映射到项目标准名 renames { 经度: lon, lon: lon, x: lon, 纬度: lat, lat: lat, y: lat, 有机质: OM, om: OM, omm: OM, 有机质_gkg: OM, pH: pH, PH: pH, ph: pH, 全氮: TN, tn: TN, N_total: TN, 速效磷: AP, ap: AP, olsen_p: AP, 速效钾: AK, ak: AK, 阳离子交换量: CEC, cec: CEC, CEC_cmolkg: CEC, } df df.rename(columnsrenames) # 看一眼有没有漏掉的字段名 print(df.columns.tolist())这里有个容易忽略的点字段名映射不能只做一层同一个含义在不同子表里可能用完全不同的缩写比如速效磷可能叫AP也可能叫Olsen_P黏粒可能叫clay也可能叫粘粒。我一般会把所有表格的列名先集合起来打印一遍肉眼扫一批再补映射字典完全靠自动匹配容易漏。字段统一之后马上要查一遍单位这是比字段名更隐蔽的坑。有机质的单位可能是g/kg、%、mg/kgpH本身无量纲还好速效磷和速效钾经常在mg/kg与ug/g之间换算CEC则在cmol/kg和mmol/kg之间差着数量级。我习惯在字段名后面直接加单位后缀比如OM_gkg、AP_mgkg这样后续融合不同来源的数据时单位差异在列名里就能发现而不是等到统计结果出现离谱值再回头查。2.3 属性字段与土壤分类体系的对应中国土壤数据集里最让人头疼的不是理化属性而是土壤类型字段。国内常用的分类体系至少有三套发生学分类黑钙土、黑土、红壤、黄壤等、中国土壤系统分类均腐土、铁铝土等、以及国际上常见的USDA土纲。同一个采样点在这三套体系里的名称完全不一样而数据集往往不告诉你某一张表用的是哪套体系。我整理过一份最小对照表覆盖大部分常见类型供做字段映射时参考发生学分类中国土壤系统分类常见代码/缩写大致分布黑土均腐土Ht / 黑土东北地区黑钙土钙积均腐土Ch / 黑钙土内蒙古东部及东北西部栗钙土干旱土K / 栗钙土内蒙古高原红壤铁铝土R / 红壤华南及西南黄壤铁铝土Y / 黄壤贵州、四川等地水稻土水耕土P / 水稻土长江中下游及华南潮土潮土F / 潮土黄淮海平原注意上表只是帮助理解不同数据集里的代码含义可能相反。有的数据集用汉语拼音首字母比如黑土写HT有的用英文译名缩写比如Phaeozems还有的直接存中文。处理时我一般单独建一个soil_type_dict字典把一张表里出现的所有类型值先value_counts()列出来再逐个映射到目标体系绝不猜。你花半小时把类型字段梳理清楚后面做面积统计、分区制图时才不会出现“红壤”和“铁铝土”被当成两个类别的低级错误。3. 把点数据变成空间数据栅格化与插值的完整落地方案3.1 从 Excel/CSV 到矢量点最稳的入坑命令土壤数据集里的理化属性很多是以“站点经纬度”的形式存放的一行一个样点。这样的表格本身不是空间数据你得先把经纬度列变成几何对象再导出成矢量文件后续才能做空间连接、栅格值提取、插值。这一步的逻辑很简单但有一个前提必须先确认经纬度字段是十进制度如116.25还是度分秒116°1500还是投影坐标如UTM的米制坐标。数据集里存度分秒的虽然少但遇到过换算比想象的更麻烦。import pandas as pd from shapely.geometry import Point import geopandas as gpd # 确保经纬度是浮点数不是字符串 df[lon] pd.to_numeric(df[lon], errorscoerce) df[lat] pd.to_numeric(df[lat], errorscoerce) # 删除经纬度缺失的行这一步必须放在构造几何之前 df df.dropna(subset[lon, lat]) # 构造 GeoDataFrame先按 WGS84 经纬度定义 gdf gpd.GeoDataFrame( df, geometry[Point(x, y) for x, y in zip(df[lon], df[lat])], crsEPSG:4326 ) # 导出为 shapefile注意编码统一 gdf.to_file(soil_points.shp, encodingutf-8)这段代码里errorscoerce会把无法转成数字的经纬度变成NaN随后的dropna把它们清理掉。有些样点经纬度填的是文本格式如“116.25E”“45.5N”pd.to_numeric会直接转成NaN所以更稳妥的是在转数值之前先做一个正则清洗去掉E/N/S/W后缀。crsEPSG:4326表示当前几何是WGS84经纬度坐标系这是最通用的声明方式如果数据说明里写的是CGCS2000建议用EPSG:4490CGCS2000地理坐标系别混。导出shp时设encodingutf-8否则属性表里的中文在ArcGIS里会乱码。3.2 栅格化与属性重分类制图分级的关键参数样点矢量做好了下一步通常是把某个属性比如有机质、pH栅格化成一张全国或区域尺度的连续分布图。这里有一个常见的误解直接用样点做栅格化得到的只是“落在该像元内的样点原值”其余像元全是空值根本没有空间连续性。要得到真正的连续土壤属性图得走“插值”这一步。我通常的做法是先用IDW或克里金把样点插值成连续表面再做栅格化或直接输出GeoTIFF。# 用 gdal_rasterize 可以把点属性写入指定像元适合做点值转栅格 gdal_rasterize -a OM_gkg -tr 1000 1000 -a_nodata -9999 \ -ot Float32 soil_points.shp soil_om_1km.tif参数说明-a OM_gkg指定要用哪个属性字段写入栅格值-tr 1000 1000设置像元大小是1000米×1000米单位与矢量数据的投影坐标系一致所以要确保此时矢量已经转成了投影坐标比如gdf.to_crs(EPSG:32650)而不是经纬度否则这个1000会被当成“1000度”输出结果完全不可用-a_nodata -9999设置无数据值为-9999这一项必须设否则空白区域默认值可能是0后续统计时会严重拉低均值-ot Float32指栅格数据类型。如果你想要的是连续分布栅格正确流程应该是先插值再把插值结果导出成tif而不是直接gdal_rasterize。插值这一步我一般用pykrige做普通克里金或者偷懒用scipy.interpolate.Rbf做径向基函数插值。重点是设置合适的分辨率比如全国尺度用1km区域尺度用100500m。这里还要注意某些数据集已经帮你做好了全国1km或250m的栅格产品那就不需要插值直接检查坐标系和nodata值即可省掉一大段工序。3.3 缺测值与异常值的处理别让 -9999 污染你的统计土壤理化数据里缺测标记是五花八门的-9999、-999、-8888、9999、-1、空字符串甚至0在某些字段里本身就是缺测。如果不先做缺测值清洗直接df.describe()你会在有机质字段里看到负数的均值或者在pH里看到最大值为99.99这种离谱数。清洗的规则我一般这么定先把所有常见缺测标记统一替换成NaN再按字段做合理性区间过滤。import numpy as np import pandas as pd # 常见缺测标记统一次替换 missing_marks [-9999, -999, -8888, 9999, -1] for col in [OM_gkg, pH, TN, CEC, AP_mgkg, AK_mgkg]: if col in df.columns: df[col] df[col].replace(missing_marks, np.nan) # 按土壤学常识做范围过滤 # pH 一般在 3.5~10.5 之间超出即为异常 df df[(df[pH] 3.5) (df[pH] 10.5)] # 有机质含量表层土壤 g/kg 一般在 0~600 之间超过 600 需要复查 df df[(df[OM_gkg] 0) (df[OM_gkg] 600)] # 速效磷和速效钾允许高值但为负值或异常大时应置为 NaN for col in [AP_mgkg, AK_mgkg]: df.loc[df[col] 0, col] np.nan这段清洗可以帮你过滤掉大部分明显错误但有一个点必须强调范围过滤只能抓物理不可能的值抓不了“合理但错误”的值。比如某样点的pH实际是5.2表格里误录成5.8这个5.8落在合理区间内任何自动清洗都发现不了。因此清洗之后我建议把每个字段按土壤类型分组看一眼均值如果某个类型的有机质均值比同类型其他区域高出几倍那大概率是单位或录入错误。清洗不是一次性动作而是要保留清洗前后的对比统计方便追溯。4. 中国土壤数据集落地应用的五个高频坑与排查清单4.1 坐标系统一问题WGS84 与 CGCS2000 混用导致分布偏移现象把数据里的样点叠加到当地行政区边界或底图上发现采样点整体偏移了几十到几百米甚至有的点落到河流或道路中间。在ArcGIS或QGIS里单独看每个图层都没问题一叠加就错位。原因数据集内不同表格或矢量文件的坐标系不一致。有的文件是WGS84经纬度有的是CGCS2000经纬度有的甚至已经是投影坐标如UTM或高斯-克吕格投影但元数据里没写清楚。CGCS2000与WGS84在大部分区域差异只有几十厘米到几米真正的偏移往往来自“把投影坐标当经纬度”或“把经纬度当投影坐标”的使用错误。解决先用ogrinfo或gdalinfo挨个查看各文件的坐标系定义凡是没写坐标系的一律视为可疑。然后以数据里占多数的坐标系为基准用gdalwarp栅格或ogr2ogr -t_srs矢量统一转换。我习惯把原始文件保留一份转换后的统一存为EPSG:4490或EPSG:4326并在文件名里加后缀如_wgs84方便区分。4.2 土壤分类体系混用发生学分类与系统分类的字段坑现象做面积统计时发现同一个采样点在两张表里被标成不同土壤类型比如一张表写“红壤”另一张表写“铁铝土”或者按类型上色制图时图例里出现大量“未分类”颜色块。原因这是数据集在整理时没有统一分类体系导致的。发生学分类是老的体系系统性分类是新的体系两张子表各用各的字段值是中文名还是代码也不一致直接合并就会产生“同一事物两种名称”的冲突。解决这一步没有捷径只能自己建对照字典。先把所有表格的类型字段value_counts()列出来把出现的每一个值都记录进字典再映射到你项目要用的目标体系。对于无法判断的类型宁可就标“未知并发层”也不要猜。我吃过一次亏某次直接把英文缩写当成系统分类代码处理后来发现那是发生学分类的拼音缩写整个分析推倒重来。4.3 栅格分辨率与重采样方法选择最邻近还是双线性现象把1km分辨率的全国土壤属性栅格重采样成250m或100m后得到的图出现明显的锯齿状边界反过来把250m栅格升尺度到1km又出现属性值被平均后失去极值的情况。原因重采样方法选错了。类别型数据比如土壤类型代码和连续型数据比如有机质含量对重采样的要求完全不同。双线性或三次卷积会把类别代码插出小数比如土壤类型1和2之间插出1.5这在语义上是无意义的而最邻近法用于连续数据会保留原始像元值但分辨率变化大的时候会产生方块效应。解决土壤类型、地质类型这类名义变量统一用最邻近nearest重采样或者配合众数重采样有机质、pH、黏粒含量这类连续变量用双线性或立方卷积。务必在重采样前先用gdalinfo确认nodata值被正确设置否则空白区域会被当成0参与插值把边界像元彻底污染。这一条如果忽略最终制图或模型输入会带上明显的条带或方块伪影。4.4 大文件解压与内存溢出.rar 的压缩方式与读取策略现象解压一个几十GB的压缩包解压过程中磁盘空间耗尽用pd.read_csv(big_table.txt)读一张几个GB的表格直接把内存占满进程被杀。原因这类数据集里的文本表格几乎都是高度可压缩的压缩包里看着不大解压后体积翻好几倍。另外土壤属性表经常是“宽表”列数多、行数也多直接全量读入内存对普通办公电脑极不友好。解决解压前先rar l看压缩包内文件大小总和确保磁盘有23倍冗余。读大表格时不要一次性全量读用chunksize分块读或者干脆先用csv模块扫一遍结构再决定只读需要的列。对于要反复使用的表格我一般会转成Parquet格式读起来比CSV快很多而且自带压缩。经验是任何超过500MB的CSV都不值得在内存里反复读。# 分块读取大表格只保留需要的列 cols_needed [lon, lat, OM_gkg, pH, soil_type] chunk_iter pd.read_csv(big_table.txt, sep\t, chunksize10000, usecolscols_needed) pieces [chunk[chunk[OM_gkg].notna()] for chunk in chunk_iter] df pd.concat(pieces, ignore_indexTrue)4.5 深度与土层字段的读取错误剖面数据不是平面数据现象用剖面数据做统计时发现同一经纬度出现多行记录直接合并时样点数量翻倍绘制土壤属性图时表层和深层数据混在一起图面一片杂乱。原因剖面数据本身是一层一行记录的是同一采样点不同深度的分层属性比如0~20cm、20~40cm、40~60cm。如果没有按深度字段区分直接当成平面点位数据使用就会把同一个点的多个土层当成多个点。解决先确认数据是“剖面结构”还是“表层结构”。如果是剖面先把目标土层筛出来比如只保留0~20cm或者按深度字段做分桶统计再合并到表层的样点表里。不要试图把多层剖面数据直接聚合到平面属性表除非你的模型明确需要多层土壤参数。我处理这类数据时会单独建一个profile文件存放所有土层信息而平面属性表只保留表层样点这样两边都不会乱。5. 用一张验证图收尾把数据集精度跑出来的最小流程数据集拿到手清洗也做完了怎么判断它能不能直接进模型我的做法是做一个最小验证流程随机抽一批样点把数据集栅格上的提取值与这批样点的实测或高精度参考值放在一起计算相关系数、均方根误差画一张1:1散点图。这张图能在一分钟内暴露坐标偏移、单位错误、缺测污染三类问题比任何目视检查都有效。import numpy as np import pandas as pd from osgeo import gdal ds gdal.Open(soil_om_1km.tif) band ds.GetRasterBand(1) gt ds.GetGeoTransform() # 从验证点表读取经纬度并转为像素坐标 vals [] for x, y in zip(vdf[lon], vdf[lat]): px int((x - gt[0]) / gt[1]) py int((y - gt[3]) / gt[5]) v band.ReadAsArray(px, py, 1, 1)[0, 0] vals.append(v if v -9990 else np.nan) vdf[extracted] vals valid vdf.dropna(subset[extracted, obs]) r2 np.corrcoef(valid[extracted], valid[obs])[0, 1] ** 2 rmse np.sqrt(((valid[extracted] - valid[obs]) ** 2).mean()) print(fR2{r2:.2f}, RMSE{rmse:.2f})这里的核心是像素坐标换算px (x - gt[0]) / gt[1],py (y - gt[3]) / gt[5]分别对应东西向和南北向的像元分辨率。注意gt[5]通常是负值不要写成gt[4]否则南北方向会翻转。如果提取值大量为NaN先怀疑坐标参考不一致再怀疑nodata设置最后才怀疑验证点本身有问题。R²在0.6以上、RMSE相对于数据自身标准差明显偏小时这个数据集就可以放心进入后续建模流程。我个人的习惯是拿到任何土壤数据集都会先跑这个验证流程连“官方整理好”的栅格产品也不跳过。某次做某流域水文模拟时我直接用了别人处理好的全国土壤栅格没做验证结果模型率定参数怎么调都不稳定查了两天才发现是栅格与流域边界坐标系相差一个像元的偏移导致每个流域的面积参数整体错位。从那以后验证图成了我的固定第一步。希望这条流程也能帮你在数据集落地的过程中少走一段弯路。本文还有配套的精品资源点击获取
