基于BASINS来创建HSPF模型

• basin的工程文件包可以自定义到任意指定位置
• 创建的hspf的控制文件固定保存在basins安装目录的modelout文件夹下面
• 运行hspf后，计算结果也在modelout文件夹下相应的hspf工程文件夹
• 输出文件夹里的.ech文件记录了hspf模型运行的过程，如发散，可在此文件里寻求解决方案

TauDEM 命令运行

• 8为计算核数
• logan为原始dem文件名，根据实际案例替换
• 运行完成后将生成的子流域栅格数据，在QGIS中转换为矢量数据
• 然后利用dem、河网、子流域数据，在Basins里选择手动流域描绘，计算流域特征参数，作为HSPF的输入数据
• dem文件不能大于4G，尽量小一点，可以将生成的矢量数据再进行拼接

mpiexec -n 8 PitRemove logan.tif
mpiexec -n 8 D8Flowdir -p loganp.tif -sd8 logansd8.tif -fel loganfel.tif
mpiexec -n 8 DinfFlowdir -ang loganang.tif -slp loganslp.tif -fel loganfel.tif
mpiexec -n 8 AreaD8 -p loganp.tif -ad8 loganad8.tif
mpiexec -n 8 AreaDinf -ang loganang.tif -sca logansca.tif
mpiexec -n 8 Aread8 -p loganp.tif -o loganoutlet.shp -ad8 loganad8o.tif
mpiexec -n 8 Gridnet -p loganp.tif -plen loganplen.tif -tlen logantlen.tif -gord logangord.tif 
mpiexec -n 8 PeukerDouglas -fel loganfel.tif -ss loganss.tif
mpiexec -n 8 Aread8 -p loganp.tif -o loganoutlet.shp -ad8 loganssa.tif -wg loganss.tif
mpiexec -n 8 Dropanalysis -p loganp.tif -fel loganfel.tif -ad8 loganad8.tif -ssa loganssa.tif -drp logandrp.txt -o loganoutlet.shp -par 5 500 10 0
mpiexec -n 8 Threshold -ssa loganssa.tif -src logansrc.tif -thresh 300
mpiexec -n 8 Streamnet -fel loganfel.tif -p loganp.tif -ad8 loganad8.tif -src logansrc.tif -ord loganord3.tif -tree logantree.dat -coord logancoord.dat -net logannet.shp -w loganw.tif -o loganoutlet.shp
mpiexec -n 8 WatershedGridToShapefile loganw.tif loganw.shp

大流域情况下，使用QGIS手动建立流域拓扑关系

• 子流域图层
  • 添加SUBBASIN字段
  • 添加SLO1字段
  • 添加AREAACRES字段，=几何面积（m2）/4046.86
  • 添加AREAMI2字段，=几何面积（m2）/2589988
• 河段图层
  • 添加LEVEL字段，=1
  • 根据子流域边界裁剪河网，每个子流域内只保留最下游的河段
  • 添加子流域下游子流域SUBBASINR字段
  • 添加河段长度LEN2字段，m
  • 添加河段local汇水区面积LAREA字段，m2，=所在子流域面积
  • 添加河段总汇水区面积TAREA字段，m2，=河段所有上游子流域面积和
  • 添加河段总汇水面积TAREAACRES字段，acres，=TAREA/4046.86
  • 添加河段总汇水面积TAREAMI2字段，square miles，=TAREA/2589988
  • 添加河段宽度WID2字段，m，=(1.29) * ((r / 1000000) ^ (0.6)) 'r is Area in Sq. m. r/1000000 is the area in Sq. Km. Area为河段总汇水区面积。
  • 添加河段深度DEP2字段，m，=(0.13) * ((r / 1000000) ^ (0.4)) 'r is Area in Sq. m. r/1000000 is the area in Sq. Km. Area为河段总汇水区面积。
  • 添加河段起点高程MINEL字段，m，=起点所在DEM值
  • 添加河段终点高程MAXEL字段，m，=终点所在DEM值
  • 添加河段坡度SLO2字段，百分比，=（MAXEL-MINEL）*100/LEN2
  • 添加河段名称SNAME字段，=自定义字符串
  • 在子流域图层添加BNAME字段，=对应河段的SNAME
