damyarou

python, GMT などのプログラム

Python 効率カーブ(表付きグラフ)

簡単な表と一体化した画像を作った。

作例

f:id:damyarou:20200516115000j:plain

このグラフにおけるTips

作図領域を上下に分ける

drawfig1() は下側のグラフ(作図領域の下側70%)、```drawfig2()''' は上川の表を描画する。

    plt.axes((0.0, 0.0, 1.0, 0.7)); drawfig1()
    plt.axes((0.0, 0.7, 1.0, 0.3)); drawfig2()

白抜き記号描画

  • ms=8 :見号の大きさ
  • color='#000080' :記号の色
  • markerfacecolor='#ffffff' :白抜き
  • markeredgewidth=1 :記号輪郭の太さ
    plt.plot(q_hmax,e_hmax,'s',ms=8,color='#000080',markerfacecolor='#ffffff',markeredgewidth=1,label='H=192.5m')

プログラム全文

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq


def reg(x,y):
    def func(param,x,y):
        a=param[0]
        b=param[1]
        c=param[2]
        return y - (a*(x-b)**2 + c)
    para0 =np.array([-1.0, 20.96, 89.65])
    res = leastsq(func, para0, args=(x,y))
    a=res[0][0]
    b=res[0][1]
    c=res[0][2]
    return a,b,c


def drawfig1():
    q_hmax=np.array([22.37,21.67,20.26,19.50,17.77,17.34,16.32,15.17,15.01,13.67,11.27])
    e_hmax=np.array([89.48,89.45,89.23,88.96,88.18,87.87,87.14,86.23,86.07,84.98,81.82])

    q_hnor=np.array([23.94,23.05,20.96,19.28,18.68,17.27,16.77,15.83,14.67,13.46,10.90])
    e_hnor=np.array([89.53,89.64,89.65,89.17,88.94,88.10,87.73,87.04,85.97,84.86,81.53])
    
    q_hmin=np.array([24.24,22.50,20.46,18.90,18.41,17.18,16.36,15.89,14.79,13.45,10.64])
    e_hmin=np.array([89.27,89.62,89.60,89.12,88.80,88.05,87.39,86.99,85.91,84.78,81.07])

    fsz=14
    xmin=5 ; xmax=25; dx=5
    ymin=75; ymax=95; dy=5
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Turbine discharge per unit Qu (m$^3$/s)')
    plt.ylabel('Combined efficiency ef (%)')
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(color='#999999',linestyle='solid')
    #plt.gca().set_aspect('equal',adjustable='box')
    #plt.title(tstr,loc='left',fontsize=fsz)

    x=q_hnor
    y=e_hnor
    a,b,c=reg(x,y)
    qmax=24.5
    qmin=qmax*0.35
    xx=np.arange(qmin,qmax,0.005)
    yy=a*(xx-b)**2+c   
    ss='ef={0:6.3f}*(Qu-{1:6.3f})**2+{2:6.3f}'.format(a,b,c)
    print(np.min(yy))
    print(np.max(yy))
    
    plt.plot(xx,yy,'-',color='#000080',lw=2,label=ss)
    plt.plot(q_hmax,e_hmax,'s',ms=8,color='#000080',markerfacecolor='#ffffff',markeredgewidth=1,label='H=192.5m')
    plt.plot(q_hnor,e_hnor,'o',ms=8,color='#000080',markerfacecolor='#ffffff',markeredgewidth=1,label='H=180.0m')
    plt.plot(q_hmin,e_hmin,'^',ms=8,color='#000080',markerfacecolor='#ffffff',markeredgewidth=1,label='H=171.5m')
    
    plt.plot([qmax,qmax],[ymin,ymax],'--',color='#000080',lw=2)
    plt.plot([qmin,qmin],[ymin,ymax],'--',color='#000080',lw=2)
    xs=qmax; ys=ymin+0.02*(ymax-ymin); ss='Qu(max)={0:.1f}m$^3$/s'.format(qmax)
    plt.text(xs,ys,ss,ha='right',va='bottom',fontsize=fsz,rotation=90)
    xs=qmin; ys=ymin+0.02*(ymax-ymin); ss='Qu(min)=Qu(max) x 0.35'
    plt.text(xs,ys,ss,ha='right',va='bottom',fontsize=fsz,rotation=90)
    plt.legend(loc='upper left',shadow=True,fontsize=fsz-2)

    
def drawfig2():
    tstr='Design head setting (Q=49.0m$^3$/s)'
    row4=['Item','RWL','TWL','Margin','head loss','Value of head']
    row3=['Maximum head','EL.284.5','EL.92.0','0.0m','0.0m', 'H=192.5m']
    row2=['Nominal head','EL.284.5','EL.92.0','1.0m','11.5m','H=180.0m']
    row1=['Minimum head','EL.275.0','EL.92.0','0.0m','11.5m','H=171.5m']
    x=np.array([0.5,1.5,2.5,3.5,4.5,5.5])
    y=np.array([0.5,1.5,2.5,3.5])    
    fsz=12
    xmin=0 ; xmax=6
    ymin=-0.5 ; ymax=5
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.axis('off')
    for i in range(len(x)):
        for j in range(len(y)):
            if j==0: _row=row1
            if j==1: _row=row2
            if j==2: _row=row3
            if j==3: _row=row4
            ss=_row[i]
            plt.text(x[i],y[j],ss,va='center',ha='center',fontsize=fsz)
    plt.text(0.5*(xmin+xmax),4.5,tstr,va='center',ha='center',fontsize=fsz)
    for yy in [0,1,2,3,4]:
        plt.plot([xmin,xmax],[yy,yy],'-',color='#000000',lw=1)
    for xx in [0,1,2,3,4,5,6]:
        plt.plot([xx,xx],[0,4],'-',color='#000000',lw=1,clip_on=False)


def main():
    fnameF='fig_efficiency.jpg'
    plt.figure(figsize=(8,6),facecolor='w')
    plt.rcParams['font.family'] ='sans-serif'
    plt.rcParams['font.size']=14
    plt.axes((0.0, 0.0, 1.0, 0.7)); drawfig1()
    plt.axes((0.0, 0.7, 1.0, 0.3)); drawfig2()
    plt.tight_layout()
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()
    
    
#---------------
# Execute
#---------------
if __name__ == '__main__': main()

以 上

Python 月次・年次データ集計(雨量)

日雨量を月・年毎に合計したもの。 各月の雨量の単位は mm/momth (月累計雨量)、 各年の雨量の単位は mm/year(年累計雨量)である。

#==============================
# Rainfall
# Dayly => Monthly, Yearly data
#==============================
import pandas as pd
import numpy as np
import datetime


def rdata():
    fnameR='xls_rdata_20190630.xlsx' # input file name
    sname='daily'
    df=pd.read_excel(fnameR, sheet_name=sname,header=0, index_col=0) # read excel data
    df.index = pd.to_datetime(df.index, format='%Y/%m/%d')
    return df    
    

def cal_y(rf,lyyy):
    rfy=np.zeros(len(lyyy),dtype=np.float)
    for i,ss in enumerate(lyyy):
        rfy[i]=rf[ss].sum() # yearly rainfall
    return rfy


def cal_m(rf,lyyy,lmmm):
    rfm=np.zeros((len(lyyy),len(lmmm)+1),dtype=np.float)
    for i,yy in enumerate(lyyy):
        for j,mm in enumerate(lmmm):
            ss=yy+'/'+mm
            try:
                rfm[i,j]=rf[ss].sum() # monthly rainfall
            except KeyError:
                rfm[i,j]=np.nan
        rfm[i,-1]=np.sum(rfm[i,:])
    return rfm


def rf_year(lyy,la,df0):
    rfy=np.zeros((len(lyy),len(la)),dtype=np.float64)
    for i,ss in enumerate(la):
        rf=pd.Series(df0[ss], index=df0.index)
        rfy[:,i]=cal_y(rf,lyy)

    ss='{0:4s}'.format('Year')
    for k in [1,2,3,4,5,6,7]:
        ss=ss+',{0:>8s}'.format('Area-'+str(k))
    ss=ss+',{0:>8s}'.format('Dam')
    ss=ss+',{0:>8s}'.format('Teromu')
    print(ss)
    for j,yy in enumerate(lyy):
        ss='{0:4s}'.format(yy)
        for i in range(len(la)):
            ss=ss+',{0:8.0f}'.format(rfy[j,i])
        print(ss)               
    ss='{0:4s}'.format('Ave.')
    for i in range(len(la)):
        ss=ss+',{0:8.0f}'.format(np.mean(rfy[:,i]))
    print(ss)


def rf_mon(lyy,lmm,rf):
    rfm=cal_m(rf,lyy,lmm)
    ss='{0:4s}'.format('Year')
    for m in ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aig','Sep','Oct','Nov','Dec','Sum']:
        ss=ss+',{0:>6s}'.format(m)
    print(ss)
    for i,yy in enumerate(lyy):
        ss='{0:4s}'.format(yy)
        for j in range(len(lmm)+1):
            ss=ss+',{0:6.0f}'.format(rfm[i,j])
        print(ss)               
    ss='{0:4s}'.format('Ave.')
    for j in range(len(lmm)+1):
        ss=ss+',{0:6.0f}'.format(np.nanmean(rfm[:,j]))
    print(ss)


def main():
    df0=rdata() # daily rainfall at each area (2000.03-2019.06)
    lyy=['2001','2002','2003','2004','2005','2006','2007','2008','2009',
         '2010','2011','2012','2013','2014','2015','2016','2017','2018']
    la=['area1','area2','area3','area4','area5','area6','area7','Dam','Teromu']
    rf_year(lyy,la,df0)
    
    lyy=['2000','2001','2002','2003','2004','2005','2006','2007','2008','2009',
         '2010','2011','2012','2013','2014','2015','2016','2017','2018','2019']
    lmm=['01','02','03','04','05','06','07','08','09','10','11','12']

    print('* Monthly rainfall at Dam')
    rf=pd.Series(df0['Dam'], index=df0.index)
    rf_mon(lyy,lmm,rf)

    print('* Monthly rainfall at Teromu')
    rf=pd.Series(df0['Teromu'], index=df0.index)
    rf_mon(lyy,lmm,rf)

        
#==============
# Execution
#==============
if __name__ == '__main__': main()

以 上

Python 月次及び年次データ集計(流量)

概要

日平均流量の月次および年次データ集計のプログラム。

ここでは流量の集計を行い、月平均・年平均を算出している。 年平均は、「月平均の平均」ではなく、「1年間全データの平均」としている。 これは365日の流況図から出した年平均と数値を合わせるため。 「月平均の平均」と「年間全データの平均」は結構違う場合があるので注意。

計算のルール

入力

データ取得開始日の月初めから、データ取得完了日の月末まで、データの有無に関わらず、全ての日付が連続的に入力されていること。

出力

  • 各月の中で1日でも欠測があれば、その月の平均は nan とする
  • 各年の平均は、各年の全データ(365日あるいは366日)の平均とする。その年に1日でも欠測があればその年の平均は nan とする。
  • 各月の集計平均(最終行)は、欠測を含まない月の全データの平均値とする。欠測がある月のデータは月平均の計算に用いない。

出力例

  • 年別に各月の平均日流量を表示している。
  • 最終列はその年の平均日流量。
  • 最終行は各月の年をまたいだ平均日流量。
  • Wordに貼り付けて数表を作るため、画面にCSVイメージで出力する。
  • 各月の中で、1日でも欠測があれば月平均は、nanとしている。
  • nan の月を含む年の平均はもちろん nan である。
Year,   Jan,   Feb,   Mar,   Apr,   May,   Jun,   Jul,   Aig,   Sep,   Oct,   Nov,   Dec,  Ave.
2013,  70.6,  83.3, 118.3, 124.3, 117.4, 115.6,  93.4,  53.2,  41.1,  28.4,  76.7, 113.9,  86.3
2014,  46.9,  36.5,  59.2,  71.9,  42.1,  86.1,  79.0,  50.6,  28.1,   nan,  23.0,   nan,   nan
2015,  47.5,  95.7, 101.0,  90.9,  74.9,  37.8,  32.0,  27.8,  15.4,  14.0,   9.1,  37.0,  48.3
2016,  38.2,  72.3, 125.8, 145.9,  72.2,  80.1,  49.4,  54.0,  49.7,  92.2, 112.1,  54.1,  78.7
2017,  50.7,  63.0,  99.2, 108.1, 128.7, 105.6,  86.3,   nan,   nan,   nan,   nan,   nan,   nan
2018,   nan,   nan,   nan,   nan,  84.1,   nan,   nan,   nan,  31.9,  32.6,  37.4,  64.7,   nan
2019,  55.2,  91.4,  68.4, 105.9,  88.5,   nan,   nan,   nan,   nan,   nan,   nan,   nan,   nan
Ave.,  51.5,  73.7,  95.3, 107.8,  86.8,  85.0,  68.0,  46.4,  33.2,  41.8,  51.7,  67.4,  71.1

プログラム

#====================================
# Data aggregation (monthly & yearly)
#====================================
import pandas as pd
import numpy as np
import datetime


def rdata():
    fnameR='df_qq2.csv' # input file name
    df=pd.read_csv(fnameR, header=0, index_col=0) # read excel data
    df.index = pd.to_datetime(df.index, format='%Y/%m/%d')
    return df    
    

def cal_m(rf,lyyy,lmmm):
    qqm=np.zeros((len(lyyy)+1,len(lmmm)+1),dtype=np.float64)
    kda=np.zeros((len(lyyy)+1,len(lmmm)+1),dtype=np.float64)
    for i,yy in enumerate(lyyy):
        for j,mm in enumerate(lmmm):
            ss=yy+'/'+mm
            try:
                qqm[i,j]=rf[ss].sum() # sum in a month
                kda[i,j]=rf[ss].count() # availavle days
                nac=np.count_nonzero(np.isnan(rf[ss])) # unavailable days
            except KeyError:
                qqm[i,j]=np.nan
                kda[i,j]=np.nan
            if 0<nac:
                qqm[i,j]=np.nan
                kda[i,j]=np.nan
        qqm[i,-1]=np.sum(qqm[i,:]) # sum in a year
        kda[i,-1]=np.sum(kda[i,:]) # sum in a year of available days
    for j in range(len(lmmm)+1):
        qqm[-1,j]=np.nansum(qqm[:,j])
        kda[-1,j]=np.nansum(kda[:,j])
    qmean=qqm/kda
    return qmean


def qq_mon(lyy,lmm,qq):
    qqm=cal_m(qq,lyy,lmm)
    ss='{0:4s}'.format('Year')
    for m in ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aig','Sep','Oct','Nov','Dec','Ave.']:
        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(lmm)+1):
            ss=ss+',{0:6.1f}'.format(qqm[i,j])
        print(ss)               


def main():
    df0=rdata() # daily rainfall at each area (2000.03-2019.06)
    lyy=['2013','2014','2015','2016','2017','2018','2019']
    lmm=['01','02','03','04','05','06','07','08','09','10','11','12']

    print('* Monthly discharge at Teromu GS')
    qq=pd.Series(df0['q_tot'], index=df0.index)
    qq_mon(lyy,lmm,qq)

        
#==============
# Execution
#==============
if __name__ == '__main__': main()

以 上

Python バージョンを入れ替える(Mac)

2020年4月27日、Pythonは3.8.1を入れていたが、3.8.2にしてみた。 Pythonと関連ライブラリのインストールはpyenvとpipのみで行っている。

目的はPythonのバージョンを変えるだけなので、.zshrcに記載しているパスや、Jupyterの設定は変更しない。もちろんこれらのファイルはいじらない。

関連事項として以下を参照。

パスの設定:.zshrc

PROMPT="%# "

export PYENV_ROOT=${HOME}/.pyenv
export PATH=${PYENV_ROOT}/bin:$PATH
eval "$(pyenv init -)"

alias brew="env PATH=${PATH/\/Users\/${USER}\/\.pyenv\/shims:?/} brew"

ファインダーに隠しファイルを表示するスクリプトa_dota.sh

defaults write com.apple.finder AppleShowAllFiles TRUE
killall Finder

ファインダーで隠しファイルを非表示にするスクリプトa_dotk.sh

defaults write com.apple.finder AppleShowAllFiles FALSE
killall Finder

以下に手順を示す。

まずは、.pyenv を含む隠しファイルをファインダーに表示する。

% sh a_dota.sh
% brew uninstall pyenv
Uninstalling /usr/local/Cellar/pyenv/1.2.18... (695 files, 2.5MB)

上記操作により、homebrew で pyenv を削除後、ディレクト.pyenvをファインダーで削除。

homebrew で pyenv を再インストール。 以下の画面では、homebrew をアップデート後、pyenvをインストールしている。

% brew install pyenv
Updating Homebrew...
==> Auto-updated Homebrew!
Updated 2 taps (homebrew/core and homebrew/cask).
==> New Formulae
autodiff                   hcxtools                   libcpuid
==> Updated Formulae
ammonite-repl              fastlane                   nushell
azcopy                     fetchmail                  open-mesh
(途中省略)
fantastical                qiyimedia
==> Deleted Casks
go2shell                                 wineskin-winery

==> Downloading https://homebrew.bintray.com/bottles/pyenv-1.2.18.catalina.bottl
Already downloaded: /Users/kk/Library/Caches/Homebrew/downloads/cd594a21b931d84b2c968913ea285837c751aa1c0db7d47a935f60c49794d6cd--pyenv-1.2.18.catalina.bottle.tar.gz
==> Pouring pyenv-1.2.18.catalina.bottle.tar.gz
🍺  /usr/local/Cellar/pyenv/1.2.18: 695 files, 2.5MB

pyenvのインストールが無事完了。

欲しいPythonのバージョンを確認。

% pyenv install -l
Available versions:
  2.1.3
  2.2.3
(途中省略)
  3.7.7
  3.8.0
  3.8-dev
  3.8.1
  3.8.2    # <= これが欲しい!
  3.9.0a5
  3.9-dev
  activepython-2.7.14
  activepython-3.5.4
  activepython-3.6.0
  anaconda-1.4.0
  anaconda-1.5.0
(途中省略)
  stackless-3.4.7
  stackless-3.5.4

Python 3.8.2 をインストール。

% pyenv install 3.8.2
python-build: use openssl@1.1 from homebrew
python-build: use readline from homebrew
Downloading Python-3.8.2.tar.xz...
-> https://www.python.org/ftp/python/3.8.2/Python-3.8.2.tar.xz
Installing Python-3.8.2...
python-build: use readline from homebrew
python-build: use zlib from xcode sdk
Installed Python-3.8.2 to /Users/kk/.pyenv/versions/3.8.2
% which python 
/Users/kk/.pyenv/shims/python

Pythonは.pyenvの下にあるのだが、呼び出そうとするとシステムのPythonが呼び出される。

% python3
Python 3.7.7 (default, Mar 10 2020, 15:43:33) 
[Clang 11.0.0 (clang-1100.0.33.17)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> quit()

pyenvのglobalで3.8.2を宣言。

% pyenv global 3.8.2
% python3
Python 3.8.2 (default, Apr 27 2020, 09:48:02) 
[Clang 11.0.0 (clang-1100.0.33.17)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> quit()

これでPythonを呼ぶと3.8.2が立ち上がるようになる。

% python3
Python 3.8.2 (default, Apr 27 2020, 09:48:02) 
[Clang 11.0.0 (clang-1100.0.33.17)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import numpy
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ModuleNotFoundError: No module named 'numpy'
>>> quit()

numpyをimportしてみるがもちろん入っていない。 以降はpipで欲しいライブラリを入れていく。

まずはpipのアップデート。

% which pip3
/Users/kk/.pyenv/shims/pip3
% pip3 install --upgrade pip
Collecting pip
  Using cached https://files.pythonhosted.org/packages/54/0c/d01aa759fdc501a58f431eb594a17495f15b88da142ce14b5845662c13f3/pip-20.0.2-py2.py3-none-any.whl
Installing collected packages: pip
  Found existing installation: pip 19.2.3
    Uninstalling pip-19.2.3:
      Successfully uninstalled pip-19.2.3
Successfully installed pip-20.0.2

pip のアップデート完了後、自分が欲しいライブラリをインストール。私の場合、以下のもの。

pip3 install numpy        # 数値計算用ライブラリ
pip3 install scipy        # 数値解析用ライブラリ
pip3 install matplotlib   # グラフ作成用ライブラリ
pip3 install pandas       # データ加工支援用ライブラリ
pip3 install pillow       # 画像処理用ライブラリ
pip3 install openpyxl     # エクセルデータ読み書き
pip3 install xlrd         # エクセルデータ読込用
pip3 install xlwt         # エクセルデータ書込用
pip3 install scikit-learn # 機械学習ライブラリ
pip3 install seaborn      # グラフ作成ライブラリ
pip3 install sympy        # 記号計算用ライブラリ
pip3 install jupyter      # Jupyter

上記で入れたものを確認。依存関係にあるものは自動的にインストールされている模様。

% pip3 list
Package            Version     
------------------ ------------
appnope            0.1.0       
attrs              19.3.0      
backcall           0.1.0       
bleach             3.1.4       
cycler             0.10.0      
decorator          4.4.2       
defusedxml         0.6.0       
entrypoints        0.3         
et-xmlfile         1.0.1       
ipykernel          5.2.1       
ipython            7.13.0      
ipython-genutils   0.2.0       
ipywidgets         7.5.1       
jdcal              1.4.1       
jedi               0.17.0      
Jinja2             2.11.2      
joblib             0.14.1      
jsonschema         3.2.0       
jupyter            1.0.0       
jupyter-client     6.1.3       
jupyter-console    6.1.0       
jupyter-core       4.6.3       
kiwisolver         1.2.0       
MarkupSafe         1.1.1       
matplotlib         3.2.1       
mistune            0.8.4       
mpmath             1.1.0       
nbconvert          5.6.1       
nbformat           5.0.6       
notebook           6.0.3       
numpy              1.18.3      
openpyxl           3.0.3       
pandas             1.0.3       
pandocfilters      1.4.2       
parso              0.7.0       
pexpect            4.8.0       
pickleshare        0.7.5       
Pillow             7.1.2       
pip                20.0.2      
prometheus-client  0.7.1       
prompt-toolkit     3.0.5       
ptyprocess         0.6.0       
Pygments           2.6.1       
pyparsing          2.4.7       
pyrsistent         0.16.0      
python-dateutil    2.8.1       
pytz               2019.3      
pyzmq              19.0.0      
qtconsole          4.7.3       
QtPy               1.9.0       
scikit-learn       0.22.2.post1
scipy              1.4.1       
seaborn            0.10.1      
Send2Trash         1.5.0       
setuptools         41.2.0      
six                1.14.0      
sympy              1.5.1       
terminado          0.8.3       
testpath           0.4.4       
tornado            6.0.4       
traitlets          4.3.3       
wcwidth            0.1.9       
webencodings       0.5.1       
widgetsnbextension 3.5.1       
xlrd               1.2.0       
xlwt               1.3.0       
% 

これにて作業完了。

以 上

流況曲線作図(改訂版)

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

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

以 上

VS code と BasicTeX 2020 をインストール

2020年4月12日、最近話題のVS codeをインストール。また、BasicTeX 2020もインストール。

BasicTeXは、2019がインストールしてあったので、まずこれを削除。

また、which ghostscript しても、なぜだか「見つからない」と言われるので、再インストール。

BasicTeX のインストールは、ここに従ってインストール。

具体的には以下のとおり。

brew cask install basictex

/Library/TeX/texbin にパスを通す。 パスの確認は echo $PATH でできる。 (BasicTeX 2019のためにパスを通してあったので、実際にはいじらない)

sudo tlmgr update --self --all
sudo tlmgr install latexmk
sudo tlmgr install collection-langjapanese
sudo tlmgr install courier
sudo tlmgr install newtx txfonts helvetic fontaxes boondox
sudo tlmgr install kastrup tex-gyre

mktexlsr は自動でやってくれる。

TeX文書のプレアンブルは以下の通り。ただし英語用。

\documentclass[10pt]{article}
\usepackage[a4paper,top=25mm,bottom=20mm,left=25mm,right=25mm]{geometry}
\usepackage[T1]{fontenc}
\usepackage{newtxtext,newtxmath,bm}
%\usepackage{amsmath,amssymb,bm}
\usepackage{float}
\usepackage{fancyvrb}
%\usepackage[titletoc,title]{appendix}
\usepackage[dvipdfmx]{graphics,graphicx}

\let\olditemize\itemize
\renewcommand{\itemize}{
\olditemize
\setlength{\itemsep}{1.2pt}
\setlength{\parskip}{0pt}
\setlength{\parsep}{0pt}
}

\begin{document}

......

\end{document}

\usepackage{newtxtext,newtxmath,bm} か、 \usepackage{amsmath,amssymb,bm} のどちらかを入れないと数式を認識してくれないので注意。

VS codeは少し触ってみたが、評判どおりとても軽快である。

以 上

matplotlibで使える色一覧

色一覧

f:id:damyarou:20200402074423p:plain:w700

画像作成プログラム

import matplotlib
import colorsys
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams['font.family'] = 'Ricty Diminished'


def penc(hval):
    r=int(hval[1:3],16)
    g=int(hval[3:5],16)
    b=int(hval[6:8],16)
    hsv=colorsys.rgb_to_hsv(r/255, g/255, b/255)
    h=hsv[0]
    s=hsv[1]
    v=hsv[2]
    col='#000000'
    if 0.5<h: col='#ffffff'
    if v<0.75: col='#ffffff'
    return col
    

cdic=matplotlib.colors.cnames
lkey=[]
lval=[]
for key, val in zip(cdic.keys(), cdic.values()):
    lkey=lkey+[key]
    lval=lval+[val]
    
xmin=0
xmax=4
k=-1
fnameF='fig_col_mpl.png'
ymin=0
ymax=37
plt.figure(figsize=(6,12))
plt.xlim([xmin,xmax])
plt.ylim([ymax,ymin])
plt.axis('off')
fsize=8 # fontsize
k=-1
for i in range(int(ymin),int(ymax)):
    for j in range(0,4):
        k=k+1
        if k<len(lval):
            xs=float(j)
            xe=xs+1.0
            ys=float(i)
            ye=ys+1.0
            xx=[xs,xe,xe,xs]
            yy=[ys,ys,ye,ye]
            plt.fill(xx,yy,color=lval[k])
            text1=lkey[k]
            text2=lval[k]
            xg=0.5*(xs+xe)
            yg=0.5*(ys+ye)
            col=penc(lval[k])
            plt.text(xg,yg-0.25,text1,rotation=0,ha='center',va='center',fontsize=fsize,color=col)
            plt.text(xg,yg+0.25,text2,rotation=0,ha='center',va='center',fontsize=fsize,color=col)
plt.savefig(fnameF, dpi=300, bbox_inches="tight", pad_inches=0)

以 上