pythonで構造体的なものを

しばらく何も書いていなかったけれど、久々に投稿。
最近はpythonでコードをよく書いている。だけどまだまだ勉強し始めたばかりで、使いこなせていないのが現状だけど、だんだん使い慣れて気持ちよくコードが書けるようになってきた。

pythonの良いところはいろいろある。例えば、気に入っているのが、

• perlと比べるときれいにコードを書ける（だけどperlは大好き）
• numpy, scipy：数値計算が簡単
• ipython：変数の中身を簡単に見れたりとか
• ctypes：cとかfortranを自由に呼び出せる。okadaプログラムDC3D.f（http://www.bosai.go.jp/study/application/dc3d/DC3Dhtml_J.html）も簡単にwrapperを作れる。
• generatorとか、変な機能があってわくわくする。

きりがないけれど、そもそもpythonを使い始めた理由が数値計算にも結構使えそうだったため。なので、c、fortranで良く使う構造体が使えないかな、と思い始めた。classで実装すればいいけれど、それでは遅そうだし（未確認）。

ということでnumpy.dtypeがそれっぽそう？なので試してみる。とりあえず国土地理院（terras.gsi.go.jp）から、GEONETの座標値を落としてきて、plotさせてみた。


なんと、data['date']、data['x']のようにしてリスト？が出せる。さらに、添字が連続しているようで、data[i][0]、data[i][1]のようにも使える。これはかなり良いかもしれない。速度はどれくらい違うんだろう。numpy、すごいなあ。

import numpy as np
import datetime as dt
import matplotlib.pyplot as plt

def main():
    posf = "950255.11.pos"

    # structure?
    pos_t = np.dtype([('date', dt.datetime),
                      ('x', np.float), 
                      ('y', np.float), 
                      ('z', np.float), 
                      ('lat', np.float), 
                      ('lon', np.float), 
                      ('h', np.float), 
                      ])

    with open(posf, 'rt') as f:
        for line in f:
            if line[1:6] == "EPOCH":
                nd = int(line[72:75]) # number of data
            if line[0:5] == "+DATA":
                break

        if nd > 0:
            data = np.ndarray([nd], dtype=pos_t)

        i = 0
        for line in f:
            if line[0:1] == "*":
                continue
            elif line[0:5] == "-DATA":
                break
            else:
                date_string    = line[1:20]
                x1,   y1,   z1 = line[21:38], line[39:56], line[57:74]
                lat1, lon1, h1 = line[75:92], line[93:110], line[111:128]

                #data[i][0] = dt.datetime.strptime(date_string, "%Y %m %d %H:%M:%S")
                #data[i][1], data[i][2], data[i][3] = np.float(x1), np.float(y1), np.float(z1)
                #data[i][4], data[i][5], data[i][6] = np.float(lat1), np.float(lon1), np.float(h1)
                data[i]['date'] = dt.datetime.strptime(date_string, "%Y %m %d %H:%M:%S")
                data[i]['x'], data[i]['y'], data[i]['z'] = np.float(x1), np.float(y1), np.float(z1)
                data[i]['lat'], data[i]['lon'], data[i]['h'] = np.float(lat1), np.float(lon1), np.float(h1)

                i+=1

    plt.plot_date(data['date'], data['x'])
    plt.savefig("test.png")
    plt.show()

    return data

if __name__ == '__main__':
    data = main()

ところで、プロットさせてみた座標値（x）がこちら。東日本大震災で大きく観測点が移動している。x座標値で10cmくらい。石川県なんだけどな、この観測点・・・。