• 流域出口图层
  • 每个河段创建一个流域出口
  • 添加ID字段
  • 添加PCSID字段
  • 添加Xpr字段，河段终点x坐标
  • 添加Ypr字段，河段终点y坐标
• 计算各参数值
  • 子流域属性
    • 面积：直接在field calculator，$area，再换算单位
    • SLO1：待定
  • 河段图层属性：
    • SUBBASIN：与子流域图层相交分析，join属性
    • SUBBASINR：结合拓扑赋值
    • LEN2：直接在field calculator，$length
    • LAREA：与子流域图层相交分析，join属性
    • TAREA：根据SUBBASIN、SUBBASINR字段在field calculator里写递归算法计算
    • TAREAACRES、TAREAMI2、WID2、DEP2根据TAREA计算
    • MINEL、MAXEL: 首先利用Extract specific vertices工具生成河段端点（0，-1），然后利用Point Sampling Tool插件结合DEM提取端点高程，然后在field calculator里编写处理代码
    • SLO2：由MINEL、MAXEL字段计算

递归计算流域面积, 核查如果计算结果有误，可以将表单拷出为csv文件，直接写python程序来处理，处理后用编号关联属性

from qgis.core import *
from qgis.gui import *

@qgsfunction(args='auto', group='Custom', referenced_columns=[])
def cal_upstream_area(value1,value2, feature, parent):
    """
    value1： SUBBASIN 字段值
    value2： LAREA 字段值
    计算流域出口总汇水面积
    """
    res=0
    layer = qgis.utils.iface.activeLayer()
    prov = layer.dataProvider()
    field_names = list(layer.attributeAliases().keys())
    f1 = layer.fields().indexFromName('SUBBASIN')
    f2 = layer.fields().indexFromName('SUBBASINR')
    f3 = layer.fields().indexFromName('LAREA')
    features = layer.getFeatures()
    sbnm = value1
    area = value2
    area = area+calculate_sum(sbnm,features,f1,f2,f3)
    return area

def calculate_sum(code, data,f1,f2,f3):
    # 找到所有与给定编码关联的项
    related_items = [x for x in data if x.attributes()[f2] == code]
    
    if not related_items:
        return 0  # 找不到项时返回0
    
    # 计算所有上游项的第三列和
    upstream_sum = 0
    for item in related_items:
        upstream_sum += item.attributes()[f3]
        upstream_sum += calculate_sum(item.attributes()[f1], data,f1,f2,f3)

    return upstream_sum

河段端点高程处理

from qgis.core import *
from qgis.gui import *

@qgsfunction(args='auto', group='Custom', referenced_columns=[])
def seq_up_dwn(value1, value2, feature, parent):
    """
    将高程高的作为线段起点
    """
    res=0
    layer = qgis.utils.iface.activeLayer()
    prov = layer.dataProvider()
    field_names = list(layer.attributeAliases().keys())
    f1 = layer.fields().indexFromName('SUBBASIN')
    f2 = layer.fields().indexFromName('jsj_yn_bas')
    features = layer.getFeatures()
    attrs = feature.attributes()
    current_id = feature.id()
    if current_id<layer.featureCount():
        for next_f in features:
            next_f_attrs = next_f.attributes()
            if next_f_attrs[f1] == attrs[f1]:                
                h1 = attrs[f2]
                h2 = next_f_attrs[f2]
                h=[]
                h.append(h1)
                h.append(h2)
                res= max(h)
                break
    return res

存储到csv文件处理

import numpy as np

def calculate_sum(code, data,f1,f2,f3):
    # 找到所有与给定编码关联的项
    related_items = [x for x in data if int(x[f2]) == code]
    
    if not related_items:
        return 0  # 找不到项时返回0
    
    # 计算所有上游项的第三列和
    upstream_sum = 0
    for item in related_items:
        upstream_sum += item[f3]
        upstream_sum += calculate_sum(item[f1], data,f1,f2,f3)

    return upstream_sum

data = np.genfromtxt(r'input.csv', delimiter=',')
od = []
for d in data:
    area = [x for x in data if int(x[0]) == int(d[0])][0][2]
    area =area+ calculate_sum(int(d[0]),data,0,1,2)
    od.append([int(d[0]),area])
np.savetxt(r'output.csv', od, delimiter=',')