简介SOTER数据库是一套面向国土规划、农业区划与生态建模的土壤和地形空间数据包适用于GIS从业者、土壤学及生态学研究者。内容涵盖全国地貌、岩性与土壤类型的分布图其中地貌分为平原、丘陵、山地、内陆水域等8类岩性含酸性火成岩、变质岩、沉积岩等58种组合土壤类型覆盖低活性强酸土、火山灰土、冲积土、铁铝土、潜育土等数十类可支撑多尺度环境因子叠加分析与区域对比统计。资源以rar压缩包发布约84.51MB便于直接下载与解压引用。目前已有154人学习使用。数据可直接导入主流GIS平台作为中国区域土壤-地形数字化底图用于生态区划、水土保持、农作物适宜性评价、土壤侵蚀模拟等研究场景也可配合空间分析工具完成专题制图与成果输出。1. 什么是SOTER数据库为什么传统土壤图满足不了地形与土壤的一体化分析SOTERSoil and Terrain Database中文通常叫“土壤和地形数据库”全国尺度的这套数据不是简单把土壤图叠在DEM上而是把地形、母质、土壤剖面按统一的空间单元组织起来形成一套能直接进模型、能按单元统计、能跨区域对比的数据库。做过水土流失评估或耕地质量评价的人应该都有体会手里有一张1∶100万土壤图却查不到某块地的坡度、坡向、有效土层厚度更查不到这个多边形对应的母质和剖面分层只能临时去翻纸质剖面记录再手工和DEM叠加。SOTER就是冲着这个痛点来的。本文适合正在做土壤普查成果转化、区域资源环境承载力评价或者准备把土壤数据接进SWAT、RUSLE等模型的人读完可以知道这套数据结构长什么样、怎么从零搭起来、用的时候哪些参数最容易被坑。2. SOTER数据库的底层逻辑单元划分与属性表设计2.1 为什么不用直接叠加SOTER单元的“地形-岩石-土壤”三级体系常见的土图与DEM叠加是把土壤图多边形和坡度栅格转成的坡级多边形做空间求交得到的每个图斑只是“某类土壤某个坡段”。这种做法有两个问题一是坡度分级阈值不同结果图斑数量会差出两三倍二是图斑内部的地貌部位、坡形、母质可能完全不一致比如同一个“平缓坡”图斑里既能出现在山顶平台也能出现在河谷阶地后续分析很难解释。SOTER的划分思路是先把地形本身做“内部分异”。一个完整的地形单元要按“地貌类型—坡位—坡度等级—坡形—起伏强度”逐层拆分。我常用的做法是先用地貌位置分类提取山脊、坡肩、背坡、坡脚、谷底再用坡度等级细分得到一组地形单元然后把岩性图或母质图重分类后与地形单元叠加形成“地形—母质”复合单元最后再叠土壤图得到真正的SOTER单元。这样每个SOTER单元内部的地形与母质条件相对均一土壤组合也稳定做出来的统计结果不会因为换一种坡度阈值就面目全非。这里的关键参数是坡度分级和地貌部位的判定。国内许多建库项目参照地貌制图规范把坡度分成0°2°、2°5°、5°15°、15°25°、25°35°、35°六档地貌部位用基于DEM的正负地形和相对位置判别比如用曲率区分凸形和凹形坡用相对坡位指数区分脊、坡、谷。大家抄作业的时候建议把这套阈值写进元数据因为后面任何分析都跟它绑死改一次就得重来一遍。2.2 属性表结构从SOTER单元到地形组成、土体剖面SOTER数据库的空间对象是SOTER单元多边形属性则分层存在几张表里。最常见的是三张核心表soter_unit主表、terrain_component地形组成表、pedon_profile土体剖面表。主表的每个多边形有一条记录存单元ID、区域代码、面积、地貌类型、母质类型和几何字段。terrain_component表则描述这个单元内部的微地形因为一个SOTER单元里可能包含几个坡段组分而pedon_profile表记录的是单元内的代表性土壤剖面及土层分层数据。用SQL建表的骨架大致如下字段名是各家常见的写法实际项目里可以根据数据库规范改名但主键与外键关系建议保持一致CREATE TABLE soter_unit ( soter_id VARCHAR(20) PRIMARY KEY, region_code VARCHAR(12), unit_area_ha NUMERIC(10,2), landform_type VARCHAR(20), parent_material VARCHAR(30), geom GEOMETRY(POLYGON, 4610) ); CREATE TABLE terrain_component ( soter_id VARCHAR(20) REFERENCES soter_unit(soter_id), component_id INTEGER, slope_class VARCHAR(10), slope_length_m NUMERIC(8,1), relief_type VARCHAR(20), surface_roughness VARCHAR(10) ); CREATE TABLE pedon_profile ( pedon_id VARCHAR(20) PRIMARY KEY, soter_id VARCHAR(20) REFERENCES soter_unit(soter_id), horizon_no INTEGER, top_depth_cm NUMERIC(5,1), bottom_depth_cm NUMERIC(5,1), texture_class VARCHAR(20), organic_carbon_pct NUMERIC(5,2), ph_h2o NUMERIC(4,2), cec_cmol_kg NUMERIC(6,2), gravel_pct NUMERIC(5,1) );SQL里的soter_id是全局唯一标识建议用“CN行政代码顺序号”例如CN350122001terrain_component里的component_id表示该单元内第几个地形组分pedon_profile里的horizon_no表示第几层土壤。这样设计的好处是一个单元可以对应多个组分、多个剖面而剖面数据不必每个图层都复制一遍逻辑清楚数据冗余也少。2.3 主键与外键的设计细节如何保证连得上、查得全做全国尺度的SOTER数据库最怕的就是属性表之间对不上。主表里明明有3000个单元剖面表却只有2500个剩下500个是“没有剖面”还是“ID不匹配”这就要靠外键约束来兜底。上面SQL里已经写了REFERENCES建库时如果数据库引擎支持一定要打开外键约束像PostgreSQL配合PostGIS就很合适。如果用的是空间文件格式比如Shapefile加dbf表外键只能靠人工维护那就更得在录入阶段反复检查。我一般会在入库后跑一条验证SQL把孤儿记录查出来SELECT p.pedon_id, p.soter_id FROM pedon_profile p LEFT JOIN soter_unit s ON p.soter_id s.soter_id WHERE s.soter_id IS NULL;另外单元与剖面的关系是一对多一个SOTER单元里如果只有一条剖面分析时可以把这条剖面直接代表整个单元如果有两条以上就要考虑按面积权重或按地貌组分进行聚合。建议在属性设计时给剖面表增加一个代表级别字段标记该剖面是“典型”还是“辅助”避免后续处理时把所有剖面都当成独立样本。3. 从零构建一份全国尺度的SOTER数据库数据准备与实操流程3.1 数据源、投影与预处理构建全国尺度的SOTER数据库至少需要四类数据DEM、土壤图、岩性图/母质图、地质地貌辅助资料。DEM一般用30m或90m分辨率全国尺度我建议先统一到90m计算量小坡度信息也够用如果做省市级换30m更好。土壤图优先用已数字化的1∶100万或1∶25万土壤类型图没有矢量版就扫描底图后在GIS里地理配准和矢量化这一步最耗时但也是后面不翻车的前提。岩性图可以用1∶50万或1∶100万地质图的岩性属性重分类。所有数据必须在同一投影下操作。全国尺度我常用Albers等积圆锥投影中央经线105°E双标准纬线25°N和47°N坐标系用CGCS2000。等积投影保证面积统计不变形坡度、坡长这类几何量也能保持相对准确。用GDAL转换DEM的示例命令如下gdalwarp -t_srs projaea lat_125 lat_247 lat_00 lon_0105 datumCGCS2000 \ dem_raw.tif dem_aea.tif这里的参数中lat_1和lat_2是双标准纬线lon_0是中央经线如果不指定重采样方法gdalwarp默认用最近邻坡度计算建议加上-r bilinear或-r cubic避免DEM出现锯齿。转换完成后还要检查土壤图与DEM的边界是否大致吻合差几百米是常事需要用控制点做配准这一步属于纯手工活不能偷懒。3.2 地形单元提取与叠加一个Python示例地形单元提取的完整流程是先由DEM计算坡度再按分级表把坡度栅格重分类同时用地貌分类算法提取坡位两者加岩性图一起叠加。这里给一个直接用rasterio和geopandas处理的Python脚本示例适合已经准备好坡度栅格、岩性重分类栅格、土壤图矢量的情况import rasterio import numpy as np import geopandas as gpd from rasterio.features import polygonize from shapely.ops import unary_union # 读取坡度栅格 with rasterio.open(slope_aea.tif) as src: slope src.read(1).astype(float32) transform src.transform crs src.crs # 坡度分级: 0,2,5,15,25,35,35 classes np.zeros_like(slope, dtypenp.uint8) classes[(slope 0) (slope 2)] 1 classes[(slope 2) (slope 5)] 2 classes[(slope 5) (slope 15)] 3 classes[(slope 15) (slope 25)] 4 classes[(slope 25) (slope 35)] 5 classes[slope 35] 6 # 矢量化坡度分级结果 shapes polygonize(classes, transformtransform) slope_poly gpd.GeoDataFrame.from_features(shapes, crscrs) slope_poly[slope_class] slope_poly[value].astype(int) # 与岩性图叠加 litho gpd.read_file(litho_aea.shp) soter_intermediate gpd.overlay(slope_poly, litho, howintersection)这段代码的逻辑是先读入坡度栅格用布尔索引做多级分类然后polygonize把分类栅格转成矢量多边形最后与岩性图求交得到“地形—母质”复合多边形。参数上slope数组中的NaN值比如DEM有空洞会被归到classes0建议提前用邻域插值把空洞补掉否则会产生大片空值区。collect之后不要直接当最终SOTER单元还要继续叠土壤图叠加完再设置最小图斑面积比如1∶100万尺度下小于100公顷的碎斑要合并。这一步有个常见的玄学问题栅格转矢量后每条边界都是锯齿多边形数量巨大。我通常会在矢量化后对每个多边形做简化用shapely的simplify方法容差设为像元大小的一半比如90m分辨率容差45mslope_poly[geometry] slope_poly.geometry.simplify(45, preserve_topologyTrue)preserve_topologyTrue能防止简化后出现自相交。如果经费允许也可以用ArcGIS的“栅格转面消除”工具链操作步骤更直观但批量处理还是Python脚本更可靠。3.3 属性表录入与批量匹配SOTER单元多边形生成后接下来要给每个单元填属性。地形属性可以从DEM批量提取比如每个单元的平均坡度、坡向、粗糙度、起伏度土壤属性则来自剖面数据库。常见做法是先用Excel整理野外或历史剖面记录再写脚本连接。这里给出用pandas连接两个表格的代码import pandas as pd terrain pd.read_csv(terrain_attr.csv, dtype{soter_id: str}) pedon pd.read_excel(pedon_records.xlsx, sheet_nameprofile, dtype{soter_id: str, horizon_no: int}) pedon pedon.sort_values([soter_id, horizon_no]) # 每个剖面取最上层数据合并到单元表 profile_top pedon.drop_duplicates(pedon_id, keepfirst) merged terrain.merge(profile_top, onsoter_id, howleft, validateone_to_many)关键点在于dtype{soter_id: str}这一步能避免“00123”和“123”匹配不上的问题。dataframe的merge里validateone_to_many表示左侧一个单元对应右侧多条记录如果右侧出现重复soter_id且不是期望的一对多会直接报错这是好事能在入库前暴露ID错误。合并完成后还要检查合并空值比例如果超过20%说明剖面数据库覆盖不够要么补采要么用后文提到的土壤-地形规则推算。属性表最终要写回空间图层。用geopandas的话直接执行soter_final soter_gdf.merge(merged, onsoter_id)再导出为GeoPackage或Shapefile。这里提醒一句不要用Excel的vlookup去和Shapefile的dbf属性表做关联因为dbf的字符编码和长度限制很容易让中文变成乱码或截断用Python或直接SQL完成更稳。4. 用SOTER数据库算一次水土流失风险字段筛选与参数换算4.1 从TERRAIN表提取LS因子拿到SOTER数据库后一个典型的应用是计算区域水土流失风险通用水土流失方程RUSLE里最难算的是LS因子坡度坡长因子。传统算法是拿DEM逐像元算坡度和坡长但这样算出来的因子忽略土壤类型和地形单元的边界经常出现一个单元内LS值剧烈抖动。用SOTER表的优势在于每个单元已经有了平均坡度和坡长属性直接用这些值代表整个单元再参与模型计算结果符合地貌单元内部的均质性假设。我一般会写这样一个函数import numpy as np def compute_ls_from_soter(slope_pct, slope_length_m): # slope_pct: 百分比坡度, 来自terrain_component表 # slope_length_m: 坡长, 来自terrain_component表 theta np.arctan(slope_pct / 100.0) # 坡长指数 m 随坡度变化 if slope_pct 5: m 0.5 elif slope_pct 10: m 0.4 else: m 0.3 L (slope_length_m / 22.13) ** m if slope_pct 9: S 10.8 * np.sin(theta) 0.03 else: S 16.8 * np.sin(theta) - 0.5 return L * S这个函数里22.13m是RUSLE标准坡长m取0.5/0.4/0.3是经验值分别对应缓坡、中坡、陡坡。S公式在不同坡度区间有两种形式9%是切换阈值。用到SOTER数据时建议从terrain_component表里取slope_class对应的中间值或者用单元平均坡度而不是直接输入某一像元值。如果数据库里只有坡度等级没有连续坡度可以用等级的中值代替比如5°15°取10°。4.2 从PEDON表计算土壤可蚀性K因子K因子表示土壤对侵蚀的敏感性SOTER的pedon表里有质地、有机碳、pH、CEC、砾石含量这些正好能算K。精确的K值需要用土壤可蚀性诺模图但工程上常用简化公式。从SOTER表出发我更推荐先用质地分类给基础K值再用有机碳修正质地类型基础K值(t·ha·h/(MJ·mm·ha))砂壤土0.15壤土0.25粉砂壤土0.30粉砂黏壤土0.32黏壤土0.28黏土0.30比如用Python来算k_map { sandy_loam: 0.15, loam: 0.25, silt_loam: 0.30, silty_clay_loam: 0.32, clay_loam: 0.28, clay: 0.30 } def calc_k_from_pedon(texture_class, organic_carbon_pct): base_k k_map.get(texture_class, 0.25) # 有机质增加会降低K, 简化为每1%有机碳降低3% k base_k * (1.0 - 0.03 * organic_carbon_pct) return round(max(k, 0.05), 4)这个简化公式里有机碳每增加1个百分点K降低3%是经验值不是严格物理公式。真正要发论文或做工程报告最好用本地径流小区的实测数据校准否则只能当相对风险排序用。计算时注意pedon表里可能有多层土壤应该用表土层horizon_no1的质地和有机碳而不是整段土壤的平均值。算完每个单元的K和LS后再把它们连接到SOTER单元多边形栅格化或直接用多边形属性参与运算。R因子降雨侵蚀力可以从气象站降雨数据插值C因子和P因子取自土地利用最后把五个因子相乘得到年土壤侵蚀模数erosion r_factor * k_factor * ls_factor * c_factor * p_factor这里的r_factor、ls_factor都可以按单元取值c和p建议按土地利用类型赋常量。最终结果是一个每平方千米吨数t/km²·a的空间分布。4.3 出图与验证分级阈值怎么定算出侵蚀模数后不能直接拿自然断点分级因为不同地貌区的侵蚀强度差别很大。国内常用的土壤侵蚀分类标准把年均侵蚀模数分为微度200轻度2002500中度25005000强度50008000极强度800015000剧烈15000单位t/km²·a。但这个标准主要面向黄壤区华东华南的红色风化壳地区阈值要适当上调否则整个山区几乎全是剧烈侵蚀失去区分度。出图时我一般会用分省或分流域的统计表配合SOTER单元做对比。验证手段有两个一是与当地水文站的多年来沙量数据对比用流域出口输沙量除以流域面积得到一个平均侵蚀模数看模型结果的数量级是否一致二是和已有高分辨率遥感解译的侵蚀图斑做空间叠置计算Kappa系数。只要SOTER单元的LS和K因子没有明显的空间突变通常结果都比较稳。5. SOTER建库最容易翻车的4个现场现象、原因、解决5.1 叠出来的图斑又碎又破拓扑错误遍地第一次做SOTER叠加时最容易遇到的情况是土壤图有8000个多边形坡度分级栅格转矢量后有2万个碎斑叠完之后直接变成十几万个图斑打开属性表一卡一卡的。原因很简单栅格边界是锯齿状和土壤图的平滑边界相交会产生大量狭长面土壤图本身的图斑边界也可能因为配准误差错位几毫米叠加后就切出无数微多边形。解决方法是先做栅格化聚合再把坡度多边形简化最后叠加后设置最小图斑面积。我用的是两个环节一是在矢量坡度图上用eliminate工具把小于最小制图面积的多边形合并到相邻最大边界二是叠加完成后再来一轮“按面积排序逐个合并”。最小制图面积要与比例尺匹配1∶100万取100公顷1∶25万取25公顷。另外叠加前最好把所有图层的节点坐标Round到1米精度减少无效碎边。5.2 剖面数据连接后大片空白属性连接是另一个高发“翻车”现场。现象是所有剖面记录在Excel里看着都对但merge之后发现60%的soter_id匹配不上。原因往往是ID的格式不一致SOTER单元ID可能是“CN350102001”Excel里被显示成科学计数法数字后面多了个“.0”或者有人不小心在ID前输入了全角空格。还有一个坑是土壤图里单元ID是数值型剖面表里是文本型两边看起来一样实际类型不同。解决方法是统一清洗文本。pandas里这样处理pd.read_excel(..., dtype{soter_id: str})然后对ID列执行strip()去掉首尾空格再用astype(str)强制转字符串。如果发现ID里有全角字符先str.replace( , )。入库前用前文的外键检查SQL跑一遍把孤儿记录全部打印到日志文件逐条查原因再补。这个服务我做了无数遍说白了都是格式细节不是复杂算法。5.3 坡度分级改了两次整个单元统计对不上建库中期最容易产生自我怀疑的问题是调整了一次坡度分级阈值后SOTER单元数量变化了百分之三四十前面算的所有面积统计和平均属性全要重做。原因很直接SOTER单元的地形分类本质上是人为离散化阈值变化会直接导致多边形边界迁移尤其是5°15°这一档坡耕地和草地在这里极容易跨档。解决的原则是“第一次建库就定死阈值并写成文档”。把分级表、DEM分辨率、投影、最小图斑面积全部记录在元数据表里后续分析一律读取这个配置文件不要手动改。即使要调整也要另建版本新旧版本同时保留。另外地形提取的算法也要固定在同一软件版本和参数上同一个DEM在ArcGIS和GDAL里算出的坡度会有小差异绝对不能混着用。5.4 克里金插值“越插越假”有些项目剖面点很少为了把SOTER每个单元都填上土壤属性会用克里金对有机质、pH做空间插值。结果常出现pH负值、有机碳高达20%的离谱数值。原因是剖面采样点多在方便到达的缓坡和谷地山顶和陡坡几乎没有样点数据不满足克里金要求的平稳假设插值结果自然失真。解决这类问题我一般倾向用“土壤-地形关系规则”而不是空间插值。规则可以是“坡度25°的石质土单元层厚最浅有机碳低山间盆地水田土有机碳高”这些规则由土壤学家经验或已有文献提炼再按SOTER单元的母质和地形组合赋值。插值只用来生成连续变量且必须按物理范围截断比如pH限制在4.58.5有机碳限制在010%。截断代码很简单interp[pH] interp[pH].clip(4.5, 8.5) interp[OC] interp[OC].clip(0, 10)但要把这条写在数据处理文档里否则审稿人或验收专家会认为数据是“造的”实际上这是对异常插值的强制约束属于建库常识。6. 进阶给SOTER数据库加一道自动质检建库完成后如果你不想在后续分析中被各种隐藏问题反咬一口建议写一个质检脚本每次改库都跑一遍。质检表至少包括三个规则每个SOTER单元必须在terrain_component表里至少有一条记录每条pedon记录必须能关联到存在的soter_id土层深度必须自上而下递增。这些年我被ID空格坑过一次被图层投影坑过一次都是靠自动检查补回来的。import pandas as pd def run_quality_check(units, terrain, pedon): errors [] # 规则1: 单元必须至少一个地形组分 unit_have_terrain terrain[soter_id].isin(units[soter_id]) missing_terrain units.loc[~units[soter_id].isin(terrain[soter_id])] if len(missing_terrain) 0: errors.append(f缺少地形组分的单元: {missing_terrain[soter_id].tolist()}) # 规则2: 剖面必须指向存在的单元 orphan pedon.loc[~pedon[soter_id].isin(units[soter_id])] if len(orphan) 0: errors.append(f悬空剖面soter_id: {orphan[soter_id].unique().tolist()}) # 规则3: 同一剖面的土层深度必须递增 pedon[top_depth_cm] pd.to_numeric(pedon[top_depth_cm], errorscoerce) pedon[bottom_depth_cm] pd.to_numeric(pedon[bottom_depth_cm], errorscoerce) bad_depth pedon[pedon[top_depth_cm] pedon[bottom_depth_cm]] if len(bad_depth) 0: errors.append(f土层深度不递增的pedon_id: {bad_depth[pedon_id].tolist()}) return errors脚本的核心是“宁可把所有问题一次性暴露也不要让它们在下游分析时变成黑匣子”。我们建全国的库时光检查ID空格就抓出过几百条错记录后来又加了面积非负、坡度等级与坡位组合的逻辑检查误差再也没在下游应用里出现过。现在每次导出数据前我都会跑一遍这个质检脚本把错误日志存成文件然后连同数据一起交付。省下的返工时间远超过写脚本的时间希望这个习惯也能帮到你。本文还有配套的精品资源点击获取
![[Linux操作系统] 添加、修改与删除用户和用户组](http://pic.xiahunao.cn/yaotu/[Linux操作系统] 添加、修改与删除用户和用户组)
