流況曲線作図(改訂版)
流況曲線作図プログラムを(自分にとって)少しわかりやすくしたのでアップ。 対象期間内に欠測はないことを前提にしている。
import numpy as np import pandas as pd import datetime import matplotlib.pyplot as plt from scipy import interpolate def rdata(fnameR): df=pd.read_csv(fnameR,header=0,index_col=0) df.index = pd.to_datetime(df.index, format='%Y/%m/%d') df=df['2001/01/01':'2018/12/31'] return df def makedata(qq,ylist): qqq=np.zeros((101,len(ylist)+1),dtype=np.float64) # discharge ppp=np.arange(101,dtype=np.float64) # propability for i,sy in enumerate(ylist): qd=np.array(qq[sy]) qd=np.sort(qd)[::-1] ii=np.arange(len(qd)) pp=ii/(len(qd)-1)*100 ff=interpolate.interp1d(pp,qd) qqq[:,i]=ff(ppp) for j in range(len(ppp)): qqq[j,-1]=np.average(qqq[j,0:len(ylist)]) return ppp,qqq def drawfig1(ylist,ppp,qqq,tstr,fnameF): fsz=16 xmin=0; xmax=100; dx=10 ymin=0; ymax=200; dy=20 fig=plt.figure(figsize=(10,6),facecolor='w') plt.rcParams['font.size']=fsz plt.rcParams['font.family']='sans-serif' plt.xlim([xmin,xmax]) plt.ylim([ymin,ymax]) plt.xlabel('Non-exceedance probability (%)') plt.ylabel('River flow Q (m$^3$/s)') plt.xticks(np.arange(xmin,xmax+dx,dx)) plt.yticks(np.arange(ymin,ymax+dy,dy)) plt.grid(color='#999999',linestyle='solid') plt.title(tstr,loc='left',fontsize=fsz-2) for i in range(0,len(ylist)): plt.plot(ppp,qqq[:,i],'-',color='#aaaaaa',lw=1) plt.plot(ppp,qqq[:,-1],'-',color='#ff0000',lw=3,label='Average') # Q95 q95=qqq[95,-1] plt.plot([95],[q95],'o',color='#ff0000') plt.plot([95,95],[q95,q95+100],'-',color='#ff0000') xs=95 ys=q95+100 ss='Q95%:{0:6.2f}m$^3$/s'.format(q95) plt.text(xs,ys,ss,va='top',ha='right',fontsize=fsz-2,rotation=90) plt.legend() plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1) plt.show() def tblq(lyy,qym): mmq=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec','Ave'] tstr='* monthly average discharge' mmm=mmq arr=qym print(tstr) ss='{0:4s}'.format('Year') for m in mmm: ss=ss+',{0:>6s}'.format(m) print(ss) lyy=lyy+['Ave.'] for i,yy in enumerate(lyy): ss='{0:4s}'.format(yy) for j in range(len(mmm)): ss=ss+',{0:6.1f}'.format(arr[i,j]) print(ss) def emon(lyy,qq): lmm=['01','02','03','04','05','06','07','08','09','10','11','12'] # month qym=np.zeros((len(lyy)+1,len(lmm)+1),dtype=np.float64) kda=np.zeros((len(lyy)+1,len(lmm)+1),dtype=np.float64) for i,sy in enumerate(lyy): for j,sm in enumerate(lmm): ss=sy+'-'+sm qym[i,j]=qq[ss].sum() # sum of each month kda[i,j]=qq[ss].count() # availavle days qym[i,-1]=np.sum(qym[i,:]) # sum of a year of data kda[i,-1]=np.sum(kda[i,:]) # sum of a year of available days for j in range(len(lmm)+1): qym[-1,j]=np.sum(qym[:,j]) # sum of months of data kda[-1,j]=np.sum(kda[:,j]) # sum of months of available days qym=qym/kda return qym def main(): lyy=['2001','2002','2003','2004','2005','2006','2007','2008','2009', '2010','2011','2012','2013','2014','2015','2016','2017','2018'] for kkk1 in [1,2]: if kkk1==1: fnameR='df_teromu_estimated_tk.csv' df0=rdata(fnameR) qd=np.array(df0['Q1516'])*822.6/1070*3393/3417 # discharge at dam site tstr='Duration curve at dam site estimated by tank model (2001-2018)' fnameF='fig_duration_dam_tank.png' if kkk1==2: fnameR='df_teromu_estimated_ml.csv' df0=rdata(fnameR) qd=np.array(df0['Q_knn'])*822.6/1070*3393/3417 # discharge at dam site tstr='Duration curve at dam site estimated by KNN (2001-2018)' fnameF='fig_duration_dam_knn.png' qq=pd.Series(qd,index=df0.index) qym=emon(lyy,qq) tblq(lyy,qym) ppp,qqq=makedata(qq,lyy) drawfig1(lyy,ppp,qqq,tstr,fnameF) #============== # Execution #============== if __name__ == '__main__': main()
以 上