WANtaroHP (Python)

Contents

Python を始める (12 October 2014)

きっかけ

2014年9月27日より Windows版Python (64bit) をいじり始めました.きっかけは,丁度そのころ日本に一時帰国していたのですが,何か新しいことを始めよう,と思い立ちPythonをいじり始めることにしました. なぜPythonかというと,明確な理由はありませんが,ネットの情報によれば,

技術計算ができるらしい
作図もできるらしい
画像処理もできるらしい
GUI もできるらしい

ということでどんなものだかやってみたくなったのがきっかけです.

インストール

以下のファイルあるいはインストーラをダウンロードし,インストールします.

機能ProgramDownload site
Python本体python-3.4.1.amd64.msihttps://www.python.org/downloads/
技術計算用ライブラリnumpy-MKL-1.9.0.win-amd64-py3.4.exehttp://www.lfd.uci.edu/~gohlke/
技術計算用ライブラリscipy-0.14.0.win-amd-py3.4.exehttp://www.lfd.uci.edu/~gohlke/
グラフ描画ライブラリmatplotlib-1.4.0.win-amd64-py3.4.exehttp://www.lfd.uci.edu/~gohlke/
matplotlib用setuptools-5.8.win-amd64-py3.4.exehttp://www.lfd.uci.edu/~gohlke/
matplotlib用pytz-2014.7.win-amd64-py3.4.exehttp://www.lfd.uci.edu/~gohlke/
matplotlib用pyparsing-2.0.2.win-amd64-py3.4.exehttp://www.lfd.uci.edu/~gohlke/
matplotlib用python-dateutil-2.2.win-amd64-py3.4.exehttp://www.lfd.uci.edu/~gohlke/
matplotlib用six-1.8.0.win-amd64-py3.4.exehttp://www.lfd.uci.edu/~gohlke/
画像処理ライブラリPillow-2.6.0.win-amd64-py3.4.exehttp://www.lfd.uci.edu/~gohlke/
matplotlib用basemapbasemap‑1.0.8‑cp34‑none‑win_amd64.whlhttp://www.lfd.uci.edu/~gohlke/pythonlibs/

勉強用サイト

主に以下のサイトの情報で勉強しています.

サイト URL
http://matplotlib.org/gallery.html
http://turbare.net/transl/scipy-lecture-notes/intro/matplotlib/matplotlib.html
http://bicycle1885.hatenablog.com/entry/2014/02/14/023734
http://kaiseki-web.lhd.nifs.ac.jp/wiki/index.php/Matplotlib_サンプル集#ASCII.E3.83.87.E3.83.BC.E3.82.BF.E3.81.AE.E8.AA.AD.E3.81.BF.E8.BE.BC.E3.81.BF.E3.83.BB.E6.8F.8F.E7.94.BB
http://guutaranoheya.web.fc2.com/math/matplotlib_ex.html
http://librabuch.jp/2013/05/python_pillow_pil/

メモ(03 December 2014)

  • 数学関数:スカラーは math の数学関数で,numpy.array は numpy の数学関数で処理するのがベター.ものによってはエラーがでる.
  • プロットは list か array か:matplotlibでのプロットは list より numpy.array がベター.ものによってはlistだとエラーが出る.
  • グローバル変数:グローバル変数は,それを使う関数に中で,global の宣言を行う.
  • 複数戻り値:関数において,入力は引数,出力は return で.return は複数の戻り値を返せるため,これを有効に使う.
  • 半透明色重ね合わせ:同色系の濃淡は,alpha を指定し,半透明色の重ねあわせを使う手もある.結構綺麗.
  • 名前の重複を避ける:import numpy as np とした場合,np を変数として使わないよう注意.ポイント数などで np を使いたくなるがこれはダメ.
  • 最大値のインデックス:list では n=Q.index(max(Q)),array では n=numpy.argmax(P) で取得.

basemapのインストール:拡張子 whl のインストール (25 Jan 2015)

