使用Python编写集中式新安江代码（二、三水源）

首先设置模型参数

import numpy as np
import pandas as pd
import math
df = pd.read_excel(r"C:\Users\Admin\Desktop\pet.xlsx")
df = df.fillna(0)
print(df)
WM,WUM,WLM,WDM,b,K,C,FC,F,T,CG,L,CR = 120,15,85,20,0.3,0.95,0.14,15,2856,24,0.9,5,0.5    # F单位为km^2 T单位为h
Sm,KI,KG,EX,FR1,CI,Q = 20,0.2,0.5,1.5,0.4175,0.9,5.6
WMM = (1+b)*WM
U = F / (3.6 * T)
SMM = Sm * (1+EX)

我的数据是excel数据其中有流域面平均雨量和蒸发资料

首先是设置一个for循坏我一共5000多行降雨数据，设置一个for循坏

for i in range(0,5479):
# 计算土湿
    if i == 0:
        df.loc[i, "WU"] = 0
        df.loc[i,"WL"] = 2.2
        df.loc[i, "WD"] = 20
        df.loc[i, "EP"] = K * df.loc[i, "E0"]
        if df.loc[i, "WU"] + df.loc[i, "P"] >= df.loc[i, "EP"]:
            df.loc[i, "EU"] = df.loc[i, "EP"]
            df.loc[i, "EL"] = 0
            df.loc[i, "ED"] = 0
        if df.loc[i, "WU"] + df.loc[i, "P"] < df.loc[i, "EP"] and df.loc[i, "WL"] >= C * WLM:
            df.loc[i, "EU"] = df.loc[i, "WU"] + df.loc[i, "P"]
            df.loc[i, "EL"] = (df.loc[i, "EP"] - df.loc[i, "EU"]) * df.loc[i, "WL"] / WLM
            df.loc[i, "ED"] = 0
        if df.loc[i, "WU"] + df.loc[i, "P"] < df.loc[i, "EP"] and df.loc[i, "WL"] < C * WLM and df.loc[i, "WL"] >= C * (df.loc[i, "EP"] - df.loc[i, "EU"]):
            df.loc[i, "EU"] = df.loc[i, "WU"] + df.loc[i, "P"]
            df.loc[i, "EL"] = C * (df.loc[i, "EP"] - df.loc[i, "EU"])
            df.loc[i, "ED"] = 0
        if df.loc[i, "WU"] + df.loc[i, "P"] < df.loc[i, "EP"] and df.loc[i, "WL"] < C * (
                df.loc[i, "EP"] - df.loc[i, "EU"]):
            df.loc[i, "EU"] = df.loc[i, "WU"] + df.loc[i, "P"]
            df.loc[i, "EL"] = df.loc[i, "WL"]
            df.loc[i, "ED"] = C * (df.loc[i, "EP"] - df.loc[i, "EU"]) - df.loc[i, "EL"]
        df.loc[i, "E"] = df.loc[i, "EU"] + df.loc[i, "EL"] + df.loc[i, "ED"]
        df.loc[i, "PE"] = df.loc[i, "P"] - df.loc[i, "E"]
        df.loc[i, "W"] = df.loc[i, "WU"] + df.loc[i, "WL"] + df.loc[i, "WD"]
    else:
        dw = df.loc[i-1,"PE"] - df.loc[i-1,"R"] #这边需要连接下面R，否则会报错
        df.loc[i,"WU"] = df.loc[i-1,"WU"] + dw
        df.loc[i,"WL"] = df.loc[i-1,"WL"]
        df.loc[i,"WD"] = df.loc[i-1,"WD"]
        if df.loc[i,"WU"] < 0 :
            df.loc[i,"WL"] = df.loc[i-1,"WL"] + dw + df.loc[i-1,"WU"]
            df.loc[i,"WD"] = df.loc[i-1,"WD"]
            df.loc[i,"WU"] = 0
            if df.loc[i,"WL"] < 0:
                df.loc[i,"WD"] = df.loc[i-1,"WD"] + dw + df.loc[i-1,"WU"] + df.loc[i-1,"WL"]
                df.loc[i,"WU"] = 0
                df.loc[i,"WL"] = 0
                if df.loc[i,"WD"] < 0:
                    df.loc[i,"WU"] = 0
                    df.loc[i,"WD"] = 0
                    df.loc[i,"WL"] = 0
        if df.loc[i,"WU"] > WUM:
            df.loc[i,"WU"] = WUM
            df.loc[i,"WL"] = df.loc[i-1,"WL"] + dw - WUM + df.loc[i-1,"WU"]
            df.loc[i,"WD"] = df.loc[i-1,"WD"]
            if df.loc[i,"WL"] > WLM:
                df.loc[i,"WU"] = WUM
                df.loc[i,"WL"] = WLM
                df.loc[i,"WD"] = df.loc[i-1,"WD"] + dw - WLM - WUM + df.loc[i-1,"WU"] + df.loc[i-1,"WL"]
        df.loc[i,"EP"] = K * df.loc[i,"E0"]
        if df.loc[i,"WU"] + df.loc[i,"P"] >= df.loc[i,"EP"]:
            df.loc[i,"EU"] = df.loc[i,"EP"]
            df.loc[i,"EL"] = 0
            df.loc[i,"ED"] = 0
        if df.loc[i,"WU"] + df.loc[i,"P"] < df.loc[i,"EP"] and df.loc[i,"WL"] >= C * WLM:
            df.loc[i,"EU"] = df.loc[i,"WU"] + df.loc[i,"P"]
            df.loc[i,"EL"] = (df.loc[i,"EP"] - df.loc[i,"EU"]) * df.loc[i,"WL"] / WLM
            df.loc[i,"ED"] = 0
        if df.loc[i,"WU"] + df.loc[i,"P"] < df.loc[i,"EP"] and df.loc[i,"WL"] < C * WLM and df.loc[i,"WL"] >= C * (df.loc[i,"EP"]-df.loc[i,"EU"]):
            df.loc[i,"EU"] = df.loc[i,"WU"] + df.loc[i,"P"]
            df.loc[i,"EL"] = C * (df.loc[i,"EP"] - df.loc[i,"EU"])
            df.loc[i,"ED"] = 0
        if df.loc[i,"WU"] + df.loc[i,"P"] < df.loc[i,"EP"] and df.loc[i,"WL"] < C * (df.loc[i,"EP"]-df.loc[i,"EU"]):
            df.loc[i,"EU"] = df.loc[i,"WU"] + df.loc[i,"P"]
            df.loc[i,"EL"] = df.loc[i,"WL"]
            df.loc[i,"ED"] = C * (df.loc[i,"EP"] - df.loc[i,"EU"]) - df.loc[i,"EL"]
        df.loc[i,"E"] = df.loc[i,"EU"] + df.loc[i,"EL"] + df.loc[i,"ED"]
        df.loc[i,"PE"] = df.loc[i,"P"] - df.loc[i,"E"]
        df.loc[i,"W"] = df.loc[i,"WU"] + df.loc[i,"WL"] + df.loc[i,"WD"]

