damyarou

python, GMT などのプログラム

流況曲線作図(改訂版)

流況曲線作図プログラムを(自分にとって)少しわかりやすくしたのでアップ。 対象期間内に欠測はないことを前提にしている。

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()

以 上