ベクトル図

数値モデルの計算から得られた流速をベクトル図として描くのに、GMT*1 を使って描こうと思った。

0 0 90 10
0 1 45 15
1 0 90 15
1 1 135 20

みたいな “x座標 y座標 角度 大きさ” が並んだテキストファイルを “psxy” コマンドに渡して、PostScript ファイルを得る。そのほかタイトルとかを付けようと思ったらさらにその PostScript ファイルに GMT のコマンド群で追記していく。

数値モデルの結果は netCDF ファイルだから、これを適当にテキストファイルにして、それをコマンドに渡していくのだけど、いちいちそれをするのは非常に不便。それで、Python の netCDF4 ライブラリを使って一時的なテキストファイルを作って GMT コマンドを呼び出すなどしてみたんだけど、割と不快な感じだった。

GMT を使う必要がないような気がしてきて、Python の matplotlib を使うことにするとかなり快適。テキストファイルにするとか余計なことしなくて良い。netCDF から読み取った Numpy の array をそのまま渡すだけでプロットしてくれる。

適当に書いたコード

# -*- coding: utf-8 -*-

from netCDF4 import Dataset
import matplotlib.pylab as plt

filename = 'results/pom2k.nc'

# read netCDF data
root = Dataset(filename, 'r')

times = root.variables[u'time'][:].copy()
depths = root.variables[u'z'][:].copy()

u_values = root.variables[u'u'][:, :, :, :].copy()
v_values = root.variables[u'v'][:, :, :, :].copy()

root.close()


# plot vector diagram
def plot_vector(time, depth):
    
    dist = 3
    
    u = u_values[time, depth, :, :]
    v = v_values[time, depth, :, :]
    
    plt.figure(figsize=(8, 6), dpi=96)
    
    x, y = plt.meshgrid(plt.arange(600), plt.arange(400))
    
    quiver = plt.quiver(x[::dist * 2, ::dist * 2], y[::dist * 2, ::dist * 2],
                        u[::dist, ::dist], v[::dist, ::dist], scale=10)
    plt.quiverkey(quiver, 0.85, 0.02, 1, '1 m/s', coordinates='figure')
    
    plt.axis([0, 600, 0, 400])
    
    plt.title('Velocity ({t} days, {d} sigma depth)'.format(
              t=times[time], d=depths[depth]))
    
    plt.savefig('velocity-{t}days-{d}depth.png'.format(
                t=times[time], d=depths[depth]))

plot_vector(time=10, depth=0)

GMT ももう少しでリリースされそうな次の v5 は API が整備されてるし、LL のバインディングとか使えるようになりそう。そしたら便利。

追記（2012/05/24）

これ書いてから完全に忘れてたんですけど、結局 GMT は使わず matplotlib 使いました。

*1:Generic Mapping Tools