Python | MasaTANO HP

moduleのインストール方法
Anacondaのインストール
テキストデータのIO
バイナリデータのIO
Mask/NaNの処理
evalの利用
dictによるファイル名の管理
awkの利用
PyGrADS
GmtPy
Fortran modules の利用（f2py）
Fortran modules の利用（ctypes）
pandas
matplotlib設定のやり方
plot(matplotlib)のまとめ
Basemap

moduleのインストール方法

モジュールのインストール方法(Ubuntu)

$ sudo apt-get install python-numpy
$ sudo apt-get install python-matplotlib
$ sudo apt-get install python-mpltoolkits.basemap

モジュールのインストール方法(Mac)

$ pip3 install pylab

モジュールのインストール方法(Anaconda)

$ conda install pylab

念のため、インストールできているかの確認。ヴァージョンが確認できればOK!!

import numpy
print (numpy.__version__)

使った感じ、condaでインストールしていくのが良いだろう（企業が提供しているという点が気になるが）。 pythonは様々なバージョンがあり、 pipだと環境を破壊してしまう場合があるので、 condaで統一するようにする。

anacondaのインストール

Pyenvのインストール

$ brew install pyenv

export PYENV_ROOT=~/.pyenv
export PATH=$PYENV_ROOT:$PATH
eval "$(pyenv init -)"

Anacondaのインストール

$ pyenv install -l | grep anaconda
$ pyenv install anaconda3-5.0.1
$ pyenv local anaconda3-5.0.1

. $PYENV_ROOT/versions/anaconda3-5.0.1/etc/profile.d/conda.sh

$ conda --version
conda 3.19.1

仮想環境の構築

$ conda create -n py2 python=2.7.14 anaconda=5.0.1
$ conda activate py2
$ conda create -n py3 python=3.6 anaconda=5.0.1
$ conda activate py3

$ connda install ipython

$ conda activate py2
(py2)$ conda install numpy
(py2)$ conda install scipy
(py2)$ conda install pandas
(py2)$ conda install matplotlib
(py2)$ conda install -c anaconda basemap
(py2)$ conda install basemap-data-hires
conda install -c conda-forge pygrib

$ conda activate py3
(py3)$ conda install numpy
(py3)$ conda install scipy
(py3)$ conda install pandas
(py3)$ conda install matplotlib
(py3)$ conda install -c anaconda basemap
(py3)$ conda install basemap-data-hires

moduleをインストールする際に、condaを使うか、pipを使うかをはっきりさせる。
pipだとcondaの環境を破壊してしまうため、condaを使う方が良い。ちなみにcondaの環境が壊れた場合は大人しくanacondaを再インストールする。

ファイルやディレクトリ名の取得

import os
# ファイル名の取得
os.path.basename(fDat)

# フォルダ名の取得
os.path.dirname(fDat)

# 拡張子の取得
root,suf = os.path.splitext(fDat)

suffixにはドットも含まれる。

テキストデータのIO

スペース区切りのテキストデータ（全て数字データ）の読み方

from numpy import *
r1dat = loadtxt (FDAT)

バイナリデータのIO

Read

1次元バイナリデータ(単精度float32)の読み込み

import numpy as np
Dat = np.fromfile("test.bin", np.float32)

1次元バイナリデータ(倍精度floati64)の読み込み

import numpy as np
Dat = np.fromfile("test.bin", np.float64)

2次元バイナリデータの読み込み。

import numpy as np
NX, NY = 360, 180
Dat = np.fromfile("test.bin", np.float32).reshape([NY, NX])

Dat = np.fromfile("test.bin", np.float32).reshape([NZ, NY, NX])

Dat = np.fromfile("test.bin", np.float32).reshape([-1, NX])

エンディアン

Write

バイナリデータの書き込み

import numpy as np
Out.astype(np.float32).tofile("out.bin")

Mask/NaNの処理

Mask処理

Mask処理。１次元配列Datが負の値を取るとき、欠損値(-9999)を代入するという計算。
値が大きいものをmaskしたい場合はlessではなくgreaterを用いる。

import numpy as np
# やり方1
Dat = np.ma.masked_less(Dat,0.0).filled(-9999.0)
# やり方2
Dat[np.where(Dat<0)] = -9999.0

