python

moduleのインストール

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

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

Anacondaのインストール

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

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

Calculate uncertanity range

nohup実行した際にプリント出力をnohup.outに書き込む方法

マニア向けの小技

awkの利用

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のインストール

Homebrewによるインストール。

$ brew install pyenv

下記、環境変数を.bash_profileに追加。

export PYENV_ROOT=~/.pyenv

export PATH=$PYENV_ROOT:$PATH

eval "$(pyenv init -)"

Anacondaのインストール

pyenvによるインストール。

$ pyenv install -l | grep anaconda

$ pyenv install anaconda3-5.0.1

$ pyenv local anaconda3-5.0.1

下記、環境変数を.bash_profileに追加。

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

テスト。ヴァージョンが帰ってくればOK。

$ 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

初期状態でipythonとするとipython3が起動するのでipythonをインストール

$ connda install ipython

py2/py3それぞれの環境下でインストール。

$ 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

(py3)$ conda install -c anaconda h5py

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

(py3)$ conda install -c anaconda h5py

(py3)$ conda install -c anaconda netcdf4

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

ファイルやディレクトリへの基本操作

ファイル名の取得

import os

os.path.basename(fDat)

フォルダ名の取得

import os

os.path.dirname(fDat)

拡張子の取得

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

import os

root,suf = os.path.splitext(fDat)

テキストデータの読み込み

f = open('test.txt', 'r')
f.read()
f.close()

リスト形式への変換

上記の方法では、読み込みんだテキストデータに改行やスペースがあり少々扱いづらい。なので、stripやsplitを用いて扱いやすようにする。

f = open('test.txt', 'r')
f.readline().strip().split(' ')
f.close()

csvなどの整備されたデータの読み込み

整備されたデータの場合はpandasで読み混んだほうがよい。

スペース区切りのテキスト・データの読み方

全て数字データの場合のみ可能

from numpy import *

r1dat = loadtxt (FDAT)

バイナリデータのIO

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

Fortranだとreal(4)書き出されたデータに対応。

import numpy as np

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

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

Fortranだとreal(8)書き出されたデータに対応。fromfileの2つ目の引数をfloat64にするだけ。簡単。この他int32などもある。

import numpy as np

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

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

一度、fromfileで読んだあとにreshapeで整形する。配列の入れ方はfortranとは逆（Y→X）なので要注意。

import numpy as np

NX, NY = 360, 180

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

三次元以上についてもreshapeすればOK。

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

reshapeの次元数の設定において、1つの次元数はわかるが、もう1つの次元数はわからない場合、下記のように-1を指定することで自動で整形してくれる。

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

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

astypeでバイナリデータの精度を指定する。tofileで書き出すのを指示する。astypeでは、'f'とかの指定もできる。

import numpy as np

Out.astype(np.float32).tofile("out.bin")

エンディアン

'f'で単精度実数、'd'で倍精度実数で読める。'>f'とすると単精度のbig endianが読める。他の読み方については下記サイトを参照。 https://qiita.com/mski_iksm/items/bb0aa375c952c2d7b91a

Mask/NaNデータの処理方法

Mask処理

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

# やり方1

import numpy as np

Dat = np.ma.masked_less(Dat,0.0).filled(-9999.0)

# やり方2

import numpy as np

Dat[np.where(Dat<0)] = -9999.0

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

Nan処理 (numpy)

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処理 (pandas)

欠損値への代入

df = df.fillna(-9999.)

欠損値データを落とす

df = df.dropna(-9999.)

2つの変数間での統計量の計算

相関係数や一次回帰式の求め方。

回帰直線

最後の整数は回帰式の次数に相当。

import numpu as np

a, b = np.polyfit(Dat1, Dat2, 1)

RegLineX = np.array([np.min(Dat1), np.max(Dat2)])

RegLineY = [a*x+b for x in RegLineX]

相関係数、P値

from scipy.stats import pearsonr

r, p = pearsonr(Model, EMDAT)

最小２乗法による重回帰式の算定

Install in conda environment

$ conda install statsmodels

Install modules

import statsmodels.api as sm

from statsmodels.sandbox.regression.predstd import wls_prediction_std

Set regression line

従属変数Yを求める説明変数を定義する。回帰式の形はコメントを参照。切片を求めるためには値1の列を挿入する必要あり(np.repeat(1,nsample))。

nsample = x.size

Set dependent variable

Y = y

# Set independent variable

# in the case of y=a+bx+cx^2

X1 = np.column_stack((np.repeat(1, nsample), x, x**2))

# in the case of y=a+bx

X2 = np.column_stack((np.repeat(1, nsample), x))

# in the case of y=ax

X3 = x

Calculate regression line by OLS

