简介全国1:100万土壤类型图数据包面向土壤、地理、生态及环境相关专业的科研人员与GIS使用者采用土壤发生分类系统覆盖全国各类土壤及其主要属性特征可用于区域土壤制图、空间分析、教学演示与科研底图。压缩包共6个文件约52.52MB核心为shp/shx/dbf/prj/xml等标准矢量数据文件配合doc数据说明文档与土壤分类代码表便于理解图层字段、分类体系和数据来源。目前已有381人学习/下载。整套资料不仅提供可直接加载的全国土壤空间数据还附带数据来源及分类代码说明帮助使用者准确调用图例、梳理属性字段降低二次整理门槛。对于需要标准化土壤类型底图或进行区域对比研究的读者这是一份实用的基础数据参考。1. 中国土壤类型图1:100万这张图到底解决什么问题做一版全国尺度的土壤空间分析最怕的不是没数据而是拿到的第一张底图就用错比例尺。中国土壤类型图1:100万是一份经过制图综合的全国性土壤分布数据集图上1厘米相当于地面10公里图斑最小面积通常在几十公顷以上。它不像高精度地块图那样能告诉你某块田的土种却能回答“这一片区域大致是红壤还是黑土”这类宏观问题。我的经验里它最适合给陆面过程模型、生态区划、大流域水文模拟做输入如果你要做精确到村庄的施肥方案这张图就太粗了。你以为它是“不精确”的图实际上它是“尺度匹配”的图。用错尺度的代价比数据本身误差更大——很多翻车现场不是投影错就是拿1:100万去当1:5万使。这篇笔记从拿到矢量数据那一步讲起落到坐标处理、分类编码、栅格化和复核照着做能少走弯路。2. 拿到1:100万土壤类型图之后先做三件事坐标、裁剪、完整性检查2.1 先别急着画图看一眼空间参考是什么全国尺度的土壤类型图在数字库里最常见的保存格式是Shapefile但不同单位分发的版本空间参考可能完全不同。有的直接存经纬度有的做了Albers等积投影还有的历史版本用了Krasovsky椭球的Gauss-Kruger投影。你第一步不是打开ArcMap渲染颜色而是先打印crs确认它是什么。如果坐标系是经纬度而你没发现后面计算面积、栅格化时的分辨率全部都会错。全国1:100万制图和后续面积统计我一般会统一转成Albers等积圆锥投影中央经线105°E双标准纬线25°N和47°N。选择等积投影而不是等角投影是因为土壤图最终要回答“某类土壤分布多少面积”和“在哪个气候带集中”这类问题面积保真比形状保真重要得多。下面是读取和重投影的代码。import geopandas as gpd # 读取矢量土壤类型图shp/gpkg均可 soil gpd.read_file(data/soil_type_1m.shp) # 先看原始坐标系 print(原始CRS:, soil.crs) print(EPSG:, soil.crs.to_epsg() if soil.crs else None) # 若图层没有坐标系先按经纬度给它一个 if soil.crs is None: soil soil.set_crs(EPSG:4326) elif soil.crs.is_geographic: print(注意原始坐标系是经纬度直接算面积会失真) # 统一转到Albers等积投影适合全国范围 albers ( projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm no_defs ) soil soil.to_crs(albers) # 投影后看一眼图斑面积量级 soil[area_m2] soil.geometry.area print(soil[area_m2].describe())这段代码的逻辑是先判断坐标系是否存在再判断是否为经纬度最后统一投影。参数是中规中矩的全国Albers投影参数没有用某个EPSG编号是因为非标准投影经常无法被to_epsg()识别直接用Proj字符串比硬找编号稳定。转完后用geometry.area检查面积如果能正常输出数值说明投影转换成功如果报错说坐标范围不对很可能是原数据根本不是经纬度却被你设成了4326这时候应该回去查元数据。2.2 用GeoPandas做完整性检查1:100万图的历史版本常常是分省分幅数字化后拼接的拼接过程会留下各种边角问题。常见的有空几何、无效几何、完全重复的图斑、编码字段空值。这些问题不清理后面做栅格化时会出现莫名其妙的黑边和覆盖顺序错乱。不要一上来就做漂亮的地图先用代码把数据“体检”一遍。# 查看字段名确认土壤编码字段 print(soil.columns.tolist()) print(soil.head(3)) # 关键字段在不同版本里叫法不同常见有 TYCODE、SOIL_CODE、TYPE key_field TYCODE # 改成实际字段 # 几何完整性 print(要素总数:, len(soil)) print(空几何:, soil.geometry.is_empty.sum()) print(无效几何:, (~soil.geometry.is_valid).sum()) # 编码字段的空值 if key_field in soil.columns: print(编码字段空值数:, soil[key_field].isna().sum()) else: print(当前没有字段, key_field, 需要人工核对) # 完全重复的几何会导致栅格化时值互相覆盖 print(完全重复几何数:, soil.geometry.duplicated().sum()) # 面积小于等于0的图斑 print(面积0:, (soil.geometry.area 0).sum()) soil.to_file(data/soil_checked.gpkg, driverGPKG)这里的逻辑是先把几何层是否有问题查清楚再查属性层。is_valid看着可有可无但后面如果用Shapely做相交或差集遇到自相交的多边形会直接抛异常duplicated()查的是完全重复的图形重叠但不是完全相同的图形要在第5章用空间连接排查。最后存成GPKG而不是回写Shapefile因为GPKG不会把中文属性名截断成10个字节字段名安全很多。2.3 裁剪到研究区掩膜提取的常见做法如果你只研究某个流域或省份不要直接按矩形框做“图框裁剪”那样会切断图斑制造出大量细碎的假边界。正确做法是用研究区边界矢量图去clip只保留与边界相交的图斑并在裁剪后按面积阈值清理碎斑。1:100万最小图斑本来就不大阈值设在1平方千米比较公允。# 读取研究区边界并统一到土壤图的坐标系 boundary gpd.read_file(data/basin.shp).to_crs(soil.crs) # 精确裁剪保留和研究区相交的所有图斑 clipped gpd.clip(soil, boundary) # 过滤掉裁剪后小于1平方公里的细碎图斑 min_area 1e6 # 平方米即1平方千米 clipped clipped[clipped.geometry.area min_area].copy() print(裁剪前:, len(soil), 裁剪后:, len(clipped)) clipped.to_file(data/soil_clipped.gpkg, driverGPKG)这段代码参数里唯一要改的就是min_area。1:100万原始制图综合时已经移除过小于规定尺寸的图斑我们裁剪后再过滤一次只是为了清掉边界上被切出来的三角碎片不是为了简化数据。阈值设得太大会把真正成片的细碎土壤类型丢失设太小又过滤不掉边界残片。1平方千米是我试下来比较稳的经验值你可以根据自己研究区大小调整。3. 读懂图里的土壤分类体系从“土类”到“亚类”别把分类码当数值算3.1 1:100万图上的分类级别中国1:100万土壤类型图在分类体系上经历过几轮演变。早期成图多采用发生学分类也就是第二次全国土壤普查建立的老分类体系层级从土纲、土类、亚类一直到土属、土种。1:100万的比例尺很难把土属和土种画出来所以图斑属性一般表达的是土类和亚类比如砖红壤、红壤、暗棕壤、黑土、黑钙土、栗钙土。后来有些数字版改用了中国土壤系统分类但图斑编号仍沿用旧名称这就埋了坑。拿到数据后我做的第一件事是找图例文件。shp包旁边如果有一个PDF或Excel格式的图例表上面能看到每个数字编码对应的土壤名称和颜色。很多网上下载的数据包把图例文件丢了只剩下一个TYCODE数字字段这种情况下宁可牺牲一点精度用图斑的中文名字段如果有的话作为语义依据也不要去猜数字编码。毕竟土壤类型数字码在行业里并没有统一国标你把“1”猜成砖红壤另一个版本里可能指赤红壤。下面这段代码是把编码和名称拉出来看一眼避免直接拿编码往下算。# 假设字段里有 TYCODE 数字码和 NAME 中文名 sample soil[[TYCODE, NAME]].drop_duplicates() print(sample.sort_values(TYCODE))如果打印出来的结果里中文名缺失那就得从图例表补齐。注意这一步没有技巧可言唯一的办法是找到同一批次的图例然后人工核对。自动化的前提是你已经拿到了可信的对照关系。3.2 把编码映射成中文名称再映射到模型分类给模型做输入时光有中文名还不够。陆面过程模型和生态模型通常需要特定分类体系比如USDA土纲或中国土壤系统分类。我的习惯是分两步映射先把数字编码映射成中文名再把中文名映射成模型需要的类别。这样中间层可检查不会出现“用了一个数字直接硬分配给模型”的黑匣子。import pandas as pd # 对照表数字编码 - 中文土壤名称示例以自己的图例为准 code2name { 1: 砖红壤, 2: 赤红壤, 3: 红壤, 4: 黄壤, 11: 暗棕壤, 12: 棕色针叶林土, } soil[SOIL_NAME] soil[TYCODE].map(code2name) # 模型分类映射中文名 - 模型需要的土纲示例 name2usda { 砖红壤: Oxisols, 红壤: Ultisols, 暗棕壤: Alfisols, 黑土: Mollisols, } soil[SOIL_USDA] soil[SOIL_NAME].map(name2usda) # 检查没映射上的类别 missing_names soil[soil[SOIL_USDA].isna()][SOIL_NAME].unique() print(未映射类别:, missing_names) soil.to_file(data/soil_mapped.gpkg, driverGPKG)这段代码里map()是查表不是计算。映射表必须自己手工建没有哪个库能自动把你的中文名翻译成USDA土纲因为不同土壤分类体系之间的对应关系本身就存在争议。我的建议是宁可把多个中文类映射到同一个模型类也不要让任何图斑落到NaN。比如“红壤”和“赤红壤”在USDA里可能都归到Ultisols那就明确写两条映射。3.3 为什么不能对编码做数值运算土壤编码是名义变量不是连续变量。1、2、3只代表“砖红壤、赤红壤、红壤”不代表它们之间有任何数量关系。常见翻车操作是拿到栅格后在计算器里写“soil_code * 10 1”想模拟亚类编码结果得到一堆不存在的类型还有人对编码取平均值统计出“平均土壤类型”这是彻底的伪科学。正确做法是重分类或条件选择。比如你只想保留暗棕壤就判断编码是否为某个集合然后输出掩膜。import numpy as np # 错误写法对编码做加法没有任何语义 # raster_after raster 1 # 正确写法只提取目标类别 target_codes [11, 12] # 暗棕壤、棕色针叶林土 mask np.isin(raster, target_codes) extracted np.where(mask, raster, 0)这里np.isin让多个类别共享一个掩膜再用np.where把非目标区域置为0。0如果不在有效编码里就可以作为背景值。这个操作常被用在模型预处理里把土壤图转成“有没有某类土壤”的0/1图它不改变语义只筛选类别。4. 栅格化与制图综合让1:100万数据变成可计算的栅格4.1 像元大小怎么选土壤类型图是矢量图但很多模型只吃栅格。矢量转栅格第一件事是定像元大小。我的判断标准很简单1:100万图上1毫米代表1千米原始图斑的细节尺度就是千米级选1千米像元最匹配。选250米虽然边界锯齿轻一点但生产不了真实细节属于“看起来精细的假信息”选2千米则会丢掉成片的小图斑后续面积统计偏差会超过10%。如果你要做的是全国陆面模拟模型网格可能是12千米或25千米我建议先转成1千米栅格后续再做聚合或面积加权不要直接从矢量图采样到12千米网格。原因是在1:100万尺度上图斑边界已经经过平滑12千米网格会把一个图斑劈成几块直接赋值会产生很大的偏差。1千米是承上启下的折中。4.2 用rasterio将矢量图斑栅格化我常用的是rasterio.features.rasterize它接受一个“几何值”的列表然后在一个目标网格上烧录数值。关键参数除了像元大小还要设定从上边界开始的转换矩阵和输出尺寸。这段代码直接可跑只要把前面准备好的soil_mapped.gpkg读进来。import geopandas as gpd import numpy as np import rasterio from rasterio.features import rasterize from rasterio.transform import from_origin soil gpd.read_file(data/soil_mapped.gpkg) # 栅格化字段必须是非空整型 soil[RASTER_CODE] soil[TYCODE].astype(uint16) # 像元大小1千米按总范围计算输出行数和列数 cell 1000 xmin, ymin, xmax, ymax soil.total_bounds width int(np.ceil((xmax - xmin) / cell)) height int(np.ceil((ymax - ymin) / cell)) transform from_origin(xmin, ymax, cell, cell) # 构建 (几何, 值) 列表剔除空几何 shapes [ (geom, code) for geom, code in zip(soil.geometry, soil[RASTER_CODE]) if geom is not None and not geom.is_empty ] # 背景/NoData用0填充 raster rasterize( shapes, out_shape(height, width), transformtransform, fill0, dtypeuint16 ) with rasterio.open( data/soil_type_1km.tif, w, driverGTiff, heightheight, widthwidth, count1, dtypeuint16, crssoil.crs, transformtransform, ) as dst: dst.write(raster, 1) dst.update_tags(CELL_SIZEstr(cell), NO_DATA0)有几个参数是关键。from_origin(xmin, ymax, cell, cell)用左上角坐标加像元大小构建仿射变换Width和Height用ceil向上取整保证右侧和下侧不会因为除以像元大小有小数而裁掉一条。fill0表示输出栅格背景值为0前提是你的土壤编码里没有0如果原图编码里有0就换成65535。整型用uint16足够因为全国土类最多几百个但如果你用USDA土纲之类的合并类也可以保持uint16。4.3 制图综合最小图斑与边界简化矢量转栅格不是制图综合。1:100万原始数据已经做完了综合我们要做的是不去破坏它。但如果你手里的源数据是1:50万或更高精度自己缩编到1:100万就要主动做一次图斑合并。常见做法是先按土壤类型字段融合把多边形内部的细小缝隙消除再按面积阈值把小于指定面积的图斑合并到相邻最长边界的图斑。面积阈值没有统一标准我一般按图上最小可读面积反推设到几十平方千米然后人工抽查最复杂的山地区域。边界简化别过度使用Douglas-Peucker简化容差设得太大会把细长河谷土壤带拉错位置。土壤边界本身受地形影响如果DEM能指导边界修正最好结合地形断裂线没有DEM就保持原样的折线不要强行平滑。制图综合是一场取舍你的目标是让图斑数量和细节符合比例尺表达能力而不是画面好看。5. 避坑5个让土壤类型图翻车的常见问题与排查5.1 编码对不上文献图例现象你按某篇论文里的“1砖红壤”去提取类别结果提取出的图斑面积和图上明显不符或者全是空值。原因1:100万图存在多个数字化版本不同单位对图例重新编号原编号系统不一定通用。解决永远先value_counts()列出所有编码再找图例表重新映射。print(soil[TYCODE].value_counts().sort_index())看到输出后把它和手头图例表逐行核对。找不到图例表时用中文名或图斑位置反推但实在推不出就放弃这个要素不要硬猜。血泪经验是猜错的编码污染的不只是一张图整个下游模型都会带着这个错。5.2 图斑之间有重叠或缝隙现象放大图斑后发现两个多边形挤在一起或中间有一道白缝栅格化后重叠区域的值取决于遍历顺序时对时错像玄学。原因分幅数字化后没有做拓扑闭合拼接处又没有精确咬合。解决先检查重叠再修复。# 自连接找出互相相交且不是自身的图斑对 intersections gpd.sjoin( soil, soil, howinner, predicateintersects ) pairs intersections[intersections.index_left ! intersections.index_right] print(重叠图斑对数:, len(pairs))这里index_left ! index_right排除的是同一个要素和自己相交。如果输出不为0说明有重叠。修复不是靠代码一行能解决的我一般在QGIS里用Fix Geometries和Delete Duplicate Geometries处理再对缝隙做Snap Geometries to Layer。缝隙如果只是视觉问题栅格化时不会造成太大影响重叠则必须处理因为rasterize遇到重叠时默认后面的值覆盖前面的值结果会随要素顺序变化。5.3 沿海地区出现大块NoData现象转出栅格后海陆交界处有一条条深色空洞看起来像数据缺失。原因原图只覆盖陆地裁剪范围又恰好包含海岸线和岛屿海岸线边界锯齿进入栅格单元后留下了空值。解决用陆地多边形先裁剪一遍或在栅格化后置NoData为0。我习惯在矢量阶段处理因为这一步同时能清理掉海洋里的岛屿碎片。# land是陆地范围多边形和soil统一坐标 land gpd.read_file(data/land.shp).to_crs(soil.crs) soil_land gpd.clip(soil, land) # 后续栅格化全部基于soil_land注意如果岛屿也是土壤图的一部分就不要用全国陆地边界简单裁剪否则把岛屿土壤全部切掉。正确做法是用原图自带的海岸线边界或者用更细的陆地行政边界。5.4 面积计算严重失真现象在ArcGIS里勾选“计算几何”得到某土类面积和统计年鉴相比差了百分之二三十。原因图层坐标系还是经纬度用平面算法按度算面积纬度越高被压缩得越狠。解决统一用Albers等积投影再算。这段代码我每轮处理都会跑一遍确保没有遗漏。import geopandas as gpd soil gpd.read_file(data/soil_clipped.gpkg) print(当前CRS:, soil.crs) if soil.crs.is_geographic: soil soil.to_crs( projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm no_defs ) soil[area_km2] soil.geometry.area / 1e6 print(soil.groupby(SOIL_NAME)[area_km2].sum())这里的检查点是is_geographic如果坐标为经纬度就立即投影。1:100万制图规范里Albers是标准投影但不能保证所有数字版都遵守。面积统计前先看CRS养成习惯后能避免最可笑的翻车。5.5 分类体系混用现象一张全国图里“暗棕壤”在不同省区的图斑定义不一样属性字段相同但内容体系不同。原因制图时间跨度长早期图幅用发生学分类后期更新图幅采用中国土壤系统分类合并时没有统一。解决先识别每块图斑所属体系再统一映射。识别方法是看伴随的属性字段如果只有老土类名多半是发生学分类如果出现“潜育”等诊断层命名特征可能是系统分类。统一映射是持久工作建议把对照表写进CSV并版本管理。原图例 原体系 统一后 暗棕壤 发生学 淋溶土 暗棕壤 系统分类 暗沃冷凉淋溶土这张示意表不严格但说明了思路映射关系要能让审查者看到原类别和目标类别而不是用一个数字糊弄过去。分类体系混用最隐蔽它不影响几何形状却直接影响模型参数和解释结果。宁可映射得粗不要映射得假。6. 进阶验证用土壤属性和地形数据反向校验这张图的可靠性6.1 用剖面样点做抽检如果你有全国第二次土壤普查的剖面点位数据或者发表的文献里有土壤类型点可以做一次抽检。把点落到图上提取所在图斑的类型和点位记录比对算出一致率。一致率低于70%说明数据版本和你的分类映射存在系统性问题需要回到第3章重新建映射。我一般抽300个点分东南西北四个区避免所有点集中在同一个土类。6.2 用地形因子交叉验证土壤类型和地形高度相关比如暗棕壤常见于东北山地砖红壤分布在低海拔热带丘陵。如果栅格化后的图上暗棕壤跑到平原绿洲里大概率是分类或几何错误。用DEM叠加统计每个土壤类型的海拔范围能快速定位异常图斑。import rasterio import numpy as np with rasterio.open(data/dem_1km.tif) as src: dem src.read(1) soil_ras rasterio.open(data/soil_type_1km.tif).read(1) for code in np.unique(soil_ras): mask soil_ras code if mask.sum() 0: valid dem[mask] valid valid[valid ! -9999] # 排除DEM NoData print(code, 海拔均值:, round(valid.mean(), 1), 最低:, round(valid.min(), 1), 最高:, round(valid.max(), 1))这段代码按栅格类别统计海拔输出后对照常见土壤地理分布规律。如果某个类别出现两个判若鸿沟的海拔区间比如一段在海拔50米一段在海拔3000米说明这个编码里混入了不同土类需要回到矢量图分割检查。地形验证不能证明某类土壤一定对但能快速证伪明显错误成本很低。6.3 发布前的重分类一致性检查最后一个技巧是在发布前把结果重分类成三层结构农业土壤、森林土壤、高山土壤等和现有土地利用数据交叉统计。如果某个森林土壤大范围落在建设用地上边界可能错了。我习惯把每一步中间产物都输出成GPKG和TIF并打上日期这样就算两个月后发现错误也能用后悔药回到正确的处理步骤。中国土壤类型图1:100万本质上是宏观决策的底图使用它时保持尺度敬畏比追求高精度更重要。希望这些踩坑记录能帮到你。本文还有配套的精品资源点击获取