なお、やり方１と２ではobject IDの取扱いが異なる。そのため、複数回maskすると、１の方法ではmaskがうまくできなかったりする。基本的に2の方法で統一するようにする。

Nan処理

math,numpy,pandasにnanを処理できるmethodが用意されている。ここではnumpyについての取扱を。

import numpy as np
dat = np.array([1., 2., nan, 4., 5.])

# nanの置換え
dat[np.isnan(dat)] = 0.0

# nanを除いたnpmeanの計算
np.nanmean(dat)

nanの代入はnp.nanでOK。ただし、np.nanはfloatの扱いなので、nparrayがintegerだと代入できない。

evalの利用

GCM01 = ["MIROC5", 1]
NUM = "01"
a = eval ("GCM" + NUM + "[1]")

dictによるファイル名の管理

fDatName = "%(GCM)s/%(RCP)s/flddph_%(YYYY)s.dat"
Dict = {"GCM":"MIROC5", "RCP":"rcp85", "SSP":"SSP3", "YYYY":"2100"}
fDat = fDatName % Dict

awkの利用

pythonスクリプト内でのawkの利用。~~awkが使えるのはやっぱり便利。基本はcommandのgetoutputのようにする。~~と、思ったがpandasのほうが遥かに便利なので、そちらを利用する。

import commands as cmd
LIST = cmd.getoutput("awk '$1==\"SSP1\"&&$3==\"2010\"{print $2,$4}' " + FDAT)

文字列を検索する際には、ダブルクォーテーションの前にバックスラッシュをつける。

変数で条件を絞りたい場合は下記の通り。

import commands as cmd
LIST = cmd.getoutput("awk '$1==\""+SSP+"\"&&$3=="+YEAR+"{print $2,$4}' " + FDAT)

numpyに変換したい場合は下記のように。

import commands as cmd
awkcmd = "awk '"+condition+"{print "+out+"}' " + c0fdat
tmpres = cmd.getoutput(awkcmd)
tmp1 = np.array(tmpres.split()).reshape([-1,2]).astype(np.float32)

追記：python 3以降、commandsがなくなったのでsubprocessでやるしかなくなる。

oid = 4
COMMAND = ["awk", '$1=="historical" && $2=='+str(oid)+' {print $3, $4}', FVASS]
result = subprocess.run(COMMAND)

PyGrADS

Install

macでもlinuxでも共通。

$ cd ${your_work_dir}
$ wget https://sourceforge.net/projects/opengrads/files/python-grads/1.1.9/pygrads-1.1.9.tar.gz
$ tar xvf pygrads-1.1.9.tar.gz
$ cd pygrads-1.1.9
$ python setup.py install

Install先を指定したい場合は、install時にprefixを指定する。

Install and start pygrads

pygradsの起動。GaNumでもGaLabでも可能。起動時にwindowをoffにさせることできる。


from grads import GaNum, GaLab
ga = GaNum('grads -l -b')
ga = GaLab(Bin='grads',Window=False,Echo=True,Verb=0,Port=False)

Open CTL file

どちらでも開ける。sdfopenも可能。


ga.open(FCTL)
ga('open %s' % FCTL)

Set dimentions and resolution

ga.query('dims')で緯度経度、グリッド情報などを収集できる。


dim = ga.query('dims')

# X座標のグリッド数
ga('set x 1 %s' % dim.nx)

# 緯度・経度の取得
InLonMin = dim.lon[0]
InLonMax = dim.lon[1]
InLatMin = dim.lat[0]
InLatMax = dim.lat[1]

# 緯度・経度の間隔の取得
InLonInt = (InLonMax-InLonMin)/float(NX)
InLatInt = (InLatMax-InLatMin)/float(NY)

Get variables


# undef valueの設定
ga('set undef -9999.0')

# 緯度・経度グリッドデータの取得
Lon = ga.exp('lon')
Lat = ga.exp('lat')

# 変数グリッド・データの取得
Out = ga.exp(VAR)

Draw

windowを表示させている場合、"d var"で変数を出力可能。GaLabを使った例として下記を示す。

import matplotlib.pyplot as plt

ga.blue_marble("on")
ga.basemap('ortho',opts=(130.,40.))
ga.contour("temp")
del ga
plt.title("Surface Temperature (Jan1958)")
plt.savefig("galab_slp.png")