再根据自己的要求计算R

# 计算R
    a = WMM * (1- math.pow((1-df.loc[i,"W"]/WM),1/(1+b)))
    if i ==0:
        if df.loc[i,"PE"] > 0:
            if a + df.loc[i,"PE"] <= WMM:
                df.loc[i,"R"] = df.loc[i,"PE"] + df.loc[i,"W"] - WM + WM * math.pow(1-((df.loc[i,"PE"]+a)/WMM),b+1)  # df["R"][i] = df["PE"][i] + df["W"][i] - WM + WM * math.pow(1 - ((df["PE"][i] + a) / WMM), b + 1)
            if a + df.loc[i,"PE"] > WMM:
                df.loc[i,"R"] = df.loc[i,"PE"] - (WM - df.loc[i,"W"])
        if df.loc[i,"R"] < 0:
            df.loc[i,"R"] = 0
    else:
        if df.loc[i,"PE"] > 0:
            if a + df.loc[i, "PE"] <= WMM:
                 df.loc[i, "R"] = df.loc[i, "PE"] + df.loc[i, "W"] - WM + WM * math.pow(1 - ((df.loc[i, "PE"] + a) / WMM),b + 1)  # df["R"][i] = df["PE"][i] + df["W"][i] - WM + WM * math.pow(1 - ((df["PE"][i] + a) / WMM), b + 1)
            if a + df.loc[i, "PE"] > WMM:
                df.loc[i, "R"] = df.loc[i, "PE"] - (WM - df.loc[i, "W"])
        if df.loc[i, "R"] < 0:
            df.loc[i, "R"] = 0

在计算产流面积

