pygmt

はじめに

PyGMTはできることとできないこととがある。その点を意識しないとドツボにはまる。

できること/できないこと

• できること

  • • surface

    • plot

    • pscontour

    • 地形図：かなり簡単。pygmt.datasets.load_earch_reliefで好きな解像度の標高をダウンロードできる。また、取得したデータはxarray形式なので、topo.valuesとすればnumpy形式で取得可能。

    • 複数の図を並べる：xshiftで対応可能。参考になる資料１の中段を参考に。

    • 画像の保存（eps, png, tif）：ダイレクトにできるので楽。

• できないこと

  • • xyz2grd：なんとpygmtに導入されてない。どうしてもやりたいないら、外部のコマンドとして呼び出す。

    • grdvector：矢印を描画できない。外部から呼び出すのも大変なので、矢印をかく場合ははじめからシェルスクリプトで対応する。

参考になる資料

1. https://www.leouieda.com/posters/agu2019.html

モジュールのインストール方法(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で統一するようにする。

1. Pyenvのインストール

Homebrewによるインストール。

$ brew install pyenv

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

export PYENV_ROOT=~/.pyenv

export PATH=$PYENV_ROOT:$PATH

eval "$(pyenv init -)"

2. 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

3. 仮想環境の構築

$ 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

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)

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処理

１次元配列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処理

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

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

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)