参考：oceansciencehack.blogspot.jp/2011/02/galab-in-python-interface-to-grads.html

小技

GrADS形式の日付データの取得

import datetime
gadate = datetime.datetime(Year, Month, Day).strftime('%d%b%Y')

Caution

macにpygradsをインストールした場合、下記の問題が発生した。

python上でpygradsをimportできない。X-windowが立ち上がらない。
python上でpygradsをimportまではできるが、そこからさきのコマンドが全てpipe errorとなる
X-windowが一瞬だけ起動し、すぐに消える。

どの問題も、pythonからgradsを起動する際(ga=GaNum("grads -l"))のバイナリに問題があるようだ。 opengradsの中にあるDarwinのgradsを使用すれば、無事にX-windowが立ち上がった。
# ただ、困ったことに問題のgradsバイナリは、macで普通に使えていた。
macでpygradsを使用する場合はgradsのPATHに注意する必要がある。 brewでgradsをインストールしてもOK。

GmtPy

Install

インストールは本家HPを参照：http://emolch.github.io/gmtpy/http://emolch.github.io/gmtpy/

GrADS to GMT

gradsから直接GRDファイルを作成するスクリプト。


# Install and start modules
from grads import GaNum
from gmtpy import GMT

ga = GaNum('grads -l -b')
gmt = GMT()

# Open CTL file
ga.open(FCTL)

# Set dimentions and resolution
dim = ga.query('dims')
ga('set x 1 %s' % dim.nx)
dim = ga.query('dims')
InLonMin = dim.lon[0]
InLonMax = dim.lon[1]
InLatMin = dim.lat[0]
InLatMax = dim.lat[1]
InLonInt = (InLonMax-InLonMin)/float(NX)
InLatInt = (InLatMax-InLatMin)/float(NY)
Resol = (InLonInt, InLatInt)
Range = (OutLonMin, OutLonMax, OutLatMin, OutLatMax)

# Set Ranges
X1 = int((OutLonMin-InLonMin)/InLonInt)
X2 = int((OutLonMax-InLonMin)/InLonInt)
Y1 = int((OutLatMin-InLatMin)/InLatInt)
Y2 = int((OutLatMax-InLatMin)/InLatInt)
ga('set x %s %s' % (X1, X2))
ga('set y %s %s' % (Y1, Y2))

# Get variables
ga('set undef -9999.0')
Lon0 = ga.exp('lon')
Lat0 = ga.exp('lat')
Out0 = ga.exp(VAR)

# Modify for GMT input format
# ここで欠損値を省くことで、xyz2grdのIOおよび計算量が少なくなるのでやったほうが良い。
Lon0 = Lon0.reshape(-1)
Lat0 = Lat0.reshape(-1)
Out0 = Out0.reshape(-1)
index = np.where( Out0!=-9999. )[0]
Lon = Lon0[index]
Lat = Lat0[index]
Out = Out0[index]
del Lon0, Lat0, Out0

# Make GRD file
# out_discard=Trueはxyz2grdをする上で必須。
# out_discard=Falseは"-K -O"をxyz2grdコマンドに追加するのと同じ意味。
gmt.xyz2grd( in_columns=[Lon,Lat,Out], R=Range, I=Resol, G=FGRD, out_discard=True)

Fortran moduleの利用（f2py）

これまで開発してきたFortran modulesをpythonで利用する方法についてのまとめ。 pythonから呼ぶmoduleとしてはctypesとf2pyとがある。ctypesはFortran modules に配列を渡す際に不一致(配列が0からはじまるか、1からはじまる)が生じるため、 Fortran側のmodulesを修正する必要があるため、比較的修正が少ないf2pyの利用方法についてまとめる。

Test programe

Fortran側のテスト・モジュール(f2py_test)。2つの配列の和を返す。欠損値p0misがある場合は欠損値を返す。

module f2py_test
  implicit none
contains

  subroutine add(r1out, r1dat1, r1dat2, p0mis)
  ! Arguments
  real(4),         intent(inout)         :: r1out(:)
  real(4),         intent(in)            :: r1dat1(:)
  real(4),         intent(in)            :: r1dat2(:)
  real(4),         intent(in)            :: p0mis
  ! Local
  integer                               :: n0l,  i0l
  real(4)                               :: r0dat1,  r0dat2

  n0l = size(r1dat1)
  do i0l = 1,  n0l
    r0dat1 = r1dat1(i0l)
    r0dat2 = r1dat2(i0l)
    if ( r0dat1 /= p0mis .and. r0dat2 /= p0mis ) then
      r1out(i0l) = r0dat1 + r0dat2
    else
      r1out(i0l) = p0mis
    endif
  enddo

  end subroutine