# 计算流域产流面积
    if i == 0:
        if df.loc[i,"R"] > 0 :
            df.loc[i,"FR"] = df.loc[i,"R"] / df.loc[i,"PE"]
            if df.loc[i,"FR"] >1:
                df.loc[i,"FR"] =1

        else:
            df.loc[i,"FR"] = FR1
    else:
        if df.loc[i, "R"] > 0:
            df.loc[i, "FR"] = df.loc[i, "R"] / df.loc[i, "PE"]
            if df.loc[i, "FR"] > 1:
                df.loc[i, "FR"] = 1
        else:
            df.loc[i, "FR"] = df.loc[i-1,"FR"]

这边是二水源代码

# 划分二水源
#     if i == 0:
#         if df.loc[i, "R"] > 0:
#             df.loc[i, "FR"] = df.loc[i, "R"] / df.loc[i, "PE"]
#             if df.loc[i,"FR"] > 1:
#                 df.loc[i,"FR"] = 1
#             if df.loc[i, "FR"] > 0:
#                 if df.loc[i, "PE"] > FC:
#                     df.loc[i, "RS"] = (df.loc[i, "PE"] - FC) * df.loc[i, "FR"]
#                     df.loc[i, "RG"] = df.loc[i, "R"] - df.loc[i, "RS"]
#                 if df.loc[i, "PE"] <= FC:
#                     df.loc[i, "RS"] = 0
#                     df.loc[i, "RG"] = df.loc[i, "R"]
#             if df.loc[i, "FR"] < 0:
#                 df.loc[i, "FR"] = 0
#                 df.loc[i, "RS"] = 0
#                 df.loc[i, "RG"] = 0
#         if df.loc[i, "R"] == 0:
#             df.loc[i, "RS"] = 0
#             df.loc[i, "RG"] = 0
#     else:
#         if df.loc[i,"R"] > 0 :
#             df.loc[i,"FR"] = df.loc[i,"R"] / df.loc[i,"PE"]
#             if df.loc[i,"FR"] > 0:
#                 if df.loc[i,"PE"] > FC:
#                     df.loc[i,"RS"] = (df.loc[i,"PE"] - FC) * df.loc[i,"FR"]
#                     df.loc[i,"RG"] = df.loc[i,"R"] - df.loc[i,"RS"]
#                 if df.loc[i,"PE"] <= FC:
#                     df.loc[i,"RS"] = 0
#                     df.loc[i,"RG"] = df.loc[i,"R"]
#             if df.loc[i,"FR"] < 0:
#                 df.loc[i,"FR"] = 0
#                 df.loc[i,"RS"] = 0
#                 df.loc[i,"RG"] = 0
#         if df.loc[i,"R"] == 0:
#             df.loc[i,"RS"] = 0
#             df.loc[i,"RG"] = 0

三水源

