python生成芝加哥雨型的简单方法

import pandas as pd
import numpy as np
import os
import math
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker

os.chdir(r'E:\HANproject\SWMMproject')

#雨强计算定义为函数，后续需要反复计算，公式参数来源于北京雨强设计标准
def intensity(t,p):
    a = 1602*0.36*(1+1.032*math.log(p,10))
    b = 11.593
    c = 0.681
    its = a*((1-c)*t+b)/math.pow(t+b,c+1)
    return its

t = np.arange(0,120)#时间序列
itAr = np.zeros((120,5))#雨强列表
xa = [] #具有时间格式的时间序列
for i in t:
    if i<60:
        if i<10:
            temp = '00:0'+str(i)
        else:
            temp = '00:'+str(i)
        xa.append(temp)
    elif i>=60:
        if i<70:
            temp = '01:0'+str(i-60)
        else:
            temp = '01:'+str(i-60)
        xa.append(temp)

#合成芝加哥雨型，雨峰位置在降雨开始后48 min，即雨峰系数为0.4
def rainCalc(t,p,j):
    for i in range(1,49):
        itAr[i-1,j] = intensity((48-i)/0.4,p)/60
    for i in range(48,121):
        itAr[i-1,j] = intensity((i-48)/0.4,p)/60
    plt.gca().xaxis.set_major_locator(ticker.MultipleLocator(20))
    plt.plot(xa, itAr[:,j], label=str(p)+'a')
    plt.margins(0.05,0.1)

#重现期列表，循环生成多种重现期降雨
pa = (2,5,20,50,100)
for i in range(0,5):
    rainCalc(t,pa[i],i)

#绘图展示雨型
plt.xlabel('Time')
plt.ylabel('Rain intensity/mm/min')
plt.legend()
plt.show()

rainIntense = pd.DataFrame(itAr)
rainIntense.insert(0,'time',xa)
rainIntense.insert(0,'date','01/01/2024')
ar = np.zeros((120,5))#npzeros有俩括号
for i in range(0,5):
    for j in t:
        ar[j,i] = itAr[:j,i].sum()
#对第i列，第j个及之前的元素求和，得到cumulative rainfall输入模型

rainGage = pd.DataFrame()
name = ['2a','5a','20a','50a','100a']
for i in range(0,5):
    rainGage.insert(i,name[i],ar[:,i])#iloc()列索引
rainGage.insert(0,'date',rainIntense.iloc[:,0])
rainGage.insert(1,'time',rainIntense.iloc[:,1])
rainGage.to_csv('rainGage.csv',float_format='%.2f')#写入文件，格式保留两位小数


    
    
    

关于芝加哥雨型相关公式推导的前置知识，可以自行搜索~