end module

コンパイル


$ ifort -c -fPIC f2py_test.f90
$ f2py -m f2py_test -h f2py_test.pyf f2py_test.f90
$ f2py --fcompiler=intelem --f90flags='-O2 -fPIC' -c f2py_test.pyf f2py_test.o

pythonからの呼び出し方

先にFortran modulesがあるディレクトリにPYTHONPATHを通しておく。 modulesがあるディレクトリが~/testだとすると、~/.bash_profileに下記の行を追加して、sourceすればOK。

export PYTHONPATH=$PYTHONPATH:~/test

次にpython側での呼び出し方。Fortranは引数の配列の精度を気にするのでFortran側に合わせる。 (Fortranの方は、はじめから倍精度でまとめておくほうがよいかも。)

import f2py_basic as mt
mt.f2py_test.add(result, r1dat1, r1dat2, p0mis)

どのようなFortran modulesを扱うべきか

pythonはdat[il]=dat[il]*frc[il]のように, 配列を指定した計算がかなり苦手。自分が扱う分野だと単純なアップスケール・ダウンスケール、それとmapデータから国ごとに値を集計するなど。このような計算は時間がかかるのでFortranに任せた方がよい。

ちなみに、pythonとFortranとではIOの時間に大きな差は見られなかった。 pythonでレコードを指定したデータの読み書きもできるが、個人的にはFortranの方が楽に感じる。

Fortran moduleの利用（ctypes）

ctypesによる方法まとめ。ここがすごく参考になった（http://d.hatena.ne.jp/ignisan/20121017）

Test module

ctypesの場合、引数の配列数を明示的に指定する必要がある。 Fortran同士でメインプロジェクトからモジュールに配列を与える場合、モジュールへintent （in）で配列を渡すと配列サイズも自動的に与えてくれるが、 pythonからFortranのモジュールを呼び出した場合、配列サイズが０になってしまう。そのため、どうしても引数に配列サイズを取る必要がある。以下、例。

$ cat fmodule_add.f90
subroutine fmodule_add(n0l, d1a, d1b)
  implicit none
  integer, intent(in)  :: n0l
  real(8), intent(in)  :: d1a(0:n0l-1)
  real(8), intent(out) :: d1b(0:n0l-1)

  d1b = d1a + d1b

end submodule

コンパイル

成功するとfmodule_add.so ファイルが作成される。

$ ifort -shared -fPIC -o fmodule_add.so fmodule_add.f90

pythonからの呼び出し

pythonプログラムは下記のとおり。呼び出しが非常にめんどい。コンパイルは楽なんだけど。

from ctypes import *
import numpy as np

add_np = np.ctypeslib.load_library("fmodule_add.so",".")

add_np.fmodule_add_.argtypes = [ POINTER(c_int64), np.ctypeslib.ndpointer(dtype=np.float64), np.ctypeslib.ndpointer(dtype=np.float64)]

add_np.fmodule_add_.restype = c_void_p
a = np.arange(0.,10.,dtype=np.float64)
b = a*2
size = byref(c_int64(b.size))

add_np.fmodule_add_(size, a, b)

print b

pandas

Import modules

import pandas as pd
import matplotlib.pyplot as plt

データフレームの作成

データフレームを直接作成する。

df = pd.DataFrame({
 'A' : [1,2,3,4,5],
 'B' : [6,7,8,9,10]
})

numpyからデータフレームを作成する。

data = np.fromfile(fDat).reshape([-1, 2])
df = pd.DataFrame(data, columns=("ID","Out"))

データの読み込み

コンマ区切りのファイルの読み込み

df = pd.read_csv(fDat)

スペース区切りのファイルの読み込み

df = pd.read_csv(fDat,delim_whitespace=True)

スペース区切りのファイルの読み込み

df = pd.read_excel('analysis.xlsx', sheetname='anl1')

データフレームの書き出し

スペース区切りのファイルの書き出し

