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

Question

前提・実現したいこと

数値計算データを可視化するために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要素あたり4byte×202*481*3=4×291486=1165944byte

# %% 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=imjm3として、uvhpyを３次元配列とreshapeして描かせたら、グラフらしきものは出てきました。

Answer 1

多分 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))
でも同じ事だと思う。

Answer 2

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

Comments

[YoshitoIshii]

まず、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

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

[y_waiwai]

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

Answer 3

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')
で正しかったのでは？

Comments

[YoshitoIshii]

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

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

[mkgrei]

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

[YoshitoIshii]

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

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

[mkgrei]

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

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

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