質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

90.35%

  • Python

    9163questions

    Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

  • NumPy

    508questions

    NumPyはPythonのプログラミング言語の科学的と数学的なコンピューティングに関する拡張モジュールです。

  • FORTRAN

    61questions

    FORTRAN(フォートラン)は科学時術計算に向いた手続き型プログラミング言語です。 並列計算の最適化が行いやすい特性上、数値予報および気候モデルなどの大規模な計算を行う分野のスーパーコンピュータで使われています。

Fortranのバイナリ出力データ(unformatted)をPythonで読む方法が分かりません

受付中

回答 3

投稿 編集

  • 評価
  • クリップ 0
  • VIEW 667

YoshitoIshii

score 2

 前提・実現したいこと

数値計算データを可視化するためにanacondaで描画プログラムを書いています。
以前は地球流体電脳ライブラリ(DCL)を用いて、fortranで計算したデータから以下のようなグラフを可視化できておりました。イメージ説明
今回はfortranで計算したデータをpythonを用いて読み込ませようとしているのですが、配列のデータの切り方が間違っているため、正しい描画ができません。イメージ説明

fortranで一定ステップ計算した後、下記のサブルーチンから(200回程度)出力させているのですが、20回分のデータをpythonで読み込み作図しようとしているがうまくいきません。
pythonのソースコードは以下のリンクを参照して作りました。
http://ig.hateblo.jp/entry/2014/05/30/225607

 該当のソースコード

  subroutine writefile(u,v,h)
    implicit none
    real(8),dimension(1:im ,1:jm ,1:3 ),intent(in)::u,v
    real(8),dimension(1:im ,0:jm ,1:3 ),intent(in)::h
    real(4),dimension(1:im ,1:jm ,1:3 ) ::uvh
    integer::i,j,m

      open(11,file='uvh.dat', status='unknown', form='unformatted')

      do j=1,jm
        do i=1,im
          uvh(i,j,1)=real(u(i,j,2))
          uvh(i,j,2)=real(v(i,j,2))
          uvh(i,j,3)=real(h(i,j,2))
        enddo
      enddo
      write(11) (((uvh(i,j,m),i=1,im),j=1,jm),m=1,3)
    return
  end subroutine writefile


*pythonのデータ読み取り箇所

import numpy as np

im = 202
jm = 481
nstep = 20

# %% 0.fromfileを用いたファイル読み込み
f90 = open('uvh.dat', 'rb')
ary = np.fromfile(f90, np.float32,count=im*jm*3*20)

# %% 1.読み込んだn番目の要素を取り出し、fortran形式の4次元配列に
uvhpy = ary.reshape(jm,im,3,nstep, order='F')  # 1要素あたりbyte×202*481*3=4×291486=165944byte

# %% 2.コンターを作図する
import matplotlib.pyplot as plt

x = np.linspace(0,50 , im)
y = np.linspace(-60, 60, jm)
X, Y = np.meshgrid(x, y)

for i in range(nstep):
    Z = uvhpy[:, :,2,i]
    plt.contour(X,Y,Z)
    plt.show()

 試したこと

pythoでのforループをなしにして、ary内のcount=im*jm*3として、uvhpyを3次元配列とreshapeして描かせたら、グラフらしきものは出てきました。

  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

回答 3

0

program test
  implicit none
  real(4),dimension(1:10, 1:10, 1:3) :: u=0d0
  integer :: i,j,k,m
  open(11, file='test.dat', form='unformatted')
  do m=1,20
    do i=1,10
      do j=1,10
        u(i,j,1) = 1
        u(i,j,2) = 2
        u(i,j,3) = 3
      enddo
    enddo
    write(11) (((u(i,j,k), i=1,10), j=1,10), k=1,3)
  end do
  return
end program test
import numpy as np
from scipy.io import FortranFile

im = 10
jm = 10
nstep = 10

f90 = open('test.dat', 'rb')
ary = np.fromfile(f90, np.float32, count=nstep*im*jm*3)

uvhpy = ary.reshape(im, jm, 3, 10, order='F')  # byte×202*481*3=4×291486=1165944
print(uvhpy[:, :, 2, 0])

#----------OR----------
#uvhpy = ary.reshape(10, 3, jm, im)  # byte×202*481*3=4×291486=1165944
#print(uvhpy[0, 2, :, :])

事前に.transpose(0,1,3,2)を入れてみてはいかがですか? for文で抽出する際のnumpy.arrayの次元を疑っています。 shapeが(n,1)になっていませんか?


