matplotlib 河川流況表示（１）
概要
河川の流況を示す３種類の作図プログラムを紹介します。
流況図
縦軸に日平均流量、横軸に超過確率（厳密には超過率）をとった年間の流況を示す図です。 通常、日本では横軸は「日」ですが、海外（東南アジア？）では横軸は「確率」としています。
ここでは、１年毎に８年分の線を引き、赤い太線で８年間の平均を示しています。 図の下に、各年の平均・最小・最大値を画面出力したものをのせています。
------------------------- Year ave min max ------------------------- 1992 59.7 11.2 216.4 1993 81.8 8.4 235.1 1994 21.2 2.7 106.8 1995 68.7 15.3 231.5 1996 36.6 12.6 161.6 1997 27.8 10.0 194.7 1998 40.5 15.1 152.2 1999 21.1 1.7 133.4 ------------------------- Ave. 44.7 1.7 235.1
月平均流量ヒートマップ
各年の月平均流量をマトリックス状に配置したヒートマップです。 ヒートマップの解説は、以下のURLに詳しいです。
月平均流量時系列図
月平均流量の時系列図を年毎にマトリックス状に示したものです。 流量、雨量、気温など季節変動するものは、横長のグラフを作るよりは、年ごとに示したほうが季節変動の傾向がわかりやすいと思います。
流況図作図プログラム
元データはエクセルファイルであり、以下に示すように、年ごとにシートに格納されていますが、１年の値は横軸に月、縦軸に日付としたマトリックスに日平均流量が格納されています。 月の日数は月によって異なるので、データが無い日付の欄は空欄にしています。
データを読み込む関数：rdata() では、上記エクセルファイルから読み取ったデータを以下のような、日付と流量値を格納した Series として出力します。 元データの格納書式が異なっても、プログラム内で同一書式に変えておけば、統計処理や作図プログラムは、同じものを流用することができます。
1992-01-01 118.3 1992-01-02 112.5 1992-01-03 134.2 ..... .... 1999-12-29 15.7 1999-12-30 13.3 1999-12-31 12.3
プログラムではscipy.interpolateをインポートしています。
これは、横軸を決める際、１年365日あるいは366日分の計測値を、0%（計測最大値）から100%（計測最小値）まで1%ピッチで当てはめるために用いています。
こうすることにより、うるう年の扱いにも困りませんし、プロット期間全年の平均を取る場合も各年の計測値が同じパーセンテージに当てはめられているので楽ちんです。
プログラム全文
import numpy as np import pandas as pd import datetime import matplotlib.pyplot as plt from scipy import interpolate # global tstr='Duration Curve at Kelaena Dam Site (1992-1999 by FS-report)' ylist=['1992','1993','1994','1995','1996','1997','1998','1999'] # sheet name def rdata(): fnameR='xls_qq_dam.xlsx' # input file name q=[] # discharge d=[] # date k=0 # counter for date dini=datetime.datetime(1991,12,31) # initial date for iii,sname in enumerate(ylist): df=pd.read_excel(fnameR, sheet_name=sname) # read excel data for i in range(1,13): # loop of month for j in range(0,31): # loop of date if df.iloc[j,i]==df.iloc[j,i]: # ignore nan q=q+[float(df.iloc[j,i])] k=k+1 dnow=dini+datetime.timedelta(days=k) d=d+[dnow] qq=pd.Series(np.array(q),index=d) return qq def makedata(qq): qqq=np.zeros((101,len(ylist)+1),dtype=np.float64) # discharge ppp=np.arange(101,dtype=np.float64) # probability sta=np.zeros((len(ylist)+1,3),dtype=np.float64) # ave, min, max of a year 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) sta[i,0]=np.average(qd) sta[i,1]=np.min(qd) sta[i,2]=np.max(qd) sta[-1,0]=np.average(qq) sta[-1,1]=np.min(qq) sta[-1,2]=np.max(qq) for j in range(len(ppp)): qqq[j,-1]=np.average(qqq[j,:]) return ppp,qqq,sta def drawfig(ppp,qqq): 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) lb=ylist ls=['-','--','-.','-','--','-.','-','--'] col=[ '#000000', '#000000', '#000000', '#0000ff', '#0000ff', '#0000ff', '#00ff00', '#00ff00' ] for i,sy in enumerate(ylist): plt.plot(ppp,qqq[:,i],linestyle=ls[i],color=col[i],lw=2,label=lb[i]) plt.plot(ppp,qqq[:,-1],'-',color='#ff0000',lw=3,label='Average') plt.legend() fnameF='fig_duration1.jpg' plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1) plt.show() def main(): qq=rdata() ppp,qqq,sta=makedata(qq) drawfig(ppp,qqq) print('-------------------------') print('{0:>4s} {1:>6s} {2:>6s} {3:>6s}'.format('Year','ave','min','max')) print('-------------------------') for i,sy in enumerate(ylist): print('{0:4s} {1:6.1f} {2:6.1f} {3:6.1f}'.format(sy,sta[i,0],sta[i,1],sta[i,2])) print('-------------------------') print('{0:4s} {1:6.1f} {2:6.1f} {3:6.1f}'.format('Ave.',sta[-1,0],sta[-1,1],sta[-1,2])) #============== # Execution #============== if __name__ == '__main__': main()
月平均流量ヒートマップおよび月平均流量時系列図作図プログラム
元データは流況図用のものと同じです。 エクセルより読み取ったデータを以下のような pandas の DataFrame に書き直し、この数値をそのままヒートマップおよび時系列図にプロットしています。
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov \
Year
1992 61.6 82.0 113.9 128.1 34.1 46.8 31.5 34.7 38.2 20.9 14.4
1993 161.2 114.7 133.6 149.7 88.5 51.0 61.0 72.7 104.1 13.0 16.2
1994 26.5 17.8 41.3 48.9 33.1 30.2 16.1 8.6 9.3 5.5 12.8
1995 65.8 66.8 95.1 82.5 60.6 61.9 56.3 72.6 85.1 51.6 95.8
1996 43.1 54.0 37.9 34.5 44.3 36.0 39.4 32.8 18.8 30.2 35.3
1997 30.1 52.1 32.0 31.6 25.7 17.2 32.8 20.2 13.3 15.4 13.2
1998 36.1 44.8 35.9 41.6 48.3 36.1 35.1 30.8 35.4 40.7 47.7
1999 19.1 44.4 41.9 21.0 28.7 24.5 14.2 8.2 2.7 12.9 22.7
Ave. 55.4 59.6 66.5 67.2 45.4 38.0 35.8 35.1 38.4 23.8 32.3
Dec Ave.
Year
1992 111.5 59.7
1993 19.4 81.8
1994 4.9 21.2
1995 31.2 68.7
1996 33.7 36.6
1997 51.4 27.8
1998 54.1 40.5
1999 14.5 21.1
Ave. 40.1 44.7
プログラム全文
import numpy as np import pandas as pd import datetime import matplotlib.pyplot as plt import seaborn as sns # global ylist=['1992','1993','1994','1995','1996','1997','1998','1999'] def rdata(): fnameR='xls_qq_dam.xlsx' # input file name flist=['1992','1993','1994','1995','1996','1997','1998','1999'] # sheet name q=[] # discharge d=[] # date k=0 # counter for date dini=datetime.datetime(1991,12,31) # initial date for iii,sname in enumerate(flist): df=pd.read_excel(fnameR, sheet_name=sname) # read excel data for i in range(1,13): # loop of month for j in range(0,31): # loop of date if df.iloc[j,i]==df.iloc[j,i]: # ignore nan q=q+[float(df.iloc[j,i])] k=k+1 dnow=dini+datetime.timedelta(days=k) d=d+[dnow] qq=pd.Series(np.array(q),index=d) return qq def makedf(qq): mlist=['-01','-02','-03','-04','-05','-06','-07','-08','-09','-10','-11','-12'] n=len(ylist) m=len(mlist) qt=np.zeros((n+1,m+1),dtype=np.float64) # sum of discharges of a month ct=np.zeros((n+1,m+1),dtype=np.int) # sum of days of a manth for i,sy in enumerate(ylist): for j,sm in enumerate(mlist): sym=sy+sm qd=np.array(qq[sym]) qt[i,j]=np.sum(qd) ct[i,j]=len(qd) for i in range(n): qt[i,m]=np.sum(qt[i,0:m]) ct[i,m]=np.sum(ct[i,0:m]) for j in range(m): qt[n,j]=np.sum(qt[0:n,j]) ct[n,j]=np.sum(ct[0:n,j]) qt[n,m]=np.sum(qq) ct[n,m]=len(qq) df = pd.DataFrame({'Year': ylist+['Ave.'], 'Jan' : qt[:,0]/ct[:,0], 'Feb' : qt[:,1]/ct[:,1], 'Mar' : qt[:,2]/ct[:,2], 'Apr' : qt[:,3]/ct[:,3], 'May' : qt[:,4]/ct[:,4], 'Jun' : qt[:,5]/ct[:,5], 'Jul' : qt[:,6]/ct[:,6], 'Aug' : qt[:,7]/ct[:,7], 'Sep' : qt[:,8]/ct[:,8], 'Oct' : qt[:,9]/ct[:,9], 'Nov' : qt[:,10]/ct[:,10], 'Dec' : qt[:,11]/ct[:,11], 'Ave.': qt[:,12]/ct[:,12]}) df=df.set_index('Year') return df def fighmap(df): wid=0.5*(len(ylist)+1) fsz=14 plt.figure(figsize=(12,wid),facecolor='w') plt.rcParams['font.size']=fsz plt.rcParams['font.family']='sans-serif' # cm='YlGnBu' cm='Blues' sns.heatmap(df,annot=True,fmt='.1f',vmin=0,vmax=160,cmap=cm,cbar=False,linewidths=0.5) plt.xlabel('Month') plt.ylabel('Year') plt.yticks(rotation=0) plt.title('Monthly Average Discharge (m$^3$/s)',loc='left') fnameF='fig_hmap1.jpg' plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1) plt.show() def figts(df): # arrange of subplot # ------------------ # i=0 i=1 i=2 i=3 # i=4 i=5 i=6 i=7 # ------------------ icol=4 irow=2 fsz=11 xmin=0;xmax=12 ymin=0;ymax=200 x=np.array([1,2,3,4,5,6,7,8,9,10,11,12]) x=x-0.5 mlist=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'] emp=['']*12 plt.figure(figsize=(12,12/icol*irow),facecolor='w') plt.rcParams['font.size']=fsz plt.rcParams['font.family']='sans-serif' plt.suptitle('Monthly Average Discharge in Each Year', x=0.5,y=0.92,fontsize=fsz+4) plt.subplots_adjust(wspace=0.07,hspace=0.07) for i,sy in enumerate(ylist): num=i+1 plt.subplot(2,4,num) plt.xlim([xmin,xmax]) plt.ylim([ymin,ymax]) if i==0 or i==4: plt.ylabel('Discharge (m$^3$/s)') else: plt.yticks([]) if 4<=i: plt.xticks(x,mlist,rotation=90) else: plt.xticks(x,emp,rotation=90) y=df.loc[sy][0:12] qa=df.loc[sy][12] plt.plot(x,y,'o-',color='#000080',lw=1) plt.plot([xmin,xmax],[qa,qa],'--',color='#ff0000',lw=2) plt.text(8.5,qa,'Qave',fontsize=fsz,va='bottom',ha='center', color='#ff0000') s1='Year: '+sy s2='Qave: {0:.1f}'.format(qa) textstr=s1+'\n'+s2 props = dict(boxstyle='square', facecolor='#ffffff', alpha=1) xs=xmin+0.95*(xmax-xmin); ys=ymin+0.95*(ymax-ymin) plt.text(xs,ys, textstr, fontsize=fsz+2,va='top',ha='right', bbox=props) fnameF='fig_series1.jpg' plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1) plt.show() def main(): qq=rdata() df=makedf(qq) pd.options.display.precision = 1 print(df) fighmap(df) figts(df) #============== # Execution #============== if __name__ == '__main__': main()