matplotlibで地図を出力するためのライブラリをインストール. ここ(http://www.lfd.uci.edu/~gohlke/pythonlibs/)からこれ(basemap‑1.0.8‑cp34‑none‑win_amd64.whl)をとってくるまではよかったが, whlファイルをどうやってインストールするかさっぱりわからず. 色々調べたが結局pipなるものを使うことが分かり,Python のディレクトリ(私の場合はc:\Python34の下)のディレクトリ\Scriptsにある, pip3.4.exeを使って,「c:\Python34\Scripts\pip3.4 install basemap-1.0.8-cp34-none-win_amd64.whl」によりインストール完了. Python,複雑!

こちら(http://matplotlib.org/basemap/)のサイトに色々な図法による描画事例があります.

こちら(http://matplotlib.org/basemap/api/basemap_api.html)にオプションの解説があります.

現在時刻を取得する (12 October 2014)

py_dtime.py ソースコード

# coding: utf-8
import datetime

mdic={
1:'January',\
2:'February',\
3:'March',\
4:'April',\
5:'May',\
6:'June',\
7:'July',\
8:'August',\
9:'September',\
10:'October',\
11:'November',\
12:'December',}
d = datetime.datetime.today()
print(d)
print('%02d %s %4d at %02d:%02d:%02d\n'\
% (d.day,mdic[d.month],d.year,d.hour,d.minute,d.second))

出力

2014-10-12 17:22:50.682100
12 October 2014 at 17:22:50

tkinter を用いてマンデルブロ集合を描く (12 October 2014)

概要

Python本体付属のtkinterでの作図を試みました. GUI なのでイベントドリブンです.また自分で作図用の座標を定義しないとならないのでコードが煩雑になります. プログラム最後の mainloop() がないと何も起こりません. しかし,いわゆるグラフではなくこのような描画には向いているのかもしれません.

py_mandel.py ソースコード

# coding: utf-8
from tkinter import *
import sys

def COORDX(xx,kxi,kxf,xmin,xmax):
    return kxi+int((xx-xmin)*(kxf-kxi)/(xmax-xmin))
def COORDY(yy,kyi,kyf,ymin,ymax):
    return kyf-int((yy-ymin)*(kyf-kyi)/(ymax-ymin))

def pr():
    d=c.postscript(file="fig_mandel.eps")
def quit():
    sys.exit()

f=Frame()
b1=Button(f,text='write canvas to fig_mandel.eps', command=pr)
b1.pack(side='left',expand=YES,fill='x')
b2=Button(f,text='quit',command=quit)
b2.pack(side='right',expand=YES,fill='x')
f.pack(fill='x')

ks=600
c=Canvas(width=ks,height=ks,background='white')
c.pack()

kxi=0
kyi=0
kxf=ks
kyf=ks
xmin=-2.0
xmax= 2.0
ymin=-2.0
ymax= 2.0
dx=(xmax-xmin)/float(ks)
dy=(ymax-ymin)/float(ks)
col_list=['red','yellow','lime','cyan','blue','magenta']
kmax=100
for i in range(0,ks):
    for j in range(0,ks):
        ar=xmin+dx*float(i)
        ai=ymin+dy*float(j)
        zr=0.0
        zi=0.0
        k=0
        while k<kmax:
            k=k+1
            rp=zr*zr-zi*zi+ar
            ip=2.0*zr*zi+ai
            if rp*rp+ip*ip>4.0:
                icol=k%6
                kxx=COORDX(ar,kxi,kxf,xmin,xmax)
                kyy=COORDY(ai,kyi,kyf,ymin,ymax)
                c.create_rectangle(kxx,kyy,kxx+1,kyy+1,\
                fill=col_list[icol],outline=col_list[icol])
                break
            zr=rp
            zi=ip

mainloop()

描画事例

pic
表示された画像のスクリーンショット

Pillow を用いてデジカメの画像を縮小する (12 October 2014)

概要

簡単な画像を扱うサンプルとして,デジカメの画像(約3MB)を250kB程度に縮小するコードを作って見ました. ここでは,画像の縦横比を保ちながら幅を600ピクセルに縮小しています.また私のデジカメの画像は拡張子が大文字の「.JPG」で作成されるので,縮小画像名ではこれを小文字に変えています. ついでに縮小画像を html で表示する部分も付け加えてあります.

facebookなどへの投稿の際,あまり大きな画像を送るのも気がひけるので,これまでは Visual C# で作ったプログラムを用いていましたが,Visual シリーズとはお別れしたいので,今回は Python で作って見たわけです. また縮小画像を html で表示することにより,アップしたい画像を比較的簡単に探し出すことができます. 使い方は,元画像のあるフォルダにプログラム py_html.py をコピーし,実行するだけです. 現時点では GUI の利用までは行っておらず,表示はコーディングが簡単な html 文書内としました.

py_html.py ソースコード

# coding: utf-8
from PIL import Image
import os

# Resize images
files = os.listdir()
for file in files:
    if 0<file.find('.JPG'):
        width=600
        imgname=str(width)+'_'+file.replace('.JPG','.jpg')
        img = Image.open(file, 'r')
        height=int(float(width)/float(img.size[0])*float(img.size[1]))
        img.thumbnail((width, height), Image.ANTIALIAS)
        img.save(imgname, 'JPEG', quality=100, optimize=True)
        print(file,imgname)

# Write html document
files = os.listdir()
imgn=[]
i=-1
for file in files:
    i=i+1
    if 0==file.find('600_') and 0<file.find('.jpg'):
        imgn=imgn+[file]
fnameW=open('small.html',"w")
fnameW.write('<html>\n')
fnameW.write('<body>\n')
fnameW.write('<table>\n')
col=4
width=200
k=-1
for i in range(0,int(len(imgn)/col)+1):
    fnameW.write('<tr>\n')
    for j in range(0,col):
       k=k+1
       if len(imgn)-1<k:
           fnameW.write('<td></td>\n')
       else:
           fnameW.write('<td align="center"><img src="%s" alt="pic" width="%d"><br>%s</td>\n' % (imgn[k],width,imgn[k]))
    fnameW.write('</tr>\n')
fnameW.write('</table>\n')
fnameW.write('</body>\n')
fnameW.write('</html>\n')
fnameW.close()

描画事例

html内では,tableタグにより画像を配置しています.ここではtableは4列とし,行数は画像の数により調整するようにしています. 画像縮小は,resize を使う方法と thumbnail を使う方法がありますが,ネットで調べていたところ thumbnail のほうが縮小画像が綺麗という情報であり,プログラムでは thumbnail を用いています. (実際に比較して見ましたがそのとおり.thumbnailのほうが断然綺麗です)

pic
表示されたhtml文書のスクリーンショット

matplotlib によるグラフ作成 (12 October 2014)

いくつかの描画事例を示します.あくまでもコーディングの事例なのでデータは添付していません. 結構かっこいいグラフが描けそうです.また TeX の書式で数式を直接書き込めるのも便利です.

ソース

FilenameDescription
py_fig_acc.txt地震加速度描画(3成分)
py_fig_spc.txt加速度応答スペクトル描画(3成分)
py_fig_hist.txtヒストグラム描画(乱数モデルデータ)
py_fig_time_his_TR.txt前橋市の月平均気温と月間降水量の推移

描画グラフ

pic
地震加速度描画(3成分)
pic
加速度応答スペクトル描画(3成分)
pic
ヒストグラム描画(乱数モデルデータ)
pic
前橋市の月平均気温と月間降水量の推移

画像縮小とhtmlによる表示 (22 October 2014)

概要

「Pillow を用いてデジカメの画像を縮小する (12 October 2014)」から少し発展させて, 「メニューでダイアログを開きフォルダを指定してフォルダ内のデジカメ画像ファイルを縮小,html で表示する」 というプログラムにしました.

Python の tk で縮小画像を表示する方法も試しましたが,結局スクロールバーの表示やクリックしたら拡大画像を表示する機能の実装がうまくいかず, これらを簡単に実現できる html での表示になっています.

Image のインポート「from PIL import Image」はtkinterのimportのあとに書かないと,何かしら衝突しているのか, 「AttributeError: type object 'Image' has no attribute 'open'」というエラーが出てしまうので注意.

簡単な解説

  • メニューによりデジカメ画像が保存されているフォルダを選定
  • メニューにより画像縮小とhtml作成・表示を実行
  • 連続処理する場合を考慮し,1回のプログラム起動で処理したフォルダ名をリストボックスに表示
  • デジカメ画像の拡張子は大文字「JPG」
  • 縮小後の画像サイズは横600ピクセルで固定
  • 縮小画像は横600ピクセルのもの1種類とし,html表示では横200に縮小,ファイル名クリックにより横600ピクセルの画像を表示する
  • 縮小画像および表示用htmlは,元画像と同じフォルダに保存される
  • 縮小画像ファイル名は元画像ファイル名の先頭に「600_」がついている
  • html文書のファイル名は「_small.html」に固定
メニュー 機能
Open フォルダ選択.もし選択したフォルダ内に,ファイル「_small.html」が存在すれば画像をブラウザで表示
Run & Show選択したフォルダ内の全JPG画像(拡張子大文字)を縮小しブラウザで表示
Exit プログラム終了

ソース

FilenameDescription
py_html_menu.txt画像縮小・html表示プログラム

二分法 (23 October 2014)

Pythonを始めて比較的初期に作ったプログラムであるが,アップしていなかったので載せておく. この程度の規模のプログラムであれば何で作っても速度・精度とも大差はないようである.

C や Fortran と比べ変数宣言が無い分,コードはあっさりしている.主要部分はほとんど同じコーディング量である. awk と比べると,awk は C の変数宣言が無いプログラムというイメージであるが,Python は独特の感じである.

py_bisection.py ソースコード

# coding: utf-8
from math import *

def FUNC(x):
    return cos(x)*cosh(x)+1.0

eps=1.0e-10
m=10
n=100
print("{0:>5} {1:>5} {2:^15} {3:^15}"\
.format('k','i','xi','fi'))
for k in range(0,m):
    x1=float(k)*pi
    x2=x1+pi
    for i in range(0,n):
        xi=0.5*(x1+x2)
        f1=FUNC(x1)
        f2=FUNC(x2)
        fi=FUNC(xi)
        if abs(x2-x1)<eps:
            break
        if f1*fi<0.0:
            x2=xi
        if fi*f2<0.0:
            x1=xi
    print("{0:5d} {1:5d} {2:15.7e} {3:15.7e}"\
    .format(k,i,xi,fi))

出力

    k     i       xi              fi       
    0    35   1.8751041e+00  -3.1298519e-11
    1    35   4.6940911e+00   1.0887963e-09
    2    35   7.8547574e+00   2.1079174e-08
    3    35   1.0995541e+01  -6.8795489e-07
    4    35   1.4137168e+01   2.3681529e-06
    5    35   1.7278760e+01  -3.3951512e-04
    6    35   2.0420352e+01   3.8182270e-03
    7    35   2.3561945e+01  -1.7217838e-01
    8    35   2.6703538e+01  -8.0416410e+00
    9    35   2.9845130e+01  -2.0823992e+02

色見本 (09 November 2014)

最近 html で作成している色見本をよく見るので,Python でも同様のものを作って見ました.

やっていることは,以下の2つ.

(1) 前半で16進での色コードを作成.
(2) 後半では作図.
    x 軸を 0 から 6,y 軸を 0 から 36 として,「長さ 1 × 高さ 1」の矩形を描きながら,その中心にコードを描画していく.
    矩形描画にはRectangleを使用.「左下座標,長さ,高さ」を指定して矩形を形成する.

py_misc_color.py ソースコード

from matplotlib import rcParams
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Consolas']
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle

h_list=['00','33','66','99','cc','ff']

scolor=[['']*6 for i in range(0,36)]
ii=-1
for c1 in h_list:
    for c2 in h_list:
        ii=ii+1
        sh=[]
        for c3 in h_list:
            sh=sh+['#'+c1+c3+c2]
        scolor[ii][:]=sh
        del sh

sx=0.02
sy=0.1
fig = plt.figure(figsize=(6,9))
#ax = fig.add_subplot(111)
ax = plt.gca()
for i in range(0,36):
    for j in range(0,6):
        col=scolor[35-i][j]
        x=float(j)+sx
        y=float(i)+sy
        dx=1.0-2.0*sx
        dy=1.0-2.0*sy
        xg=x+0.5*dx
        yg=y+0.5*dy
        cn='black'
        if j<4:
            cn='white'
        ax.add_patch(Rectangle((x,y), dx, dy, color=col))
        plt.text(xg,yg,col,ha='center',va='center',fontsize=9,color=cn)
plt.xlim([0,6])
plt.ylim([0,36])
plt.axis('off')
plt.savefig("fig_color.png", dpi=200)
plt.show()

上のコードではRectangleで矩形を描画しているが,Polygonを用いるなら以下のようにすることができる. 座標として,矩形を構成する4点を与える.

........
from matplotlib.patches import Polygon
........
p1=[x1,y1]
p2=[x2,y2]
p3=[x3,y3]
p4=[x4,y4]
ax.add_patch(Polygon([p1,p2,p3,p4], closed=True,fill=True,color=col)

出力

pic
色見本

x-yプロットの事例 (09 November 2014)

x-yプロットの事例を 3 つ紹介します. Python プログラムは,いずれも作図のみのプログラムで,数値計算は Fortran でやっています. 結構かっこいいグラフが描けた?

プログラム作成には,以下のサイトを参考にさせていただいています.

aspect ratio : http://d.hatena.ne.jp/bettamodoki/20120605/1338863706
custom dashes: http://stackoverflow.com/questions/14710221/python-matplotlib-dash-dot-dot-how-to
font : http://matplotlib.org/examples/api/font_family_rc.html

いくつかのテクニック

フォントの指定

png 出力では特に感じませんでしたが,eps 出力する場合,フォントを指定しておかないと serif 系のショボいフォントになってしまうので, Arial を明確に指定しています.

from matplotlib import rcParams
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['arial']
import numpy as np
import matplotlib.pyplot as pl

読み込みたいデータ列の指定

Numpy の loadtxt において,usecols で読み込む列を指定しています. ここで用いているデータは,空白区切り.usecols=(3,0) のように,順番を入れ替えることもできるようです.

data1 = np.loadtxt('surface1.txt', usecols=(0,3))
data2 = np.loadtxt('surface2.txt', usecols=(0,3))
data3 = np.loadtxt('surface3.txt', usecols=(0,3))

グラフの縦横比の指定

グラフの縦横比をaspectで指定しています. 下の事例では,縦軸長さを横軸の半分にしています.

aspect = 0.5*(ax.get_xlim()[1] - ax.get_xlim()[0]) / (ax.get_ylim()[1] - ax.get_ylim()[0])
ax.set_aspect(aspect)

グラフ内に強調区間を入れる

下の事例では,x軸(横軸)で 0 から 83 の範囲を灰色に,98 から 145 の範囲を薄い黄色で塗りつぶしています. ラベルもつけることができ,凡例に反映されます.

pl.axvspan(0, 83, facecolor='#d3d3d3', edgecolor='#d3d3d3',alpha=0.5,label='Dam')
pl.axvspan(98, 145, facecolor='#ffffbb',edgecolor='#ffffbb', alpha=0.5,label='Power Station')

線種の指定

線種は '-.' などにより簡単に設定できますが,あまり綺麗ではないので,set_dashes を用いて設定しています. 偶数個の数値で指定する必要があるようであり,set_dashes([5,2]) なら,「5 だけ実線,2 だけ空白」という形で線種を指定します.

line1=pl.plot(data1[:,0],data1[:,1]*1000,color='red'  ,linewidth=2.0,label='D=1000MPa')
line2=pl.plot(data2[:,0],data2[:,1]*1000,color='blue' ,linewidth=2.0,label='D=2000MPa')
line3=pl.plot(data3[:,0],data3[:,1]*1000,color='green',linewidth=2.0,label='D=3000MPa')
line1[0].set_dashes([5, 2])
line2[0].set_dashes([8,4,2,4])
line3[0].set_dashes([8, 4, 2, 4, 2, 4])

凡例の線の長さ

pl.legendで凡例を描画しますが,handlelengthで凡例の線の長さを指定します. 短すぎると一点鎖線や二点鎖線うぃ表示する場合中途半端になってしまうので,そのような時に使用できます.

pl.legend(shadow=True, prop={'size':12},loc='lower right',handlelength=3)

ソース

FilenameDescription
py_plot_dis.txt基礎地盤変位分布作図
py_plot_mn.txtRC断面スラブ耐荷力作図(軸力-モーメント関係図)
py_plot_m0.txtRC断面スラブ耐荷力作図(鉄筋間隔と抵抗モーメント)

描画グラフ

pic
基礎地盤変位分布
pic
RC断面スラブの耐荷力曲線(軸力-モーメント関係図)
pic
RC断面スラブの耐荷力曲線(鉄筋間隔と抵抗モーメントの関係)

主応力ベクトル図,主応力分布図 (09 November 2014)

2次元FEM応力解析結果から,主応力ベクトル図および主応力分布図(主応力の大きさによる要素の色分け)作成のためのPythonスクリプトと実行用バッチファイルを掲載しています. 図の作成はGMTにて行っています. また4枚の図はそれぞれ独自に作成していますが,ImageMagickにより結合させ,1枚のpng画像にしています.

Pythonでこの図まで作りたいのですが,手間取っており,結局はPythonはFEM出力からのデータ読み込みとGMT作図用データ出力までとなっています.

GMTによる処理の概要

主応力ベクトル図

要素重心座標,主応力の大きさ,方向(角度)を読み込み,色分けする数だけファイルを作成し,これをGMTで読み込ませています. ベクトルは小さいので,矢印描画オプションを用いますが,矢を省略した線分のみの表示にしています. GMTの場合,polygonに着色する場合,1色1コマンドですので,色の数だけファイルを作成して処理します.

主応力の大きさ分布図

1色1コマンドの考え方で,主応力の大きさから何色にするかを分類し,同じ色で構成される要素の座標(4点づつ)をファイルに書き込みます. GMTの中では,読み込んだ座標で囲まれる範囲を指定された色で塗りつぶす操作を行います.

Pythonによる処理の概要

主応力ベクトル図

未実施

主応力の大きさ分布図

挑戦中.

1要素4節点をPolygonで定義し,主応力の大きさに応じてそのPolygonに色を与えます. GMTのようにいくつもファイルを作成して読み込む必要がなく,手続きは比較的簡単. しかしながら,座標軸の設定や凡例をどのように表示するかを思案中. 要素の色分けまでは特に問題なくできている. ただしGMTと比べると遅い. また大きな図の中に小さな図を挿入する場合,GMTでは元の大きなグラフの描画範囲を何cmx何cmという形で管理しており, この範囲(外への拡大も可能)で完全に重ねあわせが可能であり,まさに実際の紙に切り張りするイメージで埋め込み(重ね合わせ)ができるのに対し, Python の matplotlib では全体の座標管理のイメージがまだつかめておらず,net 検索した事例ではうまくできていても実際にやろうとするとうまくいかないという状態が続いている.

現在の結論

  • 通常の 2 次元グラフなら Python のほうが凡例の設定が楽であり GMT より便利である( GMT も凡例作成のコマンドはあるのだが,使い方がさっぱり分からず,自分で gawk でスクリプトを組み凡例出力を行っている). またデータ作成から作図まで1つのプログラムで完結できるところが良い.
  • Fortran プログラムの出力リストからのデータ読み込みやその加工は,今まで使っていた gawk より自由度も高く十分使い物になる. gawk は大きな配列をとると暴走することがあったので,FEM出力など,大量データの処理は別途 Fortran で作ったプログラムで行うこともあったが,この部分は Python があればわざわざコンパイラを起動してアーダコーダやらずに済ますことができそう.
  • 2 次元座標中の Polygon の色分けは 3 次元プロットにも繋がるものであり,ここらへんは今の私の理解度では GMT のほうが格段に使いやすい.

ソース

FilenameDescription
py_fem_resul.txt作図用データ作成
a_resul.txtGMT作図制御用バッチファイル
a_draw.txtGMT作図実行用バッチファイル

描画グラフ

pic
GMT作図:主応力分布図(方向および大きさ)

FFT - numpy & scipy (10 November 2014)

numpy と scipy の fft および ifft を使って見ました.やっていることは以下のとおり.複素数計算になっていることに注意.

  • 波形データ (x) の作成(データ数: nd)
  • FFT にかけるデータ数が 2 の累乗になるようデータ数を調整.ここでは元データ数 10 に対し FFT 用データ数を 16 とし,元データの後ろに 0 を配置している
  • 関数 fft 呼び出しデータ数を調整したデータ (xx) を入力.この時点では実数波形データのみを入力
  • 関数 fft の出力は複素数.出力されたベクトルをデータ数 (nn) で割ることにより複素フーリエ係数 (sp) を求める
  • 逆変換にかけるためには複素数 (sp) をデータ数倍した値 (sp * nn) を ifft に入力
  • 逆変換の出力として得られた複素数 (wv) は,誤差の範囲で入力データと同じ実数部 (wv.real) が得られる.また当然虚数部 (wv.imag) は 0 となる

py_fft.py ソースコード

import numpy as np
from scipy import fftpack

def FFT_N(xx):
    sp=np.fft.fft(xx)/nn
    spa=np.sqrt(sp.real**2+sp.imag**2)
    wv=np.fft.ifft(sp*nn)
    return sp,spa,wv

def FFT_S(xx):
    sp=fftpack.fft(xx)/nn
    spa=np.sqrt(sp.real**2+sp.imag**2)
    wv=fftpack.ifft(sp*nn)
    return sp,spa,wv


#x=np.array([0.998, 0.567, 0.966, 0.748, 0.367, 0.481, 0.074, 0.005,
#            0.347, 0.342, 0.218, 0.133, 0.901, 0.387, 0.445, 0.662])
x=np.random.random_sample(10)
nd=len(x)
k=0
nn=0
while nn<nd:
    k=k+1
    nn=2**k
xx=np.zeros(nn)
xx[0:nd]=xx[0:nd]+x[0:nd]
sp1,spa1,wv1=FFT_N(xx)
sp2,spa2,wv2=FFT_S(xx)

print('%7s %7s %7s %7s %7s %7s'\
       %('x','Re(Ck)','Im(Ck)','Abs(Ck)','Re(x)','Re(x)'))
for i in range(0,nn):
    print('%7.3f %7.3f %7.3f %7.3f %7.3f %7.3f'\
           %(xx[i],sp1.real[i],sp1.imag[i],spa1[i],wv1.real[i],wv1.imag[i]))

print('%7s %7s %7s %7s %7s %7s'\
       %('x','Re(Ck)','Im(Ck)','Abs(Ck)','Re(x)','Re(x)'))
for i in range(0,nn):
    print('%7.3f %7.3f %7.3f %7.3f %7.3f %7.3f'\
           %(xx[i],sp2.real[i],sp2.imag[i],spa2[i],wv2.real[i],wv2.imag[i]))

出力

      x  Re(Ck)  Im(Ck) Abs(Ck)   Re(x)   Re(x)
  0.616   0.327   0.000   0.327   0.616  -0.000
  0.683  -0.001  -0.139   0.139   0.683   0.000
  0.445   0.072  -0.062   0.095   0.445   0.000
  0.878  -0.029  -0.015   0.033   0.878   0.000
  0.166   0.040  -0.030   0.050   0.166   0.000
  0.667   0.026   0.011   0.028   0.667  -0.000
  0.083   0.032  -0.017   0.037   0.083   0.000
  0.653   0.062  -0.072   0.095   0.653   0.000
  0.387  -0.115   0.000   0.115   0.387   0.000
  0.659   0.062   0.072   0.095   0.659   0.000
  0.000   0.032   0.017   0.037  -0.000   0.000
  0.000   0.026  -0.011   0.028   0.000   0.000
  0.000   0.040   0.030   0.050  -0.000  -0.000
  0.000  -0.029   0.015   0.033   0.000  -0.000
  0.000   0.072   0.062   0.095   0.000  -0.000
  0.000  -0.001   0.139   0.139   0.000   0.000

重回帰分析 - np.linalg.solve(a,b) の利用 (15 November 2014)

概要

連立方程式の解法の使用例として,重回帰分析のプログラムを作って見ました. 計算のメイン部分は以下の5行です.[xd]は入力データが格納された行列, ベクトル{c}が求める偏回帰係数, {ye}は入力データに偏回帰係数を乗じて計算した回帰推定値, rrが重相関係数です.

a=np.dot(xd.T,xd)
b=np.dot(yd,xd)
c=np.linalg.solve(a,b)
ye=np.dot(xd,c)
rr=np.corrcoef(yd,ye)[0][1]

重相関係数を自前で計算するなら,以下のようになるか.上のコードではこれを相関行列を求める関数 corrcoef() 1個で処理しています. ここに{yd}は入力された観測値,{ye}は回帰推定値です.

yem=np.average(ye)
ydm=np.average(yd)
y1=np.sum((yd-ydm)*(ye-yem))
y2=np.sum((yd-ydm)**2)
y3=np.sum((ye-yem)**2)
rr=y1/sqrt(y2*y3)

使い方

python py_sta_MRA.py fnameR fnameW
fnameR入力ファイル名(csv)
fnameW出力ファイル名(csv)

ソース

FilenameDescription
py_sta_MRA.txt重回帰分析プログラム
inp_sta_MRA.txt入力データサンプル
out_sta_MRA.txt出力データサンプル

描画グラフ

pic pic
matplotlib作図:観測値と回帰推定値の関係matplotlib作図:データ間相関係数

主成分分析 - np.linalg.eig(a) の利用 (15 November 2014)

概要

固有値解析の使用例として,主成分分析のプログラムを作って見ました. 計算のメイン部分は以下の3行です.[xd]は入力データが格納された行列, [a]は相関行列,{ev}および[vec]が固有値計算で求まる固有値と固有ベクトルです. [xx]は,入力データに固有ベクトルを乗じたものであり,スコアになります.

a=np.corrcoef(xd.T)     # correlation matrix
ev,vec=np.linalg.eig(a) # Eigenvalue analysis
xx=np.dot(xd,vec)       # score

使い方

python py_sta_PCA.py fnameR fnameW
fnameR入力ファイル名(csv)
fnameW出力ファイル名(csv)

ソース

FilenameDescription
py_sta_PCA.txt主成分分析プログラム
inp_sta_PCA.txt入力データサンプル
out_sta_PCA.txt出力データサンプル

描画グラフ

pic
matplotlib作図:スコア相関図

クラスター分析 - k-means++ (16 November 2014)

概要

行列・ベクトル演算の練習として,クラスター分析のプログラムを作って見ました. クラスタの初期値選定を改良したk-means++法を用いており,距離計算はMahalanobis距離を用いています. 分類結果は,主成分分析結果で用いたスコア相関図において,分類されたグループ別にプロットを色分けして表示しています.

k-means法の手順

(1) 各データにクラスタを割り振る
(2) 割り振ったデータでクラスタの中心を計算する
(3) 各データと各クラスタ中心の距離を計算し,最短距離となるクラスタに各データを割り当てなおす
(4) 上記を繰り返し,分類が変化しなくなったら終了する

k-means++法における初期値の与え方

(1) クラスタ数 k を定め,3個の距離を格納するベクトル {D1},{D2},{DR} を準備します.
(2) 1個目の中心点 c1 を乱数を用いて設定します.
(3) ベクトル {D1} に大きな数値をダミーとしてセットします.
(4) 中心点 c1 から全点までの距離を計算し,それらをベクトル {D2}に格納します.
(5) {D1} と {D2} の要素を比較し,小さいほうの値を {DR}に格納します.
(6) 2個目の中心点 c2 を {DR}の中で最大の距離を有する点として設定します.
(7) {D1}={DR} とおいて,(3) から (6) までのステップをk個の中心点を得るまで繰り返します.

使い方

python py_sta_KMP.py kk mds itmax fnameR fnameW
kk クラスタ数(csv)
mds クラスタ内最小要素数(csv)
itmax 最大繰り返し数(csv)
fnameR入力ファイル名(csv)
fnameW出力ファイル名(csv)

ソース

FilenameDescription
py_sta_KMP.txtクラスター分析プログラム
inp_sta_KMP.txt入力データサンプル
out_sta_KMP.txt出力データサンプル

描画グラフ

pic
matplotlib作図:スコア相関図

主応力ベクトルと簡易応力コンター (21 November 2014)

概要

Python による主応力ベクトルと簡易応力コンター作図に成功したので,プログラムと作図事例を掲載しておきます. 「簡易」応力コンターとしたのは,やっていることが,FEMの要素を応力の大きさに応じて色分けして塗りつぶしただけだからです. 主応力ベクトル図と合わせて眺めると,応力状態をイメージすることができると思います.

これらの作図のためのいくつかのテクニックも掲載しておきます. なお,作図時は,以下のモジュールは import しておく.

from matplotlib import rcParams
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['arial']
import matplotlib.pyplot as plt
from math import *
import sys
import numpy as np

上と右の軸の消去

たかがこれだけのことをやるのに,相当な時間を費やした. モジュールを import し,以下により実現.

from mpl_toolkits.axes_grid.axislines import Subplot
......
fig = plt.figure()
ax1 = Subplot(fig,111)
fig.add_subplot(ax1)
......
ax1.axis["right"].set_visible(False)
ax1.axis["top"].set_visible(False)

任意の方向・長さの直線描画(ベクトル図)

ベクトル図は,矢印は用いず,青線と赤線を用いている. 1つの応力ベクトルを plot を用いて表示している.x と y の順番に注意.

ax1.plot([x1,x2],[y1,y2],color=col[icol],lw=0.5)

任意形状四角形の塗りつぶし(コンター図)

コンター図では,FEMの1要素づつを塗りつぶしているが,これにはPolygonを用いている. モジュールのimportを忘れずに.

from matplotlib.patches import Polygon
......
p1=[x[n1],y[n1]]
p2=[x[n2],y[n2]]
p3=[x[n3],y[n3]]
p4=[x[n4],y[n4]]
ax1.add_patch(Polygon([p1,p2,p3,p4], closed=True,fill=True,color=col[icol],lw=0.1))
ax1.add_patch(Polygon([p1,p2,p3,p4], closed=True,fill=False,color='#aaaaaa',lw=0.1))

tickmarkの消去(コンター図)

コンター図で,凡例は棒グラフを加工して使っているが,tickmarkを消したいため,以下によりそれを実現.

ax2.tick_params(top="off")
ax2.tick_params(bottom="off")
ax2.tick_params(right="off")
ax2.tick_params(left="off")

図の消去

これを行わないと,繰り返し描画を行う場合に,alloc: invalid block: ... というエラーが出る場合がある.

plt.close()

ソース

FilenameDescription
py_fem_vect.txt主応力ベクトル図作成プログラム
py_fem_cont.txt簡易応力コンター図作成プログラム

描画グラフ

pic pic
matplotlib作図:第一主応力ベクトルmatplotlib作図:第二主応力ベクトル
pic pic
matplotlib作図:第一主応力コンターmatplotlib作図:第二主応力コンター

軸対称モデルの主応力ベクトルと簡易応力コンター (23 November 2014)

概要

Python による主応力ベクトルと簡易応力コンター作図の軸対称モデル版です. 作図例には,ImageMagickで結合した画像を掲載しています.画像結合バッチファイルの事例を下に示します.

convert fig_cont1.png -trim -bordercolor white -border 20x20 _1.png
convert fig_cont2.png -trim -bordercolor white -border 20x20 _2.png
convert fig_cont3.png -trim -bordercolor white -border 20x20 _3.png
montage -tile 3x1 -geometry +0+0 -quality 100 _1.png _2.png _3.png fig_cont0.png

ソース

FilenameDescription
py_asnt_vect.txt主応力ベクトル図作成プログラム
py_asnt_cont.txt簡易応力コンター図作成プログラム

描画グラフ

pic
matplotlib作図:主応力ベクトル
pic
matplotlib作図:主応力コンター

Flood Routine Analysis (23 November 2014)

概要

Flood Routine Analysis (洪水時貯水池水位解析) プログラムを Python で作って見ました. 解はFortran版と一致することを確認.Fortranでは配列を先に寸法宣言してあるので,データの読み込み量がそれを超えると回りますが結果がおかしくなるので注意. Python版ではその心配はなしです.

以下に簡単なテクニックを書いておきます.

画面に計算ステップを更新しながら出力する

収束計算を行っているので,プログラムがちゃんと動いているかを監視するため,以下のコードで画面に計算ステップを表示するようにしている. これは,こちら (http://qiita.com/airtoxin/items/bb445529c94d3cd871f3) のページを参考にさせていただきました. ついでに,計算の最後に情報として,計算結果の最大値・最小値なども表示するようにしました.

for i in range(0,nd-1):
......
    sys.stdout.write('\r%d / %d' % (i+1,nd-1))
    sys.stdout.flush()
fout.close()
sys.stdout.write('\n')
sys.stdout.write('Time: %10.3f\n'% max(ti))
sys.stdout.write('EL  : %10.3f %10.3f\n'% (min(pEL),max(pEL)))
sys.stdout.write('Qin : %10.3f %10.3f\n'% (min(q_in),max(q_in)))
sys.stdout.write('Qout: %10.3f %10.3f\n'% (min(pQo),max(pQo)))

グラフに最大値の数値を加える

グラフ内に最大値の数値を記載しています. これには,条件を満足するリスト内の要素の位置を返す関数 index を用いています. ここでは,リスト Qi で最大値が出現する位置を n1 として取得し,この位置に text() により最大値を表示するようにしています. このようなことを少ないコードで実現できるところが Python のいいところなのでしょうね.

n1=Qi.index(max(Qi))
plt.text(ti[n1],Qi[n1],'max:%.0f'%max(Qi),fontsize=10,color='black',ha='left',va='bottom')

ソース

FilenameDescription
py_floodrr.txt解析プログラム
py_fig_floodrr.txt作図プログラム
inp_ut_hv.txt入力データ(1):貯水池H-V
inp_ut_pmf1.txt入力データ(2):流入波形
out_ut_pmf1.txt出力事例

描画グラフ

pic
matplotlib作図:流量・水位経時変化図

Simple Viewer (23 November 2014)

概要

Python による簡単な画像表示プログラムです.

最も単純なViewer

以下の3行のプログラムで,指定した画像ファイルを表示してくれます. ただし,ViewerとしてはWindowsの標準が使われているようであり,私の環境では「Windowsフォトビューアー」というのが立ち上がり画像表示します.

from PIL import Image

img=Image.open('fig_hist.png')
img.show()

少しかっこつけたViewer

以下のプログラムは,メニューにより表示したい画像ファイルを選択できるようにしたものです. 画像表示そのものは,上述と同じように,Windowsの標準ビューアーが使われます. なお,かっこつけついでに,windowにListboxを貼り付け,画像表示の履歴とファイサイズなどの情報を表示するようにしています. こういうことに凝りだすと,コードの行数がどんどん増えてきてしまいます.

from tkinter import *
import tkinter.filedialog
from PIL import Image

def open_callback():
    file_name=tkinter.filedialog.askopenfilename(initialdir='C:/')
    img=Image.open(file_name)
    text='info: '+img.format+' ('+str(img.size[0])+','+str(img.size[1])+') '+img.mode
    listbox.insert(END, file_name)
    listbox.insert(END, text)
    print(file_name)
    print(text)
    img.show()

#Menu
root = Tk()
menu = Menu()
root.config(menu=menu)
file_menu = Menu(menu, tearoff=0)
file_menu.add_command(label='Open', command=open_callback)
file_menu.add_separator()
file_menu.add_command(label='Exit', command=root.destroy)
menu.add_cascade(label='File', menu=file_menu)

#Listbox
listbox = Listbox(root,width=50)
listbox.pack()

mainloop()

3次スプライン補間 - scipy.interpolate (interp1d) の利用 (26 November 2014)

概要

Python による3次スプライン補間の事例プログラムです.

ソース

FilenameDescription
py_eng_spline.txt解析プログラム
inp_eng_spl.txt入力データ事例
out_eng_spl.txt出力データ事例

描画グラフ

pic
matplotlib作図:3次スプライン補間事例

水文統計処理練習 (30 November 2014)

概要

114年の日最大雨量データから,確率密度関数あるいは確率分布関数の母数推定を行うものです. 推定母数からQ-Qプロットを作成するプログラムおよびBootstrap法により確率年雨量を求めるプログラムを掲載しました. 用いている確率分布は,対数正規分布(3変数),対数ピアソンIII型分布,一般化極値分布です. Bootstrap法では対数ピアソンIII型分布を用いています.

いくつかのプログラムのポイントを示しておきます.

math モジュールと numpy の数学関数

数学関数による処理は,以下のように使い分けるといいようです.

スカラー : math の数学関数
numpy.array: numpy の数学関数

モジュール

ガンマ分布・正規分布およびガンマ関数を使うため,scipyを利用しています.

import numpy as np
from math import *
import sys
from scipy.stats import gamma,norm
from scipy.special import gamma as Gamma
import matplotlib.pyplot as plt

非超過確率の作成

def NEP(n,alpha):
    ii=np.arange(1,n+1,1)
    return np.array((ii-alpha)/(float(n+1)-2.0*alpha))

平均値・標準偏差の計算

平均値・標準偏差はnumpyの関数で一発で計算できますが,ここでは定義どおりの計算をしています.

mz=np.sum(z)/float(n)                    # mz=np.average(z)
sz=np.sqrt(np.sum((z-mz)**2)/float(n-1)) # sz=np.std(z,ddof=1)
_s=np.sqrt(np.sum((z-mz)**2)/float(n))   # _s=np.std(z,ddof=0)

標準正規分布の%点計算

pp=NEP(n,alpha)
ww=np.array(norm.ppf(pp,0,1))

ガンマ分布の%点計算

pp=NEP(n,alpha)
ww=np.array(gamma.ppf(pp,bb))

文字を相対位置に表示 (transform=ax.transAxes)

ax.text(0.05,0.95,str_t,fontsize=12,fontweight='bold',color='black',ha='left',va='top',transform=ax.transAxes)

2点を結ぶ直線

ax.plot([xmin,xmax],[ymin,ymax],'-',color='black')

行数不明データの入力 (while 1:)

1行目を読み飛ばし,2行目以降をリストx_orgに格納しています.

x_org=[]
fin=open(fnameR,'r')
text=fin.readline()
while 1:
    text=fin.readline()
    if not text: break
    text=text.strip()
    text=text.split()
    x_org=x_org+[float(text[0])]
fin.close()
n=len(x_org)

ヒストグラム重ね合わせ

alphaを小さくして半透明にし,重ね合わせています.

ax1.hist(data1,bins=ndx,range=(xmin,xmax),alpha=0.2, histtype='stepfilled', color='red',label='100years')
ax1.hist(data2,bins=ndx,range=(xmin,xmax),alpha=0.2, histtype='stepfilled', color='yellow',label='1000years')
ax1.hist(data3,bins=ndx,range=(xmin,xmax),alpha=0.2, histtype='stepfilled', color='blue',label='10000years')

鉛直線を書く (axvline)

ax2.axvline(x=np.median(data1),linestyle='--',color='red')

文字を回転して表示 (rotation)

ax2.text(np.median(data1),0.95,txt1,fontname='monospace',fontsize=12,color='red', ha='right',va='top',rotation=90)

凡例の線の高さ調整 (handleheight)

本来ヒストグラムのための凡例ですので,完全に線にすることはできないようです.

ax2.legend(fontsize=10,title='LP3',loc='lower right',handleheight=0.1)

ソース

FilenameDescription
py_HFA_qq.txtQ-Qプロット作成プログラム
py_HFA_bs.txtBootstrap法プログラム
inp_RF_M.txt入力データ事例(前橋雨量データ)
out_RF_M.txt出力データ事例(Q-Qプロットプログラム出力)

描画グラフ

pic
matplotlib作図:Q-Qプロット(LN3)
pic
matplotlib作図:Bootstrap法による推定値分布図(LP3)

サージング解析 (30 November 2014)

概要

Python によるサージング解析プログラムです.Fortran版を書き換え,matplotlibによる作図を追加しました. 計算時間は1ケース10秒程度以内であり,使い物になりそうです.

ソース

FilenameDescription
py_eng_surge.txt解析・作図プログラム
inp_KN1_Ha.txt入力データ事例(1):導水路調圧水槽負荷遮断
inp_KN1_Hb.txt入力データ事例(2):導水路調圧水槽半負荷急増
inp_KN1_Hc.txt入力データ事例(3):導水路調圧水槽入力遮断
inp_KN1_Ta.txt入力データ事例(4):放水路調圧水槽負荷遮断
inp_KN1_Tb.txt入力データ事例(5):放水路調圧水槽半負荷急増
inp_KN1_Tc.txt入力データ事例(6):放水路調圧水槽入力遮断

描画グラフ

pic pic
pic pic
pic pic
matplotlib作図:導水路調圧水槽水位変化図 matplotlib作図:放水路調圧水槽水位変化図

河川断面等流解析 (02 December 2014)

概要

任意形状の河川断面での等流解析プログラムです. 河川勾配と断面形を入力し,計算する水位はプログラム内で指定するようにしています.

numpyas np で呼び出す場合,変数名に np を使わないよう注意のこと!

等流計算では,水位を与えて相当する流量を求めていますが,下に示すように,水位計算範囲の最小値,最大値,増分を定め, np.arange で array をつくり,これを for ループに適用するのが Python らしいように思います.

ymin=72
ymax=85
dy=0.1
yy=np.arange(ymin,ymax+dy,dy)
for el in yy:
    ....

断面図の作図では,「fill_between」により,指定した水面以下を薄い青で塗るようにして見ました. 下のスクリプトで,xy は地形の x-y座標を格納したリスト(データ入力時にリストを作成している), z1z2 は水面標高を表すリストで,地形を示す座標と同じ要素数のものを準備します. これらを array で numpy の配列 xxyyel1el2 に変換してプロットに用います. また「fill_between」に「interpolate=True」を入れることにより,指定した水位線と地形との交点を探して結んでくれます.

z1=[pmf for i in range(0,nn)]
z2=[ttf for i in range(0,nn)]
xx=np.array(x)
yy=np.array(y)
el1=np.array(z1)
el2=np.array(z2)
xmin=-60.0
xmax=100.0
ymin=70.0
ymax=100
fig = plt.figure()
ax=fig.add_subplot(111)
ax.plot(xx,yy,color='black',lw=1.0,label='ground surface')
ax.fill_between(xx,yy,el1,where=yy<=el1,facecolor='cyan',alpha=0.3,interpolate=True)
ax.fill_between(xx,yy,el2,where=yy<=el2,facecolor='cyan',alpha=0.3,interpolate=True)

断面図作成プログラムの png 出力は,余白が大きいので,以下の ImageMagick コマンドで余白を縮小(一端削除し20ピクセルの余白を追加)しています.

convert fig_lt_HQ.png -trim -bordercolor white -border 20x20 fig_HQ.png
convert fig_sec_XS76.png -trim -bordercolor white -border 20x20 fig_XS76.png
convert fig_sec_XS77.png -trim -bordercolor white -border 20x20 fig_XS77.png

ソース

FilenameDescription
py_eng_uniflow.txt解析プログラム
py_eng_uni_figHQ.txtH-Q作図プログラム
py_eng_uni_figsec.txt解析断面作図プログラム
inp_uni_XS76.txt入力データ事例(1)
inp_uni_XS77.txt入力データ事例(2)

描画グラフ

pic
matplotlib作図:水位-流量関係図
pic
pic
matplotlib作図:解析断面図

マーカーサンプル (11 December 2014)

概要

色々なマーカーの出力サンプルを作って見ました. データはこちら(http://matplotlib.org/api/markers_api.html)のページより取得しました. まず,プロット記号と説明を対にしたタプルを作成します. これを enumerate にいれて for でまわしプロットしていきます.

コードは以下のとおり.

from matplotlib import rcParams
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Arial']
import matplotlib.pyplot as pl

ml=(
('.','point'),
(',','pixel'),
('o','circle'),
('v','triangle_down'),
('^','triangle_up'),
('<','triangle_left'),
('>','triangle_right'),
('1','tri_down'),
('2','tri_up'),
('3','tri_left'),
('4','tri_right'),
('8','octagon'),
('s','square'),
('p','pentagon'),
('*','star'),
('h','hexagon1'),
('H','hexagon2'),
('+','plus'),
('x','x'),
('D','diamond'),
('d','thin_diamond'),
('|','vline'),
('_','hline'),
('$\sum$','TeX')
)
fig=pl.figure()
pl.xlim(0,6)
pl.ylim(25,0)
for i,k in enumerate(ml):
    pl.plot(1,i+1,linestyle='None',marker=k[0],markersize=10,markerfacecolor="w",markeredgecolor="b",markeredgewidth=1,label=k[1])
pl.legend(shadow=True,title='Marker',loc='upper right',ncol=2)
#pl.legend(shadow=True,title='Marker',loc='upper right', prop={'size',14},ncol=2)
pl.savefig('fig_marker.png',dpi=200)
pl.show()

下に示す行は,2行目をコメントにしています. 1行目は動くのですが,2行目は prop={'size',14} でエラーが出て動きません.なぜだろう?

pl.legend(shadow=True,title='Marker',loc='upper right',ncol=2)
#pl.legend(shadow=True,title='Marker',loc='upper right', prop={'size',14},ncol=2)

描画グラフ

pic
matplotlib作図:マーカーサンプル

文字コード変換 (13 December 2014)

概要

指定したディレクトリ以下にある全ての html ファイル,css ファイル,txt ファイルの文字コードを utf-8 に変更し同じファイル名で上書きします.

主な処理の流れは以下のとおり.

  • os.walk を用いて処理するファイルのフルパスのリストを作成します.os.walk は,ディレクトリの変更は行わないので,指定ディレクトリは python プログラムがある場所を除いてドライブを含めて指定します.
  • ファイルの入力は fin=codecs.open(f, 'r', 'shift_jis') というように,文字コードを指定して行っています.
  • 下の事例では html の head の charset および python ソース ( txt ファイル) の文字コード指定も utf-8 に書き換えています.
  • 基本は,shift_jis から utf-8 への変換用ですが,ソースが既に utf-8 になっている場合や euc-jp の場合に対応するため try: ..... excep: で試行錯誤しながらソースの文字コードを設定しています.
  • ファイルの出力は,fout=codecs.open(f, 'w', 'utf-8') というように,明示的に utf-8 で行います.

ソース

# coding: utf-8
import os
import codecs

str1s='<meta charset="shift_jis">'
str1u='<meta charset="utf-8">'
str2s='<meta http-equiv="Content-Type" content="text/html;charset=shift_jis">'
str2u='<meta http-equiv="Content-Type" content="text/html;charset=utf-8">'
str3s='#! coding: Shift_JIS'
str3u='# coding: utf-8'
str4s='# coding: Shift_JIS'
str4u='# coding: utf-8'

f_list=[]
#dir0='c:\DATA_WANtaro\Fixed_WANtaroHP'
dir0='d:\WAN1'
for root, dirs, fname in os.walk(dir0):
    for fname1 in fname:
        if '.html' in fname1: f_list=f_list+[os.path.join(root,fname1)]
        if '.css' in fname1: f_list=f_list+[os.path.join(root,fname1)]
        if '.txt' in fname1: f_list=f_list+[os.path.join(root,fname1)]

for f in f_list:
    print(f)
    try:
        fin  = codecs.open(f, 'r', 'shift_jis')
        data=fin.read()
        fin.close()
    except: pass
    try:
        fin  = codecs.open(f, 'r', 'utf-8')
        data=fin.read()
        fin.close()
    except: pass
    try:
        fin  = codecs.open(f, 'r', 'euc-jp')
        data=fin.read()
        fin.close()
    except: pass

    data=data.replace(str1s,str1u)
    data=data.replace(str2s,str2u)
    data=data.replace(str3s,str3u)
    data=data.replace(str4s,str4u)

    fout = codecs.open(f, 'w', 'utf-8')
    fout.write(data)
    fout.close()

pass

カレンダー (16,20 December 2014)

概要

ネットで Python で作成したカレンダーをアップしている方がいるのを見て,作って見たくなりました. 会社のPCにはガジェットがないので,役に立ちそうです.実際に,今日は何曜日だっけ?とか,今年は平成何年だったっけ?ということもよくありますしね. 12月16日に初アップしましたが,機能を改良して12月20日に再アップです. 機能を以下に示します.

  • プログラム起動日を含む月の月間カレンダー,もしくはプログラム起動日を含む年の年間カレンダーを表示するものです.
  • 付属としてデジタル時計もつけてみました.
  • 左上と右上の矢印を押すと,月もしくは年が変化します.
  • 真ん中の起動日の日付を押すと,起動日を含む月もしくは年のカレンダーに戻ります.

カレンダー作成の数学的準備として,ツェラーの公式による曜日の確定と,閏年の判定が必要になります.これらは wikipedia で確認できます. 特にツェラーの公式適用に当たっては,整数化処理をきちんと行わないと正確と曜日を正確に算定できないので注意. ちなみに,私の場合,月曜日から始まるカレンダーにしています.

月や年をボタンで変化させる場合,うまくいったと思ったら起動して16回目の変更で必ず挙動がおかしくなりました. このため,1回づつdestroyでカレンダーを表示しているFormを消去・再描画を繰り返すようにしています. 切り替わりの時一瞬画面が消えるような挙動になりますが,個人で使う上では問題視していません.

ソース

FilenameDescription
py_misc_calm.txt月間カレンダー表示スクリプト
py_misc_caly.txt年間カレンダー表示スクリプト

描画グラフ

pic
画面キャプチャー:カレンダー

デジタル時計と新しいwindowの表示 (16 December 2014)

概要

カレンダープログラムを改良するために必要な,基本的なコードをネット検索しました.

ソース

デジタル時計

こちらのページ (http://www.geocities.jp/m_hiroi/light/pytk07.html) 掲載のコードを参考にしました. strftime のformatはこちら (http://docs.python.jp/2/library/time.html) のページで確認できます.

# coding: utf-8
from tkinter import *
from time import *

def show_time():
    buff.set(strftime('%H:%M:%S (%a %d %b %Y)'))
    root.after(1000, show_time)

root = Tk()
buff = StringVar()
buff.set('')
Label(root,textvariable=buff,font='consolas 14 bold',relief=FLAT).pack()
show_time()
root.mainloop()

ボタンを押すと新しいwindowが表示される

# cording: utf-8
#ボタンを押すとwindowとメッセージが表示される
from tkinter import *

def show_win():
    sub_win = Toplevel()
    Label(sub_win, text = u'サンプルプログラムです').pack()

root = Tk()
root.option_add('*font', ('FixedSys', 14))
button = Button(root, text = "Button %d" % 1, command = show_win)
button.pack()

root.mainloop()

コマンドプロンプト画面を最小化してバッチファイルを実行 (17 December 2014)

概要

ファイラーなどから Python を実行する場合,実行中コマンドプロンプトの画面(黒画面)が立ち上がっています. カレンダープログラムの実行中など,黒画面が画面の大きな面積を占めるのはかっこよくありません. また手動で最小化するのもなんとなくかっこよくない,ということで,黒画面を最小化してコマンドを実行する方法を探しました.

2つのバッチファイルを用いる

下に掲載のバッチファイルは,こちらのページ (http://piyopiyocs.blog115.fc2.com/blog-entry-731.html) に掲載されているものです. start /MIN %~dp0a_digit.bat %* の a_digit.bat が実行するバッチファイルです. a_digit.bat の中で Python のファイル py_digit.py を実行しています.Python ファイル実行後,exit でコマンドプロプトを終了します.

@echo off
start /MIN %~dp0a_digit.bat %*
python py_digit.py
exit

1つのバッチファイルを用いる

試して見たところ,以下の1バッチファイルでコマンドプロンプト最小化実行が可能なようです. 結局は,dosコマンド「start /min」にPythonの実行コマンドを適用しただけです.

start /min python py_misc_calm.py
exit

選択したディレクトリ内のjpg・png画像表示 (26 December 2014)

概要

前から気になっていたテクニックで,ようやくやりたいことができた気がします. やっていることは以下のとおり.

  • プログラム起動でメニューを含むダイアログを表示
  • メニューを開き,ディレクトリを選択
  • 選択されたディレクトリ内の jpg および png ファイルを検索
  • 発見された全ての jpg および png ファイルのサムネイル(最大サイズ 150 x 150)を作成し,ボタンに貼り付けて window に表示
  • ボタンを押すと拡大された画像を表示.画像サイズは最大 600 x 600 になるよう調整

このスクリプト作成のため,以下のサイトが大変参考になりました.

フレームにスクロールバーを付けるhttp://effbot.org/zone/tkinter-autoscrollbar.htm
イメージのラベルへの表示 http://effbot.org/tkinterbook/label.htm
ボタン位置を変数で扱う https://www.daniweb.com/software-development/python/threads/424169/tkinter-center-a-grid-of-buttons

フレームにスクロールバーを付ける

簡単にできそうで,できない.結局上記メージに紹介されているコードをそのまま用いました.

イメージのラベルへの表示

これも単純にはいかない.このスクリプトではボタンニイメージを貼り付けていますが,やり方はラベルの場合と同じ.

ボタン位置を変数で扱う

上述 Webpage の方法

これも簡単にできそうで,できない.各ボタンに通し番号を付け,ボタンのオブジェクトと番号を辞書として保存し, これを用いてボタンを押したときに目的の画像ファイルを開くのに用いる.

lambda を使う方法

色々調べて結局,ommand=lambda kk=k: SHOW_PIC(kk) で関数SHOW_PICにボタンの番号を示す数値kkを渡すことに成功. これを使っています.この方法でも,ボタンをリストに保存しておくことが必要.

ソース

FilenameDescription
py_img_view.txt画像表示スクリプト

描画画像

pic
画面キャプチャー:実行時スクリーンショット

地下発電所概略形状描画 (28 December 2014)

概要

こういうのは,作ったらすぐアップしておかないとスクリプトがどこにいったか分からなくなってしまうのでアップ. Cross Section(1) は,矢印の描画法の勉強などを含め,これ1枚描くのに丸1日かかった. その他の3枚は,慣れてきたせいか,半日で完了. 両矢印のうまい描画方法がわからず,半分から右側と左側,というように2回の描画を行っている.もっと合理的にできないものか。。。

小技

3点を通る円の方程式

平面図を描くとき,妻壁を描くには3点を通る円の方程式を知らなければならない. 円の方程式は

x2 + y2 + c[0]*x + c[1]*y + c[2] = 0
c[0]*x + c[1]*y + c[2]*1 = - x2 - y2

とおけるから, 既知の3点を (x1,y1),(x2,y2),(x3,y3)として,連立方程式を作成・解くことにより,係数 C[0]、C[1]、C[2] を求めています. これらの係数が求まれば,2次方程式の解の公式により,(ここでは)y を与えることにより同一円上の x を求めることができます. このようなことを,簡単にできるところが Python のいいところですかね.

def CALX(c,y,pm):
    return (-c[0]+pm*sqrt(c[0]**2-4*(y**2+c[1]*y+c[2])))/2

def CALC(x1,y1,x2,y2,x3,y3):
    a=np.array([[x1,y1,1],[x2,y2,1],[x3,y3,1]])
    b=np.array([-x1**2-y1**2,-x2**2-y2**2,-x3**2-y3**2])
    c=np.linalg.solve(a,b)
    return c

array を連結して行ベクトルを作る

下の事例では,arange を用いて y1,y2,y3,y4 を作り,これを連結し,最後に [339.400] という数値を加えて yy を作っています.

y1=np.arange(339.400,320.900,-dy)
y2=np.arange(320.900,314.8-dy,-dy)
y3=np.arange(314.8,320.900,dy)
y4=np.arange(320.900,339.400,dy)
yy=np.r_[y1,y2,y3,y4,[339.400]]

ソース

FilenameDescription
py_ps_sec1.txt断面図表示スクリプト(1)
py_ps_sec2.txt断面図表示スクリプト(2)
py_ps_sec3.txt縦断面図表示スクリプト
py_ps_sec4.txt平面図表示スクリプト

描画画像

pic pic
Cross Section (1) Cross Section (2)
pic
Longitudinal Section
pic
Plan

画像の一部を切り取り保存 (01 January 2015)

概要

facebook などに投稿する際,画像の一部だけを切り取りたい場合がある. この場合,今までは「ペイント」か ImageMagick の crop で行っていたが,せっかくなのでPythonで行って見ることにした.

やっていることは以下のとおりです.

  • ダイアログから対象とする画像ファイルを選択
  • 対象画像ファイルをmatplotlibに読み込み表示
  • 表示された画像に対し,切り取りたい部分の左上と右下をマウスでクリック
  • クリックされた部分を囲んだ画像を表示
  • 切り取られた画像を保存すると同時に表示
  • 保存する画像ファイル名は,「fig_crop_」を元画像ファイル名の先頭に加えたもの

matplotlib上でマウスの座標を取得するテクニックは,こちら (http://d.hatena.ne.jp/Megumi221/20111212/132369201)を参考にしました. ただし,Python3では掲載されているコードはすぐ終了してしまい,うまくいきません. 私は以下のように修正したものをベースに用いています.

import matplotlib.pyplot as plt
import numpy as np

def onclick(event):
    print('button=%d, x=%d, y=%d, xdata=%f, ydata=%f' %(
    event.button, event.x, event.y, event.xdata, event.ydata))

fig = plt.figure()
plt.plot(np.random.rand(10))
fig.canvas.mpl_connect('button_press_event', onclick)
plt.show()

使い方

  • (1) 起動するとファイルダイアログと空のリストブックスが現れます.
  • (2) ダイアログタから対象画像ファイルを選択するとダイアログは消え,リストボックスにファイル名・パスなどが表示されます.
  • (3) リストボックスを終了させると選択した画像ファイルが表示されるので,切り取りたい部分の左上と右下をクリックし画像を閉じます.
  • (4) 選択した範囲が赤い枠で囲まれた画像が現れますのでそのまま終了します
  • (5) 選択範囲の画像が表示されるので確認後終了します.この画像が表示された時点で既に切り取られた画像の保存は完了しています

ソース

FilenameDescription
py_img_crop.txt画像切り取り保存スクリプト

描画画像

pic pic
スクリーンショット (1)スクリーンショット (2)
pic pic
スクリーンショット (3)スクリーンショット (4)
pic
スクリーンショット (5)

クリップボードの画像データを保存する (01 January 2015)

今までは,PrtScで取得したスクリーン画像をペイントで加工していましたが,何かと不便なので,Pythonでできないか調べて見ました.

色々検索して見ましたが,結局 PIL を用いて,以下のスクリプトでできます.これはこちら (http://stackoverflow.com/questions/7045264/how-do-i-read-a-jpg-or-png-from-the-windows-clipboard-in-python-and-vice-versa) のサイトからいただきました. 一応クリップボードに画像がなかった場合のエラー停止を回避しています.

from PIL import ImageGrab
import tkinter.messagebox
im = ImageGrab.grabclipboard()
try:
    im.save('somefile.png','PNG')
    tkinter.messagebox.showinfo('showinfo', 'Image was saved as "somefile.png" successfully')
except:
    tkinter.messagebox.showerror('showerror', 'No image in clipboard')

前出「画像の一部を切り取り保存 (01 January 2015)」のスクリーンショットはこれを使って画像化しています. ただし画像サイズが少し大きかったので,以下のImageMagickコマンドでリサイズしています.

convert -resize 600x -quality 100 -verbose fig_crop1.png fig_cr1.png
convert -resize 600x -quality 100 -verbose fig_crop2.png fig_cr2.png
convert -resize 600x -quality 100 -verbose fig_crop3.png fig_cr3.png
convert -resize 600x -quality 100 -verbose fig_crop4.png fig_cr4.png
convert -resize 600x -quality 100 -verbose fig_crop5.png fig_cr5.png

画像の一部を切り出す時は,画像を見ながらクリックして範囲を定める方法が便利ですが, 機械的にフォーマットを変えるとか,リサイズするとかだと ImageMagick を用いるのが楽ですね.

棒グラフ作成 (05 January 2015)

概要

普段はめったに作らないのだが,棒グラフ作成の必要性が生じたため,雛形を作成した.

凡例が複数列になる場合,デフォルトだと縦方向が優先されて配置されますが,これを横方向優先で配置する方法を,こちら (http://stackoverflow.com/questions/10101141/matplotlib-legend-add-items-across-columns-instead-of-down) のサイトからいただきました.誰かが同じようなことをやろうとしてつまづき,誰かがそれに答えているのですね.

ソース

FilenameDescription
py_fig_bar2.txt単純な横棒グラフ作成スクリプト
py_fig_bar3.txt少し複雑な横棒グラフ作成スクリプト

描画画像

pic pic
単純な横棒グラフ少し複雑な横棒グラフ

円グラフ作成 (01 January 2015)

概要

棒グラフのついでに,円グラフのスクリプトもgetしておきました.こちら (http://matplotlib.org/examples/pie_and_polar_charts/pie_demo_features.html) のサイトを参照.

ソース

import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Arial'
plt.rcParams['font.size'] = 14

# The slices will be ordered and plotted counter-clockwise.
labels =['Frogs', 'Hogs', 'Dogs', 'Logs']
sizes = [15, 30, 45, 10]
colors = ['yellowgreen', 'gold', 'lightskyblue', 'lightcoral']
explode = (0, 0.1, 0, 0) # only "explode" the 2nd slice (i.e. 'Hogs')

plt.pie(sizes, explode=explode, labels=labels, colors=colors,
        autopct='%1.1f%%', shadow=True, startangle=90)
# Set aspect ratio to be equal so that pie is drawn as a circle.
#plt.axis('equal')
plt.gca().set_aspect('equal')
plt.title('Pie Chart', size=24)

fnameF = "fig_pie.png"
plt.savefig(fnameF)
plt.show()
  • startangle=0は水平.startangleより反時計回りに描画が進むことに注意.
  • explodeは切り取る要素とその距離.ここでは2番目の要素「Hogs」が切り取られている
  • グラフを円形とするため,aspectを指定する.

描画画像

pic

HSV指定の色見本 (12 January 2015)

概要

HSV指定の色見本をmatplotlibで作成しました.実際の描画には HSV => RGB => hex の変換を行い,16進数で色指定を行っています.

ソース

FilenameDescription
py_misc_col_hsv.txtHSV指定の色見本作成スクリプト

描画画像

pic

tkinter の Label による色見本 (12 January 2015)

概要

tkinter の Label による色見本を作って見ました.このように表のようなものを作るには Label を用いるのが楽です. ただし Label は1行しか表示できないらしい. matplitlibなら座標を指定して描画するのでどこにでも線やテキストを入れられるのだけれど.

やっていることは以下のとおり.

  • 色はHSVにて指定しますが,これをRGB => 16進数に変換,描画時は.16進数による色指定を行っています.リスト(1)のスクリプト
  • matplotlibの使用可能色名はリスト(5)のとおり.HSVで指定した色に対し,matplotlibの有する色のうち最も近いものを表示するスクリプトがリスト(2)とリスト(3)です
  • リスト(2)では,HSVで指定した色の16進表示値と,matplotlib使用可能な色名の16進値の差の絶対値が最小となるものを最も近い色と判定し,これを表示しています
  • リスト(3)では,HSVで指定した色のHSV値と,matplotlib使用可能な色名のHSV値の距離の差が最小となるものを最も近い色と判定し,これを表示しています
  • リスト(4)はリスト(3)と同様に選定した色について,色名を表示しています
  • リスト(6)は,HSVで指定した色のHSV値と,リスト(7)のHTML color name のHSV値の距離の差が最小となる色を,色名で表示しています

ソース

FilenameDescription
py_misc_col_label.txtリスト(1):HSV指定による色見本作成スクリプト
py_misc_col_dic0.txtリスト(2):HSV指定色をmatplot色名から探す(1)
py_misc_col_dic1.txtリスト(3):HSV指定色をmatplot色名から探す(2)
py_misc_col_dic2.txtリスト(4):HSV指定色をmatplot色名から探す(3):色名表示
colorname.txtリスト(5):matplotlibの保有する色名リスト
py_misc_col_dic3.txtリスト(6):HSV指定色をHTML color nameから探す(3):色名表示
col_name_html.txtリスト(7):HTML color name リスト

描画画像

pic
HSV指定による色見本:リスト(1)による
pic pic
16進数の差:リスト(2)によるHSV空間距離:リスト(3)による
pic
HSV空間距離より選定したもっとも近い色(色名表示):リスト(4)による (matplotlib color name)
pic
HSV空間距離より選定したもっとも近い色(色名表示):リスト(6)による (HTML color name)

処理時間計測 (18 January 2015)

概要

計算処理時間計測のスクリプトです.

ソース

from datetime import datetime

t_start = datetime.now()
for i in range(10000):
    for j in range(10000):
        pass
t_stop = datetime.now()

print(t_stop - t_start)
print(t_start)
print(t_stop)

出力は下のような感じ.forループ 10000 x 10000 で 5.38秒でした.

0:00:05.382009
2015-01-18 07:27:56.900580
2015-01-18 07:28:02.282589

Pythonプログラムを実行するボタン (18 January 2015)

概要

いくつかのPythonプログラムを実行するボタンを配置したものを作って見ました. 同時に2個以上のプログラムの実行はできません.

pic

ソース (py_button.py)

from tkinter import *
import subprocess

def cmd1():
    subprocess.call('python py_misc_calm.py') 
def cmd2():
    subprocess.call('python py_misc_caly.py') 
def cmd3():
    subprocess.call('python py_img_view.py') 
def cmd4():
    subprocess.call('python py_img_PrtSc.py') 
def cmd5():
    subprocess.call('python py_img_crop.py') 
def cmd6():
    subprocess.call('python py_misc_col_label.py') 

root = Tk()
f0 = Frame(root)
Button(f0,text=u'カレンダー(月)',command=cmd1).pack(fill=BOTH)
Button(f0,text=u'カレンダー(年)',command=cmd2).pack(fill=BOTH)
Button(f0,text=u'画像表示',command=cmd3).pack(fill=BOTH)
Button(f0,text=u'スクリーンショット',command=cmd4).pack(fill=BOTH)
Button(f0,text=u'画像クロップ',command=cmd5).pack(fill=BOTH)
Button(f0,text=u'HSV色見本',command=cmd6).pack(fill=BOTH)
f0.pack(fill = BOTH)

root.mainloop()

実行用バッチファイル

start /min python py_button.py
exit

HTML: &<>” 変換 (19 January 2015)

概要

webページにプログラムやhtmlタグを含む文書をアップする時,<pre>を使っても「&」や「<」を「&amp;」や「&lt;」に変換しないと予期せぬ表示となる場合があります. このため,これらを自動変換するスクリプトを作りました.今まではエディタで置換していましたが一挙に変換するものです. 使い方はだいたい以下のとおり.

  • プログラム起動.プログラム起動前にクリップボードに取り込んでいても良い
  • 変換したいスクリプトなどをクリップボードに取り込む
  • 「convert」ボタンを押すとテキストボックスが別画面で立ち上がり変換後の文字列が表示される
  • この時クリップボードには既に変換された野時列が格納されているため,そのまま必要箇所にペーストする
  • クリップボードにデータがないときはその旨がメッセージボックスで表示される
  • 変換されるのは,半角「&」「<」「>」「”」の4文字のみです

ソース (py_misc_gtconv.py)

from tkinter import *
import tkinter.scrolledtext
import tkinter.messagebox

def cmd1():
    sub_win = Toplevel()
    editArea = tkinter.scrolledtext.ScrolledText(sub_win,width=40,height=20)
    editArea.pack(padx=10, pady=10, fill=BOTH, expand=True)
    try:
        cb0 = root.clipboard_get()
        cb1 = cb0.replace('&', '&amp;')
        cb2 = cb1.replace('<', '&lt;')
        cb3 = cb2.replace('>', '&gt;')
        cb4 = cb3.replace('"', '&quot;')
        editArea.insert(INSERT,cb4)
        root.clipboard_clear()
        root.clipboard_append(cb4)
    except:
        tkinter.messagebox.showerror('showerror', 'No data')

def cmd2():
    root.clipboard_clear()

root=Tk()
f=Frame(root)
f.pack(fill='both', expand='yes')
Button(f,text='Convert',command=cmd1).pack(side=LEFT)
Button(f,text='Clipboard clear',command=cmd2).pack(side=LEFT)

root.mainloop()
inserted by FC2 system