简介上海市土壤水文分组HSG高精度栅格数据是基于USDA曲线数CN方法估算降雨径流的关键基础数据适用于SWAT水文模型建模、城市排水规划及流域水文分析。数据源自HYSOGs250m与soilGrids250m系统按省级行政区裁剪为上海市范围WGS84地理坐标系空间分辨率约250m提供A、B、C、D四类径流潜力等级并包含湿润土壤双重HSG分类信息。压缩包大小3.46MB文件总数0类型明细暂无数据可用GIS软件直接加载分析。目前已有212人学习适合水文、环境、GIS领域的研究者与工程师可直接用于计算径流曲线数、划分水文响应单元及开展区域水资源评估。1. 上海市土壤水文分组高精度栅格给产流模型一张能直接查表的底图做上海市域内涝模拟或海绵城市绩效评估的人迟早会被同一个问题卡住SCS-CN 查表需要土壤水文分组HSGA/B/C/D 四档而手头常见的是二普扫描图、HWSD 的 1 km 网格或 SoilGrids 的 250 m 插值。上海市土壤水文分组高精度栅格就是把上海地表土壤按入渗能力归入四档做成覆盖全市的 10 m 或 30 m 整幅栅格让 SWMM、InfoWorks ICM、HEC-HMS 在每个像元上直接完成 HSG 到 CN 的映射。难在上海几乎每一项地表条件都在抬升难度高差几米、地下水位贴着地表、硬化面把土体切得支离破碎而 1:50 万土壤图投影到 10 m 栅格本身就是一种假精细。下面按判定、制作、验证、对接四步把链路拆开。2. A 到 D 四档怎么定USDA 水文土壤分组的判定参数与上海土壤归属逻辑2.1 水文土壤分组不是土壤分类而是入渗能力分级中国土壤发生分类里的水稻土、潮土、滨海盐土回答的是土壤怎么来的USDA 水文土壤分组回答的是另一件事长时间湿润、无冻土、无地表封盖条件下土壤在强降雨里能吞下多少水。NRCS 在 National Engineering Handbook Part 630 中按入渗能力把土壤分成 A、B、C、D 四组A 组通常为砂质或砾质入渗极快B 组中等C 组偏慢D 组极慢。D 组有一条补充定义对上海特别致命即使质地不粘只要下挖 60 cm 内遇到持续地下水位或 50 cm 内遇到致密不透水层也判为 D。上海大部分地区地下水位在 0.51.5 m雨季贴着地表这一条直接锁定了大片 D 组。2.2 先看质地再看 Ks最后用排水条件兜底判定流程我习惯按三步走。第一步查土种的质地剖面按 USDA 质地三角图映射出候选组第二步用饱和导水率 Ks 复核实验室定水头法、现场双环入渗仪都可以第三步用地下水埋深和不透水层深度做一票否决。三者对应关系如下HSG典型质地稳态入渗率 / Ks (mm/h)排水与水位条件A砂土、壤质砂土大于 7.6排水迅速无滞水迹象B砂质壤土、壤土、粉砂壤土3.8 ~ 7.6中等排水C粘壤土、砂质粘壤土、粉砂质粘壤土1.3 ~ 3.8有弱滞水层排水较慢D粘土、重粘土或任意质地但 60 cm 内有水位小于 1.3排水差或地下水位高表里的 Ks 取的是 NRCS 常用参考区间。实际判定时看的是整层剖面中最差的那个层次不是表层——上海水稻土的耕作层透水性尚可真正决定分组的是犁底层。2.3 上海主要土种的归属判断按土类分布看南汇、奉贤沿海的滨海盐土粉砂质粘壤土加常年高水位归 D青浦、松江、金山的青紫泥、青黄土等水稻土犁底层紧实、下有潜育层归 D嘉定、宝山及黄浦江以西的灰潮土粉砂壤土到粘壤土、水位 12 m多为 C崇明东滩、长兴岛新近沉积的砂质潮土局部到 B砂性强的可到 A。老城区内的回填土和建筑渣土扰动层属性无法从二普图获得单独编码为 9 类扰动土不参与常规查表。具体归属如下土类典型分布控制性剖面特征HSG滨海盐土南汇、奉贤沿海粉砂质粘壤土水位常年小于 1 mD青紫泥 / 青黄土水稻土青浦、松江、金山粘壤土加紧实犁底层加潜育层D灰潮土嘉定、宝山、黄浦江以西粉砂壤土至粘壤土水位 12 mC砂质潮土崇明、长兴岛、古河道带砂质壤土剖面疏松A / B2.4 从粒径或质地名称自动落组字典映射就够如果数据源给的是粘粒、砂粒、粉粒百分数先按 USDA 三角图定质地类再按字典映射落到四档TEXTURE_HSG { sand: A, loamy sand: A, sandy loam: B, loam: B, silt loam: B, silt: B, sandy clay loam: C, clay loam: C, silty clay loam: C, sandy clay: D, silty clay: D, clay: D, } def texture_to_hsg(texture: str) - str: return TEXTURE_HSG.get(texture.strip().lower(), None)映射本身很朴素功夫在前处理把《上海土壤》里每个土属的典型质地剖面整理成 CSV逐层读入取剖面中最差一层的质地做最终归属。老资料里中壤重壤这类卡庆斯基制名称要先转换成 USDA 制的砂粒、粉粒、粘粒百分数再进三角图两套制式的质地边界并不重合。粒径数据充足时用 Saxton-Rawls pedotransfer 函数把粒径、容重、有机质换算成 Ks再按 2.2 节的阈值分组比纯质地映射更贴实测但该函数标定数据主要来自美国土壤对上海高粉粒、高盐基的冲积土外推时建议留 20% 左右的余量。注意城市绿地里的客土绿化工程换填的种植土资料上常标为砂壤土实际压实后入渗率往往跌到 C 组遇到这种情况以现场探测为准。3. 高精度栅格怎么造数据源分工、投影统一与 rasterio 栅格化参数3.1 先搞清楚高精度到底来自哪里1:50 万上海市第二次土壤普查数字化图的精度在制图综合层面最小图斑对应的地面范围远大于 10 m 像元HWSD 是 30 弧秒网格SoilGrids 是 250 m。直接拿任何一个单一来源做 10 m 输出都只是把粗边界切细叫高精度名不副实。常见做法是把来源按职责拆开让各自干各自擅长的事数据来源尺度在高精度栅格里的职责1:50 万上海土壤图二普数字化版上海农科院、土地档案整理图斑级主边界与土种属性《上海土壤》土种志1990 年代公开出版属性表质地、剖面、水位逐层记录ALOS PALSAR 12.5 m DEMASF / OpenTopography 分发12.5 m河漫滩、古河道、湖沼平原地貌细分上海土地利用 10 m 产品测绘遥感解译10 m城市扰动区标记、水体掩膜SoilGrids 250 mISRIC250 m粒径空间插值的辅助背景这样组合出来的栅格边界精度来自土壤图内部细分来自地形与用地信息才算真正对得起高精度三个字。3.2 投影与像元对齐统一到本地投影坐标系上海市常用 CGCS2000 基准、3° 带高斯-克吕格投影中央经线取 121.5°E单位为米。生产上第一步是把所有矢量重投影到与模板栅格一致的坐标系模板一般取土地利用栅格或 DEM。像元对齐有三条硬约束像元尺寸取 10 m 的整数倍像元原点左上角与模板完全一致行列数不能四舍五入否则整幅栅格会沿东北方向漂移半个像元与路网叠合时错位肉眼可见。最省事的做法是让模板栅格先定义好 transform后续所有栅格化和重采样都以它为唯一基准。3.3 rasterize 参数与输出写法import geopandas as gpd import numpy as np import rasterio from rasterio.features import rasterize # hsg_code: 1A, 2B, 3C, 4D, 9扰动土, -9999无数据 gdf gpd.read_file(sh_soil_50w_hsg.shp) with rasterio.open(sh_lu_10m_template.tif) as ref: gdf gdf.to_crs(ref.crs) # 统一坐标系 transform ref.transform out_shape (ref.height, ref.width) shapes [(geom, int(code)) for geom, code in zip(gdf.geometry, gdf[hsg_code])] raster rasterize( shapes, out_shapeout_shape, transformtransform, fill-9999, dtypenp.int16, all_touchedFalse, ) profile ref.profile.copy() profile.update(dtypenp.int16, count1, nodata-9999, compressdeflate) with rasterio.open(sh_hsg_10m_raw.tif, w, **profile) as dst: dst.write(raster, 1)rasterize 的三个参数值得单独说。all_touchedFalse 表示只有像元中心落在图斑内才赋值避免沿边界多占一圈像元对 10 m 栅格和细碎图斑先设 False若边界处出现大量空洞再改 True。dtype 用 int16 不是因为分组只有 4 个值而是要为 9扰动土和 -9999无数据留编码空间同时保持与土地利用栅格一致的整型语义后续 CN 查表时布尔索引的速度也更接近 C 级。fill-9999 先铺底栅格化后未覆盖区域全部是 -9999后面统一用上海市域掩膜裁掉不加掩膜的话沿海潮间带和江面里的沉积物会被当成真实土壤算进 D 组。3.4 边界修正与孔洞清理gdal_sieve 的用法原始栅格化结果有两个典型毛病河岸边土壤图与水体相交产生锯齿状伪像元细碎图斑在 all_touchedFalse 下变成孔洞。处理分两步。先做水体掩膜把黄浦江、长江口、淀山湖等水体多边形栅格化为 -9999再做连通域清理用 GDAL 自带的 sieve 工具去掉小于阈值的碎斑gdal_sieve.py -st 20 -8 sh_hsg_10m_raw.tif sh_hsg_10m_sieve.tif-st 20 表示小于 20 个像元的连通域并入周边最大类-8 表示八邻域连通判定。阈值按碎斑分布调上海市区地块被道路切得很碎20 太小会保留地块内部的图斑噪声我一般给到 3050。最后把城市建设用地边界内的像元用土地利用图覆写为 9 类扰动土避免老图斑在建成区复活。sieve 只改空间连续性不改变大类内部的真实过渡对 10 m 栅格是安全的。提示如果下游是 SWMM 这类子汇水区模型栅格分辨率取子汇水区平均尺寸的 1/4 到 1/2 即可过细只增加重采样耗时不会提高 CN 精度。4. 上海土壤水文分组栅格的精度验证样点布局、混淆矩阵与边界效应排查4.1 样点布局分层抽样比均匀撒点可靠验证样本按 HSG 和地貌单元双重分层抽取。每个 HSG 组不少于 20 个点全市 80120 个点足以支撑四类混淆矩阵。按二项分布估算期望总体精度 80%、允许误差 10%、置信水平 95% 时最少约 62 个点和这个范围吻合。点位约束有三条距图斑边界至少一个像元避开道路硬化面避开已竣工基坑。现场每点挖 60 cm 剖面读质地与层次再用双环入渗仪测稳定入渗率剖面判分组入渗率用于解释偏差来源而不是直接改判——单点入渗受压实、根系和裂隙影响很大一次读数不能代表整个图斑。4.2 混淆矩阵、总体精度与 Kappa把实测编码与栅格提取编码对齐后用 sklearn 一次算出三个指标import numpy as np from sklearn.metrics import confusion_matrix, cohen_kappa_score obs np.array([1, 1, 2, 2, 3, 3, 3, 4, 4, 4]) # 实测 HSG pred np.array([1, 2, 2, 2, 3, 4, 3, 4, 4, 4]) # 栅格提取 cm confusion_matrix(obs, pred, labels[1, 2, 3, 4]) oa np.trace(cm) / cm.sum() kappa cohen_kappa_score(obs, pred) producer np.diag(cm) / cm.sum(axis0) # 生产者精度 user np.diag(cm) / cm.sum(axis1) # 用户精度一组合格的 90 点验证结果如下预测 A预测 B预测 C预测 D合计用户精度实测 A162101984.2%实测 B318212475.0%实测 C022022483.3%实测 D001222395.7%合计1922242590生产者精度84.2%81.8%83.3%88.0%总体精度 84.4%Kappa 0.79。这个数放在 1:50 万源图衍生的 10 m 产品里是合格的真正的诊断信息在混淆结构上B 与 C 互相错分占了大头说明质地过渡带被硬切成两档的边界误差远大于整块分错类的误差。Kappa 低于 0.8 时不要急着改栅格先回到 2.3 节的剖面最差层判定逻辑复查原始属性表。4.3 三类最容易出错的边界河岸带是重灾区。黄浦江和吴淞江沿岸的砂质透镜体多是古河道摆动留下的土壤图上没有要用 DEM 的河漫滩相对高度和影像纹理补。做法是取距水系 100 m 缓冲带叠加相对高程小于 2 m 的区域把原判 C 的像元复核为 B。第二类是田埂和大棚菜地犁底层和压实层让剖面整体下移半档这类地类面积逐年缩小但残余图斑对 CN 的影响不小。第三类是建设用地边缘回填土与原地表犬牙交错栅格上表现为 9 类像元包围 D 组的椒盐噪声。直接做 majority filter 会抹掉真实信息我一般保留原值在属性里加 flag 字段模型侧对 9 类单独设 CN而不是硬并进 C 或 D。5. 把 HSG 栅格接进上海 SCS-CN 径流模拟CN 查表合成、AMC 修正与反推验证5.1 TR-55 的 HSG-CN 查表关系HSG 单独不产生水文意义意义在 SCS-CN 的查表环节。TR-55 给每个土地利用类别配了四组 CN 默认值上海的常见类别如下土地利用上海常见类别ABCD绿地 / 公园良好状态39617480低密度住宅1/3 英亩61758387高密度住宅1/8 英亩77859092商业 / 工业85% 不透水89929495旱作农田常规耕作67758387硬化路面 / 停车场98989898这张表的工程潜台词是土地利用对 CN 的影响远大于 HSG 分组。上海市区以高密度建成区为主HSG 的误差只有在绿地、农田这类透水地类上才会放大所以验证资源应当优先投到郊野和蓝绿空间。5.2 栅格级 CN 合成的 numpy 写法把 HSG 栅格与同期土地利用栅格对齐到同一像元网格后按查表逐项赋值import numpy as np import rasterio CN_TABLE { (1, 1): 39, (2, 1): 61, (3, 1): 74, (4, 1): 80, # 绿地 (1, 2): 89, (2, 2): 92, (3, 2): 94, (4, 2): 95, # 商服/工业 (1, 3): 98, (2, 3): 98, (3, 3): 98, (4, 3): 98, # 不透水 (1, 4): 67, (2, 4): 75, (3, 4): 83, (4, 4): 87, # 旱作农田 } with rasterio.open(sh_hsg_10m.tif) as hs, \ rasterio.open(sh_lu_10m.tif) as lu: hsg hs.read(1) land lu.read(1) cn np.full(hsg.shape, -9999, dtypenp.float32) for (h, l), val in CN_TABLE.items(): cn[(hsg h) (land l)] val profile hs.profile.copy() profile.update(dtypenp.float32, nodata-9999, compressdeflate) with rasterio.open(sh_cn_10m.tif, w, **profile) as dst: dst.write(cn, 1)这里有一个覆盖顺序的坑字典遍历的先后决定同名像元谁胜出。上面这张表各键互斥所以安全若引入水域、湿地等多条件叠加的类别必须把优先级最高的类别放最后赋值或用 np.select 按条件数组顺序求值。5.3 三个必调参数与一个反推验证技巧第一个必调参数是前期土壤湿度条件 AMC。上海梅雨和台风季前期土壤接近饱和AMC-III 的 CN 比 AMC-II 高 515只出一版栅格会系统性低估产流。常见做法是准备两版 CN 栅格用前 5 日降雨量判断当次取值AMC-II 到 AMC-III 的换算用 SWAT 近似式 CN3 CN2 * exp(0.00673 * (100 - CN2))。第二个是初损率 λ。TR-55 默认 0.2上海城镇流域的短历时降雨用 0.05 与实测过程线拟合更好S 的换算按 Hawkins 的 S05 1.33 * S02^1.15英寸制重算不能直接套 0.2 的 S 值。第三个是聚合重采样。10 m 栅格聚合到子汇水区时只能取众数或面积加权平均双线性插值会把四类整数插出 3.7 这样的伪类。注意凡 CN_TABLE 未覆盖的像元nodata、9 类扰动土最后都会落在 -9999 上模型侧要按无效值剔除不要用 0 兜底0 会被当成 CN100 参与计算。验证技巧挑上海近年来 35 场有完整雨量与出口流量记录的短历时降雨反推流域平均 CN再与同流域栅格 CN 的面积加权均值对比。偏差在 ±5 以内说明 HSG 栅格与土地利用编码的配合是自洽的若偏差方向一致且普遍偏大优先检查土地利用的 CN 取值和不透水率假设而不是回头改土壤分组。本文还有配套的精品资源点击获取 
长河君 · 2026/9/15 0:28:03 · 阅读 1.2k 上海土壤水文分组高精度栅格制作与SCS-CN径流模拟应用
简介上海市土壤水文分组HSG高精度栅格数据是基于USDA曲线数CN方法估算降雨径流的关键基础数据适用于SWAT水文模型建模、城市排水规划及流域水文分析。数据源自HYSOGs250m与soilGrids250m系统按省级行政区裁剪为上海市范围WGS84地理坐标系空间分辨率约250m提供A、B、C、D四类径流潜力等级并包含湿润土壤双重HSG分类信息。压缩包大小3.46MB文件总数0类型明细暂无数据可用GIS软件直接加载分析。目前已有212人学习适合水文、环境、GIS领域的研究者与工程师可直接用于计算径流曲线数、划分水文响应单元及开展区域水资源评估。1. 上海市土壤水文分组高精度栅格给产流模型一张能直接查表的底图做上海市域内涝模拟或海绵城市绩效评估的人迟早会被同一个问题卡住SCS-CN 查表需要土壤水文分组HSGA/B/C/D 四档而手头常见的是二普扫描图、HWSD 的 1 km 网格或 SoilGrids 的 250 m 插值。上海市土壤水文分组高精度栅格就是把上海地表土壤按入渗能力归入四档做成覆盖全市的 10 m 或 30 m 整幅栅格让 SWMM、InfoWorks ICM、HEC-HMS 在每个像元上直接完成 HSG 到 CN 的映射。难在上海几乎每一项地表条件都在抬升难度高差几米、地下水位贴着地表、硬化面把土体切得支离破碎而 1:50 万土壤图投影到 10 m 栅格本身就是一种假精细。下面按判定、制作、验证、对接四步把链路拆开。2. A 到 D 四档怎么定USDA 水文土壤分组的判定参数与上海土壤归属逻辑2.1 水文土壤分组不是土壤分类而是入渗能力分级中国土壤发生分类里的水稻土、潮土、滨海盐土回答的是土壤怎么来的USDA 水文土壤分组回答的是另一件事长时间湿润、无冻土、无地表封盖条件下土壤在强降雨里能吞下多少水。NRCS 在 National Engineering Handbook Part 630 中按入渗能力把土壤分成 A、B、C、D 四组A 组通常为砂质或砾质入渗极快B 组中等C 组偏慢D 组极慢。D 组有一条补充定义对上海特别致命即使质地不粘只要下挖 60 cm 内遇到持续地下水位或 50 cm 内遇到致密不透水层也判为 D。上海大部分地区地下水位在 0.51.5 m雨季贴着地表这一条直接锁定了大片 D 组。2.2 先看质地再看 Ks最后用排水条件兜底判定流程我习惯按三步走。第一步查土种的质地剖面按 USDA 质地三角图映射出候选组第二步用饱和导水率 Ks 复核实验室定水头法、现场双环入渗仪都可以第三步用地下水埋深和不透水层深度做一票否决。三者对应关系如下HSG典型质地稳态入渗率 / Ks (mm/h)排水与水位条件A砂土、壤质砂土大于 7.6排水迅速无滞水迹象B砂质壤土、壤土、粉砂壤土3.8 ~ 7.6中等排水C粘壤土、砂质粘壤土、粉砂质粘壤土1.3 ~ 3.8有弱滞水层排水较慢D粘土、重粘土或任意质地但 60 cm 内有水位小于 1.3排水差或地下水位高表里的 Ks 取的是 NRCS 常用参考区间。实际判定时看的是整层剖面中最差的那个层次不是表层——上海水稻土的耕作层透水性尚可真正决定分组的是犁底层。2.3 上海主要土种的归属判断按土类分布看南汇、奉贤沿海的滨海盐土粉砂质粘壤土加常年高水位归 D青浦、松江、金山的青紫泥、青黄土等水稻土犁底层紧实、下有潜育层归 D嘉定、宝山及黄浦江以西的灰潮土粉砂壤土到粘壤土、水位 12 m多为 C崇明东滩、长兴岛新近沉积的砂质潮土局部到 B砂性强的可到 A。老城区内的回填土和建筑渣土扰动层属性无法从二普图获得单独编码为 9 类扰动土不参与常规查表。具体归属如下土类典型分布控制性剖面特征HSG滨海盐土南汇、奉贤沿海粉砂质粘壤土水位常年小于 1 mD青紫泥 / 青黄土水稻土青浦、松江、金山粘壤土加紧实犁底层加潜育层D灰潮土嘉定、宝山、黄浦江以西粉砂壤土至粘壤土水位 12 mC砂质潮土崇明、长兴岛、古河道带砂质壤土剖面疏松A / B2.4 从粒径或质地名称自动落组字典映射就够如果数据源给的是粘粒、砂粒、粉粒百分数先按 USDA 三角图定质地类再按字典映射落到四档TEXTURE_HSG { sand: A, loamy sand: A, sandy loam: B, loam: B, silt loam: B, silt: B, sandy clay loam: C, clay loam: C, silty clay loam: C, sandy clay: D, silty clay: D, clay: D, } def texture_to_hsg(texture: str) - str: return TEXTURE_HSG.get(texture.strip().lower(), None)映射本身很朴素功夫在前处理把《上海土壤》里每个土属的典型质地剖面整理成 CSV逐层读入取剖面中最差一层的质地做最终归属。老资料里中壤重壤这类卡庆斯基制名称要先转换成 USDA 制的砂粒、粉粒、粘粒百分数再进三角图两套制式的质地边界并不重合。粒径数据充足时用 Saxton-Rawls pedotransfer 函数把粒径、容重、有机质换算成 Ks再按 2.2 节的阈值分组比纯质地映射更贴实测但该函数标定数据主要来自美国土壤对上海高粉粒、高盐基的冲积土外推时建议留 20% 左右的余量。注意城市绿地里的客土绿化工程换填的种植土资料上常标为砂壤土实际压实后入渗率往往跌到 C 组遇到这种情况以现场探测为准。3. 高精度栅格怎么造数据源分工、投影统一与 rasterio 栅格化参数3.1 先搞清楚高精度到底来自哪里1:50 万上海市第二次土壤普查数字化图的精度在制图综合层面最小图斑对应的地面范围远大于 10 m 像元HWSD 是 30 弧秒网格SoilGrids 是 250 m。直接拿任何一个单一来源做 10 m 输出都只是把粗边界切细叫高精度名不副实。常见做法是把来源按职责拆开让各自干各自擅长的事数据来源尺度在高精度栅格里的职责1:50 万上海土壤图二普数字化版上海农科院、土地档案整理图斑级主边界与土种属性《上海土壤》土种志1990 年代公开出版属性表质地、剖面、水位逐层记录ALOS PALSAR 12.5 m DEMASF / OpenTopography 分发12.5 m河漫滩、古河道、湖沼平原地貌细分上海土地利用 10 m 产品测绘遥感解译10 m城市扰动区标记、水体掩膜SoilGrids 250 mISRIC250 m粒径空间插值的辅助背景这样组合出来的栅格边界精度来自土壤图内部细分来自地形与用地信息才算真正对得起高精度三个字。3.2 投影与像元对齐统一到本地投影坐标系上海市常用 CGCS2000 基准、3° 带高斯-克吕格投影中央经线取 121.5°E单位为米。生产上第一步是把所有矢量重投影到与模板栅格一致的坐标系模板一般取土地利用栅格或 DEM。像元对齐有三条硬约束像元尺寸取 10 m 的整数倍像元原点左上角与模板完全一致行列数不能四舍五入否则整幅栅格会沿东北方向漂移半个像元与路网叠合时错位肉眼可见。最省事的做法是让模板栅格先定义好 transform后续所有栅格化和重采样都以它为唯一基准。3.3 rasterize 参数与输出写法import geopandas as gpd import numpy as np import rasterio from rasterio.features import rasterize # hsg_code: 1A, 2B, 3C, 4D, 9扰动土, -9999无数据 gdf gpd.read_file(sh_soil_50w_hsg.shp) with rasterio.open(sh_lu_10m_template.tif) as ref: gdf gdf.to_crs(ref.crs) # 统一坐标系 transform ref.transform out_shape (ref.height, ref.width) shapes [(geom, int(code)) for geom, code in zip(gdf.geometry, gdf[hsg_code])] raster rasterize( shapes, out_shapeout_shape, transformtransform, fill-9999, dtypenp.int16, all_touchedFalse, ) profile ref.profile.copy() profile.update(dtypenp.int16, count1, nodata-9999, compressdeflate) with rasterio.open(sh_hsg_10m_raw.tif, w, **profile) as dst: dst.write(raster, 1)rasterize 的三个参数值得单独说。all_touchedFalse 表示只有像元中心落在图斑内才赋值避免沿边界多占一圈像元对 10 m 栅格和细碎图斑先设 False若边界处出现大量空洞再改 True。dtype 用 int16 不是因为分组只有 4 个值而是要为 9扰动土和 -9999无数据留编码空间同时保持与土地利用栅格一致的整型语义后续 CN 查表时布尔索引的速度也更接近 C 级。fill-9999 先铺底栅格化后未覆盖区域全部是 -9999后面统一用上海市域掩膜裁掉不加掩膜的话沿海潮间带和江面里的沉积物会被当成真实土壤算进 D 组。3.4 边界修正与孔洞清理gdal_sieve 的用法原始栅格化结果有两个典型毛病河岸边土壤图与水体相交产生锯齿状伪像元细碎图斑在 all_touchedFalse 下变成孔洞。处理分两步。先做水体掩膜把黄浦江、长江口、淀山湖等水体多边形栅格化为 -9999再做连通域清理用 GDAL 自带的 sieve 工具去掉小于阈值的碎斑gdal_sieve.py -st 20 -8 sh_hsg_10m_raw.tif sh_hsg_10m_sieve.tif-st 20 表示小于 20 个像元的连通域并入周边最大类-8 表示八邻域连通判定。阈值按碎斑分布调上海市区地块被道路切得很碎20 太小会保留地块内部的图斑噪声我一般给到 3050。最后把城市建设用地边界内的像元用土地利用图覆写为 9 类扰动土避免老图斑在建成区复活。sieve 只改空间连续性不改变大类内部的真实过渡对 10 m 栅格是安全的。提示如果下游是 SWMM 这类子汇水区模型栅格分辨率取子汇水区平均尺寸的 1/4 到 1/2 即可过细只增加重采样耗时不会提高 CN 精度。4. 上海土壤水文分组栅格的精度验证样点布局、混淆矩阵与边界效应排查4.1 样点布局分层抽样比均匀撒点可靠验证样本按 HSG 和地貌单元双重分层抽取。每个 HSG 组不少于 20 个点全市 80120 个点足以支撑四类混淆矩阵。按二项分布估算期望总体精度 80%、允许误差 10%、置信水平 95% 时最少约 62 个点和这个范围吻合。点位约束有三条距图斑边界至少一个像元避开道路硬化面避开已竣工基坑。现场每点挖 60 cm 剖面读质地与层次再用双环入渗仪测稳定入渗率剖面判分组入渗率用于解释偏差来源而不是直接改判——单点入渗受压实、根系和裂隙影响很大一次读数不能代表整个图斑。4.2 混淆矩阵、总体精度与 Kappa把实测编码与栅格提取编码对齐后用 sklearn 一次算出三个指标import numpy as np from sklearn.metrics import confusion_matrix, cohen_kappa_score obs np.array([1, 1, 2, 2, 3, 3, 3, 4, 4, 4]) # 实测 HSG pred np.array([1, 2, 2, 2, 3, 4, 3, 4, 4, 4]) # 栅格提取 cm confusion_matrix(obs, pred, labels[1, 2, 3, 4]) oa np.trace(cm) / cm.sum() kappa cohen_kappa_score(obs, pred) producer np.diag(cm) / cm.sum(axis0) # 生产者精度 user np.diag(cm) / cm.sum(axis1) # 用户精度一组合格的 90 点验证结果如下预测 A预测 B预测 C预测 D合计用户精度实测 A162101984.2%实测 B318212475.0%实测 C022022483.3%实测 D001222395.7%合计1922242590生产者精度84.2%81.8%83.3%88.0%总体精度 84.4%Kappa 0.79。这个数放在 1:50 万源图衍生的 10 m 产品里是合格的真正的诊断信息在混淆结构上B 与 C 互相错分占了大头说明质地过渡带被硬切成两档的边界误差远大于整块分错类的误差。Kappa 低于 0.8 时不要急着改栅格先回到 2.3 节的剖面最差层判定逻辑复查原始属性表。4.3 三类最容易出错的边界河岸带是重灾区。黄浦江和吴淞江沿岸的砂质透镜体多是古河道摆动留下的土壤图上没有要用 DEM 的河漫滩相对高度和影像纹理补。做法是取距水系 100 m 缓冲带叠加相对高程小于 2 m 的区域把原判 C 的像元复核为 B。第二类是田埂和大棚菜地犁底层和压实层让剖面整体下移半档这类地类面积逐年缩小但残余图斑对 CN 的影响不小。第三类是建设用地边缘回填土与原地表犬牙交错栅格上表现为 9 类像元包围 D 组的椒盐噪声。直接做 majority filter 会抹掉真实信息我一般保留原值在属性里加 flag 字段模型侧对 9 类单独设 CN而不是硬并进 C 或 D。5. 把 HSG 栅格接进上海 SCS-CN 径流模拟CN 查表合成、AMC 修正与反推验证5.1 TR-55 的 HSG-CN 查表关系HSG 单独不产生水文意义意义在 SCS-CN 的查表环节。TR-55 给每个土地利用类别配了四组 CN 默认值上海的常见类别如下土地利用上海常见类别ABCD绿地 / 公园良好状态39617480低密度住宅1/3 英亩61758387高密度住宅1/8 英亩77859092商业 / 工业85% 不透水89929495旱作农田常规耕作67758387硬化路面 / 停车场98989898这张表的工程潜台词是土地利用对 CN 的影响远大于 HSG 分组。上海市区以高密度建成区为主HSG 的误差只有在绿地、农田这类透水地类上才会放大所以验证资源应当优先投到郊野和蓝绿空间。5.2 栅格级 CN 合成的 numpy 写法把 HSG 栅格与同期土地利用栅格对齐到同一像元网格后按查表逐项赋值import numpy as np import rasterio CN_TABLE { (1, 1): 39, (2, 1): 61, (3, 1): 74, (4, 1): 80, # 绿地 (1, 2): 89, (2, 2): 92, (3, 2): 94, (4, 2): 95, # 商服/工业 (1, 3): 98, (2, 3): 98, (3, 3): 98, (4, 3): 98, # 不透水 (1, 4): 67, (2, 4): 75, (3, 4): 83, (4, 4): 87, # 旱作农田 } with rasterio.open(sh_hsg_10m.tif) as hs, \ rasterio.open(sh_lu_10m.tif) as lu: hsg hs.read(1) land lu.read(1) cn np.full(hsg.shape, -9999, dtypenp.float32) for (h, l), val in CN_TABLE.items(): cn[(hsg h) (land l)] val profile hs.profile.copy() profile.update(dtypenp.float32, nodata-9999, compressdeflate) with rasterio.open(sh_cn_10m.tif, w, **profile) as dst: dst.write(cn, 1)这里有一个覆盖顺序的坑字典遍历的先后决定同名像元谁胜出。上面这张表各键互斥所以安全若引入水域、湿地等多条件叠加的类别必须把优先级最高的类别放最后赋值或用 np.select 按条件数组顺序求值。5.3 三个必调参数与一个反推验证技巧第一个必调参数是前期土壤湿度条件 AMC。上海梅雨和台风季前期土壤接近饱和AMC-III 的 CN 比 AMC-II 高 515只出一版栅格会系统性低估产流。常见做法是准备两版 CN 栅格用前 5 日降雨量判断当次取值AMC-II 到 AMC-III 的换算用 SWAT 近似式 CN3 CN2 * exp(0.00673 * (100 - CN2))。第二个是初损率 λ。TR-55 默认 0.2上海城镇流域的短历时降雨用 0.05 与实测过程线拟合更好S 的换算按 Hawkins 的 S05 1.33 * S02^1.15英寸制重算不能直接套 0.2 的 S 值。第三个是聚合重采样。10 m 栅格聚合到子汇水区时只能取众数或面积加权平均双线性插值会把四类整数插出 3.7 这样的伪类。注意凡 CN_TABLE 未覆盖的像元nodata、9 类扰动土最后都会落在 -9999 上模型侧要按无效值剔除不要用 0 兜底0 会被当成 CN100 参与计算。验证技巧挑上海近年来 35 场有完整雨量与出口流量记录的短历时降雨反推流域平均 CN再与同流域栅格 CN 的面积加权均值对比。偏差在 ±5 以内说明 HSG 栅格与土地利用编码的配合是自洽的若偏差方向一致且普遍偏大优先检查土地利用的 CN 取值和不透水率假设而不是回头改土壤分组。本文还有配套的精品资源点击获取