python
python
$ sudo apt-get install python-numpy
$ sudo apt-get install python-matplotlib
$ sudo apt-get install python-mpltoolkits.basemap
$ pip3 install pylabモジュールのインストール方法(Anaconda)
$ conda install pylab
念のため、インストールできているかの確認。ヴァージョンが確認できればOK!!
import numpy
print (numpy.__version__)
使った感じ、condaでインストールしていくのが良いだろう （企業が提供しているという点が気になるが）。 pythonは様々なバージョンがあり、 pipだと環境を破壊してしまう場合があるので、 condaで統一するようにする。
datetimeが簡単で良い。
import datetime
# 日付変数の定義
Date = datetime.datetime(2024,1,1,6)
# 日付の計算
Date2 = Date + datetime.timedelta(days=1)
Date3 = Date + datetime.timedelta(hours=6)
# grads形式の変換
Date.strftime("%HZ%d%h%Y")
# 10桁数字への変換
Date.strftime("%Y%m%d%H00")
datetime.strptime(Date, '%Y%m%d%H00')
下記のようにすれば、引数で受け取った文字列の日付を、任意の文字列に変換できる。
datetime.datetime.strptime(Date, '%Y%m%d%H00').strftime("%HZ %d%h%Y")
Homebrewによるインストール。
$ brew install pyenv
下記、環境変数を.bash_profileに追加。
export PYENV_ROOT=~/.pyenv
export PATH=$PYENV_ROOT:$PATH
eval "$(pyenv init -)"
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()
整備されたデータの場合はpandasで読み混んだほうがよい。
全て数字データの場合のみ可能
from numpy import *
r1dat = loadtxt (FDAT)
Fortranだとreal(4)書き出されたデータに対応。
import numpy as np
Dat = np.fromfile("test.bin", np.float32)
Fortranだとreal(8)書き出されたデータに対応。fromfileの2つ目の引数をfloat64にするだけ。簡単。この他int32などもある。
import numpy as np
Dat = np.fromfile("test.bin", np.float64)
一度、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
１次元配列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の方法で統一するようにする（メモリは食うが）。
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)
欠損値への代入
df = df.fillna(-9999.)
欠損値データを落とす
df = df.dropna(-9999.)
相関係数や一次回帰式の求め方。
最後の整数は回帰式の次数に相当。
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]
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
従属変数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()
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)
dictを利用した書式管理のほうが楽。
GCM01 = ["MIROC5", 1]
NUM = "01"
a = eval ("GCM" + NUM + "[1]")
ファイルやディレクトリの管理が劇的にしやすくなる。
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に書き込まれない問題が生じる。
これはflush outしてやれば解決する。
import sys
sys.stdout.flush()
マニアというか信者向けの技。自分はやらないが念の為にメモっておく。
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)