res = sm.OLS(Y, X).fit()

Check results by OLS

print res.summary()

Calculate uncertanity range

prstdの意味がよくわらからなが、iv_lとiv_uが信頼性区間の下限と上限にあたる回帰式。alphaで信頼区間を設定できる。今回は10%。

prstd, iv_l, iv_u = wls_prediction_std(res, alpha=0.09)

Mann-KendallとTeil-sen slope

import scipy.stats as sp

# Mann-Kendall

tau, prob = sp.kendalltau(x, y)

# Teil-sen slope

res = sp.mstats.theilslopes(dat)

slope = res[0] # slope

inter = res[1] # intercept

複数の変数間での相関係数

pandasだとめちゃ簡単。

df[["A", "B", "C", "C"]].corr()

ただ、p値を出したい場合は不向き。その場合はscipyなどを使って算出する。

from scipy import stats

for Var in ["B", "C", "C"]:

r, p = stats.pearsonr(df["A"].values, df[Var].values)

print ("Pb-%s: R^2 = %.3f (p-value=%.5f)" % (Var, r*r, p))

空間内挿

空間内挿その１

ctlファイルからデータを読み取り、pythonのscipy.interpolate.RegularGridInterpolaterにより、空間内挿する参考プログラム。多次元に拡張可能。

#!/usr/bin/env python

import numpy as np

from scipy.interpolate import RegularGridInterpolator

from grads import GaNum

fDatCtl = 'sst.ctl'

f1dgCtl = 'deg1x1grids.ctl' # 空間内挿後のグリッドと同じグリッドを持つctlファイル。

ga = GaNum('grads -l')

ga.open(f1dgCtl)

ga('set x 1 360')

lon1deg = ga.exp('lon').reshape(-1)

lat1deg = ga.exp('lat').reshape(-1)

ga('reinit')

grid1deg = np.array([ [y0, x0] for y0, x0 in zip(lat1deg, lon1deg)])

Out = np.zeros([12, len(lon1deg)]) # 出力したいグリッド数。

# Prepare for original grid information.

ga.open(fDatCtl)

ga('set y 1')

x = ga.exp('lon')

ga('set y 1 120')

ga('set x 1')

y = ga.exp('lat')

y[0] = -90. # 外挿できないので、端の値は-90度とする。

y[-1] = 90. # 外挿できないので、端の値は90度とする。

ga('reinit')

# Main

ga.open(fDatCtl)

for t in range(1, 13):

# Prepare interpolate function

ga('set t %s' % t)

dat = ga.exp('var')

func = RegularGridInterpolator((y, x), dat)

# Interpolate

Out[t-1, :] = func(grid1deg)

Tips

eval

dictを利用した書式管理のほうが楽。

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

辞書の追加と更新

基本のやり方。追加になる。

dicttest= {"a0":1, "a1":2}

dicttest.update({"a2":100})

→ {'a0': 1, 'a1': 2, 'a2': 100}

辞書の中に別の辞書を持つ場合は注意が必要。下記のやり方だと「更新」になってしまう。

dicttest= {"a":{"a0":1, "a1":2}}

dicttest.update({"a2":100})

→ {'a': {'a2': 100}}

「追加」したい場合はこう。

dicttest= {"a":{"a0":1, "a1":2}}

dicttest["a"].update({"a2":100})

→ {'a': {'a0': 1, 'a1': 2, 'a2': 100}}

ft00, ft24, ft48, ft72, ft96, ft 120のような書式にしたい場合、下記のように0.2iとすると良い。

A= "ft%(FT)0.2i"

グリッドデータのマスク集計

例えば、グリッドごとに算定した被害額を国ごとに集計するなど。前はFortranを噛ましてたが、numpyとpandasで簡単にできる。

import numpy as np

import pandas as pd

# DamageとMaskは一次元配列。

List = np.array([MskTmp1, DamageTmp1]).T

df = pd.DataFrame(List, columns=('ID', 'Damage'))

dfpv = pd.pivot_table(df, values=['Damage'], index=['ID'], aggfunc=np.sum).reset_index()

グリッドデータが２次元の場合は、一度.reshape(-1)として一次元化すること。また、一次元化したデータはそのままではDataFrameで取り扱えないので、４行目の"List=.."のような配列の形にしてDataFrameする。

nohup実行した際にプリント出力をnohup.outに書き込む方法

スクリプトをnohupで実行すると、計算が完了するまでnohup.outに書き込まれない問題が生じる。
これはflush outしてやれば解決する。

import sys

sys.stdout.flush()

マニア向けの小技

マニアというか信者向けの技。自分はやらないが念の為にメモっておく。

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でやる。→というかpandasでやる。

oid = 4

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

result = subprocess.run(COMMAND)