简介:这份中国1:100万土壤数据集,源自FAO与IIASA共同构建的世界和谐土壤数据库(HWSD)1.1版本,覆盖全国范围,可为地理建模、生态农业分区、粮食安全与气候变化研究提供关键输入参数,适合GIS、土壤学及相关领域的科研人员和学生作为基础底图或模型驱动数据使用。压缩包共含29个文件,整体大小约7.78MB,主要文件类型包括ArcGIS栅格数据(adf/aux)、矢量属性表与空间数据库(dbf/mdb)、说明文档(doc/pdf/txt)以及索引、日志等辅助文件,其中包含可查询的土壤属性数据库和HWSD官方说明,便于理解字段含义并直接导入常用GIS平台。已有6424人学习下载,数据组织规范,拿到后可直接进行制图、统计或作为生态、农业模型的输入参数,省去自行下载和处理多源数据的环节,适合需要规范全国尺度土壤信息的用户快速上手。
1. 中国1:100万土壤数据集,为什么我今天还在教你怎么用
中国1:100万土壤数据集,说白了就是一张把全国土壤类型画成图斑的矢量底图,比例尺的“100万”既是它的精度标签,也是它的能力边界。这套数据脱胎于历史上一次全国性土壤调查的成果整理,后来被多处科学数据平台收录并再分发,今天你在公开渠道搜这个标题,拿到的通常是一组面状图层加一份属性表。对它抱有正确预期的人会把它当区域分析的宝贝;指望它给出地块级答案的人会觉得被坑了。这篇文章就是写给前者,以及那些刚拿到数据还不知道怎么下手的GIS从业者。我会按实际处理的流程讲清楚:比例尺怎么读、属性表怎么清洗、投影怎么选、栅格化分辨率怎么定,以及最容易翻车的几个地方。
2. 读懂1:100万土壤数据集的地图语言:比例尺、分类体系与图斑逻辑
2.1 比例尺背后的精度账:图上1毫米相当于地面1公里
1:100万的比例尺意味着图上1毫米对应地面1公里。这不仅仅是一个数学换算,它直接决定了这套数据能表达什么、不能表达什么。制图综合是比例尺的天然产物:一份1:100万土壤图上,那些面积只有几公顷的小土斑通常已经被合并进相邻图斑,或者干脆被删掉。能稳定表达出来的最小图斑,往往在几平方公里到几十平方公里这个量级。所以这套数据的服务对象永远是区域格局问题:某流域的土壤类型组合、省级生态区划的底图、全国尺度的土类分布对比。
如果你要做的是某个县、某个园区的地块决策,请换更大比例尺的土壤调查数据,而不是试图在1:100万图层里硬挖细节。这不算数据缺陷,而是比例尺决定的“分辨率”。很多新人拿到数据后抱怨图斑太粗,本质上是预期错位。我一般会在项目启动时先给协作方打个预防针:这套数据回答“这片区域以什么土壤为主”,回答不了“这块地到底是什么土”。
| 应用场景 | 是否合适 | 判断依据 |
|---|---|---|
| 省级、流域级土壤类型分布制图 | 合适 | 比例尺匹配,图斑粒度恰当 |
| 全国生态分区与农业区划 | 合适 | 覆盖全域、分类口径相对统一 |
| 县级土壤资源评价 | 勉强 | 只能作背景参考,不能单独支撑决策 |
| 地块级精准农业 | 不合适 | 最小图斑远大于地块单元 |
2.2 两套分类体系:发生分类与系统分类,图斑属性怎么读
土壤图的核心不是边界,而是每个图斑属性里挂的土壤类型。1:100万数据最常见的属性组织方式是“代码+名称”,代码位数通常能看出分类层级。早期成果多采用土壤发生分类,按成土条件和地带性划分土类,常见的如黑钙土、栗钙土、红壤、黄壤,图斑属性里一般有土类、亚类、土属、土种四级;后续整理版本有的会同步挂接系统分类的诊断层信息或系统分类名称。
两种分类不是能简单互查的字典,因为划分依据不同:发生分类看重成土过程与地带性,系统分类看重诊断层与诊断特性。实际使用中我的建议是:不要试图在代码层面做自动转换,选定一套作为主分类,另一套作为辅助描述。比如做生态分区时,以土类为统计口径;做地块解释时,以系统分类的土纲为参考。属性字段名在不同版本里差异很大,有的叫“SOIL_CODE”,有的叫“T_CODE”“土类代码”,拿到手第一件事不是分析,而是先把随数据附带的说明文档通读一遍,搞清楚每一列到底挂的是哪一级分类。
2.3 拿到手的文件长什么样:图层组织、属性项与坐标底细
公开渠道拿到的数据,文件组织常见两种:一套shapefile多图层,或者一个文件地理数据库。面图层是主体,属性项通常有图斑编号、类型代码、类型名称、面积、周长,有的版本还带行政代码字段。辅助图层可能包括境界线、主要河流,这些在制图时有用,在空间分析时通常可以关掉。
最容易出问题的就是坐标底细。很多版本以经纬度存储且坐标参考信息残缺,个别版本坐标已经做过投影却仍标为地理坐标。这个“黑匣子”拿到手不要急着投影,先原样加载,叠上行政区划边界看空间范围对不对。范围对,再做后续;范围错,先解决坐标问题。坐标问题属于那种不查不知道、一查吓一跳的类型,等分析做到一半才发现错位,返工成本极高。我的习惯是把原始目录设为只读,所有操作都在副本上做,这样坐标、属性折腾坏了还能回退。
3. 把数据搬进自己的分析库:预处理、坐标转换与属性清洗
3.1 从原始数据到干净图层:先备份、再标记坐标、后投影
拿到任何土壤数据,我建议都按这个顺序做:原样读取、检查坐标参考、标记坐标、投影转换、另存新文件。不要在原始文件上做任何写操作,尤其是dissolve和保存,历史数据的几何偶尔有损坏,一旦在原文件上执行过再保存,想恢复就只能重新下载。
import geopandas as gpd from pyproj import CRS # 第一步:原样读取,先观察再动手 gdf = gpd.read_file("soil_1m_raw/soil_polygon.shp") print(gdf.head()) print(gdf.crs) # 查看是否带坐标参考信息 print(gdf.geometry.type.value_counts()) # 确认要素类型全是面 # 第二步:如果crs为空,按数据说明手工指定地理坐标 if gdf.crs is None: gdf = gdf.set_crs(CRS.from_epsg(4326)) # 第三步:投影到等积圆锥并另存,原始文件不动 gdf_proj = gdf.to_crs( CRS.from_proj( "+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 " "+lon_0=105 +datum=WGS84 +units=m" ) ) gdf_proj.to_file("soil_1m_albers.gpkg", layer="soil_polygon", driver="GPKG")这里的核心是区分set_crs和to_crs:set_crs只是给数据打上坐标参考标签,不移动任何点位;to_crs才真正做坐标转换。初学者最常见的错误就是拿set_crs当转换用,结果所有坐标全错。投影参数选用Albers等积圆锥,中央经线105°E、双标准纬线25°N和47°N,是覆盖全国范围制图统计的常用设置,面积变形小,适合后续做面积量算。
3.2 属性清洗与图斑溶解:把同义代码合并成一张完整分布
1:100万土壤数据的属性表普遍存在两类问题:代码列混入空格和短横,以及同一土壤类型存在多个写法。前者源于历史表录入格式不统一,后者源于数据多次拼接。清洗要分两步走,先规范代码,再按标准代码溶解图斑。
import geopandas as gpd import pandas as pd gdf = gpd.read_file("soil_1m_albers.gpkg") # 第一步:代码列统一转字符串并去掉首尾空格 gdf["TYPE_CODE"] = gdf["TYPE_CODE"].astype(str).str.strip() # 第二步:过滤无效值,避免空代码干扰统计 gdf = gdf[~gdf["TYPE_CODE"].isin(["", "-", "999", "0", "nan"])].copy() # 第三步:建立别名映射,把历史版本的同义代码归并到标准码 code_map = { "101": "100", "0101": "100", "A101": "100", "205": "200", } gdf["CODE_STD"] = gdf["TYPE_CODE"].map(code_map).fillna(gdf["TYPE_CODE"]) # 第四步:只保留分类代码和几何,按标准码溶解 gdf = gdf[["CODE_STD", "geometry"]].copy() gdf_dissolved = gdf.dissolve(by="CODE_STD").reset_index() # 第五步:在等积投影下计算面积,单位是平方米 gdf_dissolved["AREA_KM2"] = gdf_dissolved.geometry.area / 1e6 print(gdf_dissolved[["CODE_STD", "AREA_KM2"]].head())注意astype(str)这一步很关键。土壤代码如果按数值读取,前导0会被丢掉,比如“0101”变成“101”,再和映射表对照就永远对不上。先转字符串再处理,能避免这类“看起来一样的代码实际不匹配”的问题。dissolve按CODE_STD分组,把同一类型的碎图斑合并成连续面,合并后面积统计就清爽多了。
3.3 全国尺度的投影参数:Albers等积圆锥怎么设不踩坑
投影选型的原则很简单:做面积统计用等积投影,做距离分析用等距投影,做形态对比用等角投影。土壤数据90%的分析落在面积统计上,所以Albers等积圆锥是默认选择。
| 参数 | 全国尺度推荐值 | 说明 |
|---|---|---|
| 投影 | Albers等积圆锥 | 面积变形小,适合区域统计 |
| 中央经线 | 105°E | 全国尺度居中;分区分析改为区域中心经线 |
| 标准纬线 | 25°N、47°N | 控制南北向面积变形 |
| 椭球/基准 | WGS84或CGCS2000 | 与原始数据保持一致,不要混用 |
做分区研究时,这套参数要按研究区重新调整。比如只做某河流域,中央经线取流域中心经度,标准纬线取流域南北边界纬度,能进一步压缩面积误差。参数不是死的,但有一个底线:不要用经纬度坐标直接算面积,那会把1度经度当1度纬度处理,结果差出好几倍。换投影后,记得用gdf_proj.geometry.area检查一下面积量级是否合理,常识不够时数字会替你报警。
4. 用这套数据出活:栅格化、面积统计与边界取舍
4.1 矢量转栅格:分辨率为什么不能低于1公里
做模型输入时,经常要把土壤矢量图转成栅格。分辨率的取值逻辑不是“越细越好”,而是“匹配数据的真实精度”。1:100万土壤数据的制图精度大约对应地面1公里,栅格分辨率取1公里是合理的;取250米甚至30米,得到的是“看起来精细、实际全是相邻像元复制”的假象,还平白增加存储和计算开销。如果研究区是局部区域,可以与气象格网对齐,取和格网一致的分辨率。
import rasterio from rasterio.features import rasterize from rasterio.transform import from_origin import numpy as np # 用图斑范围和代码做栅格化,分辨率取0.008333度,约1公里 res = 0.008333 bounds = gdf_dissolved.total_bounds # (minx, miny, maxx, maxy) width = int((bounds[2] - bounds[0]) / res) height = int((bounds[3] - bounds[1]) / res) transform = from_origin(bounds[0], bounds[3], res, res) shapes = [ (geom, int(code)) for geom, code in zip(gdf_dissolved.geometry, gdf_dissolved["CODE_STD"]) ] meta = { "driver": "GTiff", "height": height, "width": width, "count": 1, "dtype": "int32", "crs": gdf_dissolved.crs, "transform": transform, } with rasterio.open("soil_code_1km.tif", "w", **meta) as dst: dst.write(rasterize(shapes, out_shape=(height, width), transform=transform, fill=0), 1)from_origin用左上角坐标计算仿射变换,原点取自total_bounds的左上角,这样栅格范围能完整覆盖所有图斑。fill=0表示无数据区域填0,后续分析中要记得把0当掩膜处理。代码列必须转成整数才能写入int32栅格,所以前面清洗时的CODE_STD这里要做int()转换。如果原始代码里还有非数字字符,这一步会报错,正好借机回头检查清洗是否彻底。
| 目标用途 | 推荐分辨率 | 说明 |
|---|---|---|
| 全国生态分区分析 | 1公里或30弧秒 | 与数据精度匹配 |
| 流域水文模型输入 | 250~500米 | 与气象格网对齐使用 |
| 省级土壤制图 | 不需要栅格化 | 矢量出图更合适 |
4.2 面积统计与占比速算:让groupby替你干活
面积统计是这套数据最高频的操作。很多人在Excel里手动算,其实投影后的GeoDataFrame直接用groupby就能出结果。土壤类型面积统计有一个原则:必须先投影、后统计。地理坐标下算出来的面积没有物理意义,只属于“能算但不能信”的范畴。
import pandas as pd # gdf_dissolved 已经在 Albers 投影下,且带 AREA_KM2 列 summary = ( gdf_dissolved.groupby("CODE_STD")["AREA_KM2"] .sum() .sort_values(ascending=False) ) share = summary / summary.sum() result = pd.DataFrame({"面积(km²)": summary, "占比": share.round(4)}) print(result)groupby按标准代码聚合所有图斑的面积,sort_values让最大的类型排在最前,一眼看出区域主导土壤。加一列占比,做报告时直接引用。这里有个细节:dissolve之后每个类型的图斑数量大幅减少,比原始碎图斑逐条统计快一个数量级。统计结果出来后,我习惯把总面积和行政区域总面积做一次粗略对比,偏差在合理范围内再往下走。
4.3 接边与碎图斑:制图综合的取舍逻辑
处理这套数据时,最考验工程判断的是碎图斑去留。1:100万原始数据里,图斑大小差异巨大,大的覆盖上千平方公里,小的不足1平方公里。做分析时,小图斑不仅拖慢计算,还会在栅格化后产生大量零碎像元,影响后续分类统计的稳定性。
我一般用两个参数控制制图综合:简化容差和合并阈值。简化容差控制在100米以内,主要去除边界上的微小锯齿,不会明显改变图斑形态;合并阈值设为1平方公里,小于这个面积的图斑并入相邻最大图斑。这样的取舍不会改变区域尺度上的类型占比,但能让下游计算的稳定性明显改善。如果你用的是PostGIS或QGIS,处理逻辑一样:先按面积筛选,再按相邻关系合并。关键在于把阈值写进数据处理记录,报告中说明“本次分析剔除了小于1平方公里的图斑”,后人复核时就不会觉得数据被篡改。
5. 避坑清单:用这套土壤数据最容易翻车的5个地方
5.1 面积统计和行政区划对不上,别急着怀疑数据
现象:按这套数据统计出的某区域土壤总面积,和统计年鉴里的行政区域面积差出10%以上。原因:1:100万制图综合时,大量小图斑被并入相邻图斑或直接删除,区域边界附近还有大量未分类图斑被归为“其他”;另一方面,统计年鉴面积是经过精确量算的行政面积,两套数字口径不同。解决:固定投影和统计口径,只做同数据内的相对比较。要和其他数据对比时,先统一到同一行政边界,再计算土壤类型的占比,而不是比较绝对面积。“对不上”是常态,不等于数据坏了。
5.2 代码列转数字就报错,先查空格和短横
现象:int(code)直接抛ValueError,或者某些图斑的代码在结果里变成NaN。原因:历史属性表录入不规范,“101”写成“101 ”或“10-1”,空格和短横混进代码列。解决:先astype(str)再str.strip(),过滤空字符串和占位符“-”“999”“0”。这个清洗步骤必须在任何映射、统计之前完成,否则错误会一路传导到最后的统计表里。
5.3 同一个土类在图层里有两个代码,合并前先建对照
现象:同一片分布的土壤,图斑A的代码是“101”,图斑B的代码是“0101”,属性名称完全一样,但统计时被当成两个类型。原因:数据来自多次拼接,早期版本和后整理版本混在一起,代码规则没对齐。解决:建立代码对照表,把别名映射到标准码。映射表写成一个单独的CSV或字典,不要直接在代码里写死,方便复核和复用。
5.4 叠加行政边界时整体错位,先怀疑坐标系
现象:土壤图斑和行政区划边界套在一起时,整体偏移几公里到几十公里,南部偏、北部也偏。原因:数据标注的是地理坐标,但坐标参考信息缺失或标注错误;也可能是把投影坐标当成地理坐标加载。解决:先原样加载,叠加一个已知坐标系正确的边界要素查看范围是否吻合。范围明显不对时,检查图层属性里的坐标参考信息,必要时用已知控制点做空间校正。“先查坐标再做投影”这个顺序不能颠倒。
5.5 栅格化后细碎图斑全部消失
现象:矢量图里有几十个小图斑,转成栅格后一个都不见,只剩下大片连续类型。原因:栅格分辨率太粗,或者图斑面积小于一个像元;rasterize默认把面积不足的图斑直接丢弃。解决:分辨率取1公里时,小于1平方公里的图斑就会被吞掉。要么接受这个精度损失并在文档中说明,要么先用合并阈值把碎图斑并入相邻图斑,这样至少保留类型归属,而不是凭空消失。
6. 把老数据用活:属性重编码与多源校验的进阶习惯
拿到1:100万土壤数据,多数场景不会直接用原始土类代码,而是按项目需求重新归并。比如做生态分区时,把几十个土属归并成森林土、草原土、荒漠土、水成土、盐成土几大类;做水文模型时,按土壤质地组重编码成产流参数。重编码规则应该写成外部文件,而不是散落在代码里,这样换数据、换分类体系时不用重写逻辑。
import pandas as pd import geopandas as gpd def recode_soil(series, rule_csv): """按规则表将原始土壤代码映射到项目分类代码""" rule = pd.read_csv(rule_csv) mapping = dict(zip(rule["old_code"], rule["new_code"])) return series.map(mapping).fillna(series) gdf = gpd.read_file("soil_1m_albers.gpkg") gdf["GROUP_CODE"] = recode_soil(gdf["CODE_STD"], "soil_group_rule.csv")规则表至少保留old_code、new_code、new_name三列,old_code未出现在表里的代码用fillna保留原值,避免数据静默丢失。重编码完成后,用野外样点或区域文献做一次校验:统计每个重分类组的面积和占比,和已知的区域土壤分布特征对照,偏差过大时优先检查映射表,而不是怀疑数据。
我早年做某区域生态评价时,拿着这套数据直接出图,边界错位被合作方当场指出来,后来每次拿到土壤数据都先叠边界再动刀。数据老没关系,流程不乱就不会翻车。希望帮到你。
本文还有配套的精品资源,点击获取