# 划分三水源
    if i == 0:     # 第一时段计算
        df.loc[i, "S1"] = 3.4525
        if df.loc[i,"PE"] > 0:
            AU = SMM * (1 - math.pow((1-(((df.loc[i,"S1"]*FR1)/df.loc[i,"FR"])/Sm)),1/(1+EX)))
            if df.loc[i,"PE"] + AU < SMM:
                df.loc[i,"RS"] = df.loc[i,"FR"] * (df.loc[i,"PE"] +(df.loc[i,"S1"] * FR1)/df.loc[i,"FR"] - Sm +
                                                 Sm *math.pow((1-(df.loc[i,"PE"] +AU)/SMM),1+EX))
            if df.loc[i,"PE"] + AU >= SMM:
                df.loc[i,"RS"] = df.loc[i,"FR"] * (df.loc[i,"PE"] + (df.loc[i,"S1"] * FR1)/df.loc[i,"FR"] - Sm)
            S = (df.loc[i,"S1"] * FR1)/df.loc[i,"FR"] + (df.loc[i,"R"] - df.loc[i,"RS"])/df.loc[i,"FR"]
            df.loc[i,"RI"] = KI * S *df.loc[i,"FR"]
            df.loc[i,"RG"] = KG * S *df.loc[i,"FR"]
            df.loc[i+1,"S1"] = S*(1-KI-KG)
        else:
            S = (df.loc[i, "S1"] * FR1) / df.loc[i, "FR"]
            df.loc[i + 1, "S1"] = S * (1 - KG - KI)
            df.loc[i, "RS"] = 0
            df.loc[i, "RG"] = KI*S * df.loc[i,"FR"]
            df.loc[i, "RI"] = KG * S *df.loc[i,"FR"]
    else: #其余时段计算
        if df.loc[i, "PE"] > 0:
            AU = SMM * (1 - math.pow((1-(((df.loc[i,"S1"]*df.loc[i-1,"FR"])/df.loc[i,"FR"])/Sm)),1/(1+EX)))
            if df.loc[i,"PE"] + AU < SMM:
                df.loc[i,"RS"] = df.loc[i,"FR"] * (df.loc[i,"PE"] +(df.loc[i,"S1"] * df.loc[i-1,"FR"])/df.loc[i,"FR"] - Sm +
                                                 Sm *math.pow((1-(df.loc[i,"PE"] +AU)/SMM),1+EX))
            if df.loc[i,"PE"] + AU >= SMM:
                df.loc[i,"RS"] = df.loc[i,"FR"] * (df.loc[i,"PE"] + (df.loc[i,"S1"] * df.loc[i-1,"FR"])/df.loc[i,"FR"] - Sm)
            S = (df.loc[i,"S1"] * df.loc[i-1,"FR"])/df.loc[i,"FR"] + (df.loc[i,"R"] - df.loc[i,"RS"])/df.loc[i,"FR"]
            df.loc[i,"RI"] = KI * S *df.loc[i,"FR"]
            df.loc[i,"RG"] = KG * S *df.loc[i,"FR"]
            df.loc[i+1,"S1"] = S*(1-KI-KG)
        else:
            S = (df.loc[i, "S1"] * df.loc[i-1,"FR"]) / df.loc[i, "FR"]
            df.loc[i + 1, "S1"] =S * (1 - KG - KI)
            df.loc[i, "RS"] = 0
            df.loc[i, "RG"] = KG *S *df.loc[i,"FR"]
            df.loc[i, "RI"] = KI * S *df.loc[i,"FR"]

二水源汇流

# # 计算二水源汇流
#     if i == 0:
#         df.loc[i, "Qs"] = df.loc[i, "RS"] * U
#         df.loc[i,"Qg"] = 5.6
#         df.loc[i, "QT"] = df.loc[i, "Qs"] + df.loc[i, "Qg"]
#         df.loc[i,"Qt"] = 5.6
#     else:
#         df.loc[i,"Qs"] = df.loc[i,"RS"] * U
#         df.loc[i,"Qg"] = CG * df.loc[i-1, "Qg"] + (1-CG) * df.loc[i,"RG"] * U
#         df.loc[i,"QT"] = df.loc[i,"Qs"] + df.loc[i,"Qg"]
#         df.loc[i,"Qt"] = CR * df.loc[i-1,"Qt"] + (1 - CR) * df.loc[i-L,"QT"]

三水源汇流

# 计算三水源汇流
    if i == 0:
        df.loc[i,"QS"] = df.loc[i,"RS"] * U
        df.loc[i,"QI"] = 1/3 * Q
        df.loc[i,"QG"] = 1/3 * Q
        df.loc[i, "QT"] = df.loc[i, "QS"] + df.loc[i, "QI"] + df.loc[i, "QG"]
    if i >0:
        df.loc[i, "QS"] = df.loc[i, "RS"] * U
        df.loc[i,"QI"] = CI * df.loc[i-1,"QI"] + (1-CI)*df.loc[i,"RI"]*U
        df.loc[i,"QG"] = CG* df.loc[i-1,"QG"]+(1-CG)*df.loc[i,"RG"] *U
        df.loc[i,"QT"] = df.loc[i,"QS"]+df.loc[i,"QI"]+df.loc[i,"QG"]
    if i >=0 and i <= L:
        df.loc[i,"Qt"] = Q
    if i >= L+1:
        df.loc[i,"Qt"] = df.loc[i-1,"Qt"]*CR+(1-CR) * df.loc[i-L,"QT"]

最后写入excel表，保留4位小数

ds = pd.DataFrame(df)
ds.to_excel(r"C:\Users\Admin\Desktop\新安江1.xlsx",float_format='%.4f')

最后三水源结果

写到这感觉自己写的太麻烦了，希望以后能够改进吧

这边蒸散发还是借鉴了一个哥们的代码，链接在这边：(26条消息) Python数据分析实例，利用Pandas建立流域三层蒸发和蓄满产流模型_三层蒸发模型代码_Joyonlyonly的博客-CSDN博客