df.to_csv(fOut, sep=' ', index=False)

データの追加/削除

データフレームへのデータの追加

df['NewVar']=df['Var1']/df['Var2']

データフレームからのデータの削除

df.drop('Var1',axis=1)

重複部分の削除

df.drop_duplicates()

欠損値NaNの除去

欠損値が含まれている行ごと除去する。

df = df.dropna()

ラベルを指定して欠損値を除去する。

df = df.dropna(subset=['ID','Value'])

条件付きデータの抽出

ある列のデータを取得する。

df['Value']

ある複数列のデータを取得する。

df[['ID','Value']]

ある条件に一致するデータを取得する。

df[df['ID']==776]

ある複数条件に一致するデータを取得する。

df[(df['ID']==776) & (df['Year']==2004)]

ある複数条件に一致するデータのうち、特定の変数を取得する。

df[(df['ID']==776) & (df['Year']==2004)] ['Value']

IDのリストに含まれているIDと同じデータを取得する。

df[df['ID'].isin(list(npcode.reshape([-1])))]

ちなみに、and(&)、or（|）、not(~)だそうな。

for ループ

縦方向にループ

for index, item in df.iterrows:
    item['Var']

groupby

グループ化した要素のカウント

dftmp = df[['ID', 'Value']]
grouped = dftmp.groupby('ID')
dftmp2 = grouped.count()

グループ化した要素同士の相関係数の算定

grouped = df.groupby(['ID'], )
dfcorr = grouped.apply(lambda x: x['A'].corr(x['B'])).reset_index()

2つのデータフレームの結合

めちゃ簡単。Excelより楽。

共通のラベルを持つデータフレームの結合

pd.merge(df1, df2)

複数条件が一致する場合のデータフレームの結合

pd.merge(df1, df2, on=['key1', 'key2'])

参考HP：http://sinhrks.hatenablog.com/entry/2015/01/28/073327

pivot table

条件付き合計や特定IDの平均値などを簡単に算定できる。

dfglobal =  pd.pivot_table(df, index='Year', margins=True, aggfunc=np.sum)

ポイントはmarginsをTrueにして、aggfuncでどのように計算するのかを指定する。marginsをTrueにしないと、条件に合うデータが全て出力される。

DataFrameの有効桁数やら書式やらの表示変更方法

整数と少数での表記

pd.options.display.float_format = '{:.2f}'.format

有効桁数の変更

pd.options.display.float_format = '{:.3g}'.format

Draw figures

散布図の描画

df.plot.scatter(x='label_name1',y='label_name2')

２つのデータを重ねて描画

df.plot.scatter(x='label_name1',y=['label_name2', 'label_name3'])

複数のグラフを描画

nrows=4 ; ncols=3
fig, axes = plt.subplots(nrows=nrows, ncols=ncols, figsize=(12, 9))

for irow in range(nrows):
    for icol in range(ncols):
        axes[irow,icol].set_xlim([-9,9])
        axes[irow,icol].set_xticks([-9,-6,-3,0,3,6,9])
        axes[irow,icol].set_yscale('log')
        axes[irow,icol].set_ylim([1,1000])
        axes[irow,icol].set_title('Test')
        dfmean.plot.scatter(ax=axes[irow, icol], x='X',y='Y')

Convert to numpy format

nparray = df.as_matrix()

matplotlib設定のやり方

ややこしい。

軸の設定

範囲指定

plt.xlim([0,10])

範囲指定(日付)

import datetime
plt.xlim([datetime.datetime(2010,1,1),datetime.datetime(2010,12,31)])

対数軸

plt.xscale('log')

軸メモリの書き換え

plt.xticks([0,1,2], ["a","b","c"])

ラベル名の設定

plt.plot(X, Y, label="LabelName")

凡例の設定

plt.legend(loc='upper top')

locに指定可能な値は次のとおり。 (best, right, center, left) (lower, center, upper) * (right, center, left)

文字サイズの変更

plt.title("Title", fontsize=18)
plt.xlabel("xlabel", fontsize=18)
plt.ylabel("ylabel", fontsize=18)
plt.legend(fontsize=18)
plt.tick_params(labelsize=18)

まとめて変更する場合。

plt.rcParams["font.size"] = 18

出力する画像ファイルの位置調整