手元でやってみたところ、再現しませんでした。
uvhpy = ary.reshape(jm,im,3,20)#, order='F')
ではなく、
uvhpy = ary.reshape(jm,im,3,20, order='F')
で正しかったのでは?

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2018/05/19 13:34

    編集した「前提・実現したいこと」に、ご指摘いただいた
    uvhpy = ary.reshape(jm,im,3,20, order='F')
    を基にしたグラフの1枚目を掲載しました。

    uvhpyの4つ目の要素は、fortranからの数値計算を何回読み込むかを意図しているつもりです。しかし複数のグラフをループで描かせたい場合は、配列の次元数を無闇に増やさず、別のデータの切り出し方を考えたほうが賢明でしょうか?

    キャンセル

  • 2018/05/19 13:58

    配列の並び方が異なるので、20回だけ欲しい場合でも、書き出し分を全部読み込んでみてください。

    キャンセル

  • 2018/05/20 16:43 編集

    了解しました。変数aryで全てのデータを読み込みさせます。

    追加で教えていただきたい点を1つだけリクエストさせてください。pythonでのreshape内部の変数を複数設定した場合(例えば(jm,im,m,nstep))、読まれるのはカッコ内の先頭からループ順番からでしょうか?
    現状の問題に即すると、fortran側から計算結果を書き出したのは、外側からm,j,iの順で3重ループにしたと思っています。これをnstep回行ったと理解しています。numpy側から1次元データとして読み、データ列からreshapeを通じて配列を復元する場合、(nstep,m,j,i)と順序を変更したらうまくいくのではと考えました。
    (この変数の順序ではreshapeできないよ、というエラーが出てきてワークしない結果を得ました。またfortran側の配列データにない、変数nstepをもって配列の次元数を変えてしまっているのも正確な配列復元の邪魔をしているのではないかなと思います。)

    キャンセル

  • 2018/05/24 06:56

    私が少し混乱していたようです。

    上のコードを2つ追記しました。
    これで挙動を確認していただければわかると思います。

    (nstep,3,j,i)か(i,j,3,nstep)+order='F'かですね。
    読み込まれた後、軸の順番が異なることに気をつけてください。

    キャンセル

0

あなたのお使いのFortrunの浮動小数点形式は4バイトなんでしょうか。
そして、その形式はPythonのそれと同一なんでしょうか。
まずそこから確認する必要があると思いますが。

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2018/05/19 12:51

    まず、fortranのソースコードをコンパイルする際にインテルfortranを用いています。書き出しているデータは単精度実数浮動小数点値で、データ型はrea(4)、記憶域は4バイトです。
    http://wwweic.eri.u-tokyo.ac.jp/computer/manual/altix/compile/Fortran/Intel_Fdoc100/main_for/mergedProjects/bldaps_for/common/bldaps_datarepov.htm

    次に、pythonでのnumpyを用いた読み込みでは、データ型はnp.float32を明示しています。これは32bitで単精度の浮動小数点数として同一の対応をさせていると思います。
    https://deepage.net/features/numpy-dtype.html

    当方、基本がよく分かっていないので的外れな回答になっているかもしれません。ご不明な点があれば都度いただければ幸甚です。

    キャンセル

  • 2018/05/19 22:27

    なら、単純な浮動小数点数をいくつか並べた1次元配列をバイナリデータとして保存し、それをpythonで読み込んで期待する結果になるか確認することですね。
    それで期待する結果になるのを確認したあと、それを2次元配列で同様にやってみる、ってことが必要なんじゃないでしょうか

    キャンセル

0

多分 open(11, 云々 ,access='stream') でデータべた書きにしないと駄目なのでは?
unformatted だけだとバイナリで書き出すがレコード区切りが入る。 stream は Fortran2003 以降の属性だが、昔風には direct access 形式でもまぁ普通は行けなくもない。

あと、DO LOOP を書かずに全配列操作を使えばいいのでは?
open(11,file='uvh.dat', status='unknown', form='unformatted, access='stream')
uvh(:,:,1)=real(u(:,:,2))
uvh(:,:,2)=real(v(:,:,2))
uvh(:,:,3)=real(h(:,:,2))
write(11) uvh
これだけでいいはず。
もしくは、いちいち代入せずに
open(11,file='uvh.dat', status='unknown', form='unformatted, access='stream')
write(11) real(u(:,:,2))
write(11) real(v(:,:,2))
write(11) real(h(:,:,2))
でも同じ事だと思う。

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 90.35%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

同じタグがついた質問を見る

  • Python

    9163questions

    Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

  • NumPy

    508questions

    NumPyはPythonのプログラミング言語の科学的と数学的なコンピューティングに関する拡張モジュールです。

  • FORTRAN

    61questions

    FORTRAN(フォートラン)は科学時術計算に向いた手続き型プログラミング言語です。 並列計算の最適化が行いやすい特性上、数値予報および気候モデルなどの大規模な計算を行う分野のスーパーコンピュータで使われています。