ここはGMTより楽でよい。

plt.savefig(fPng, dpi=300, bbox_inches="tight")

pltのリセット

何回も図を描くと重ね書きされてしまう。その場合はplt.figure()でリセット。

軸に平行な特定の線を引く。

plt.axhline(y=2, lw=2, color='red')
plt.axvline(x=2, lw=2, color='blue')

plot(matplotlib)のまとめ

散布図

基本

plt.scatter(Dat,c='red')

plt.scatter(Dat,c=value,cmap=cmap)

半透明化

plt.scatter(Dat,alpha=0.5)

エラーバー

plt.scatter(Dat,xerr=0.3,yerr=yerr)

ヒストグラム

基本

plt.hist(Dat,color='red',bins=10)

正規化

plt.hist(Dat,normed=True)

半透明化

plt.hist(Dat,alpha=0.5)

軸に平行な任意の線を引く

plt.axhline(y=10, lw=1, color='red')
plt.axvline(x=30, lw=2, color='blue')

Original cmapの設定方法

GMTよりめんどい。とりあえずは下記のHPを参考に、自動でcmapを作成するmodueを作成。ただし、今のところ不等間隔には対応できてないので要改良(2018-12-07)。 https://qiita.com/kenmatsu4/items/fe8a2f1c34c8d5676df8$

def generate_cmap(slef, colors='rainbow', inum=10):
from matplotlib.colors import LinearSegmentedColormap
values  =  range(len(colors))
vmax  =  np.ceil(np.max(values))

color_list  =  []
for v,  c in zip(values,  colors):
        color_list.append( ( v/vmax, c) )
return LinearSegmentedColormap.from_list('custom_cmap', color_list)._resample(inum)

Basemap

インストール

anacondaでインストール。こちらを参照。

Import

from mpl_toolkits.basemap import Basemap

基本設定

Basemapの設定。描画する範囲や図法など。

m=Basemap(projection='cyl',\
          resolution='h',\
          llcrnrlon=LonMin,\
          llcrnrlat=LatMin,\
          urcrnrlon=LonMax,\
          urcrnrlat=LatMax)

海岸線、国境線などの引き方

国境線

m.drawcountries(color='black',linewidth=0.5)

海岸線

m.drawcoastlines(color='k',linewidth=0.4)

河川

m.drawrivers(color='b',linewidth=0.3)

経度・緯度

LonMin, LonMax, LatMin, LatMax = -180, 180, -90, 90
m.drawparallels(range(LatMin,LatMax+1,1),\
                color='gray',\
                linewidth=0.5,\
                labels=[1,0,0,0],\
                fontsize=12)
m.drawmeridians(range(LonMin,LonMax+1,1),\
                color='gray',\
                linewidth=0.5,\
                labels=[0,0,0,1],\
                fontsize=12)

Contents

moduleのインストール方法

anacondaのインストール

ファイルやディレクトリ名の取得

テキストデータのIO

スペース区切りのテキストデータ（全て数字データ）の読み方

バイナリデータのIO

Read

Write

Mask/NaNの処理

Mask処理

Nan処理

evalの利用

dictによるファイル名の管理

awkの利用

PyGrADS

Install

Install and start pygrads

Open CTL file

Set dimentions and resolution

Get variables

Draw

小技

Caution

GmtPy

Install

GrADS to GMT

Fortran moduleの利用（f2py）

Test programe

コンパイル

pythonからの呼び出し方

どのようなFortran modulesを扱うべきか

Fortran moduleの利用（ctypes）

Test module

コンパイル

pythonからの呼び出し

pandas

Import modules

データフレームの作成

データの読み込み

データフレームの書き出し

データの追加/削除

欠損値NaNの除去

条件付きデータの抽出

for ループ

groupby

2つのデータフレームの結合

pivot table

DataFrameの有効桁数やら書式やらの表示変更方法

Draw figures

Convert to numpy format

matplotlib設定のやり方

軸の設定

ラベル名の設定

凡例の設定

文字サイズの変更

出力する画像ファイルの位置調整

pltのリセット

軸に平行な特定の線を引く。

plot(matplotlib)のまとめ

散布図

ヒストグラム

軸に平行な任意の線を引く

Original cmapの設定方法

Basemap

インストール

Import

基本設定

海岸線、国境線などの引き方