我是Python的新手,我正在从头开始用Python编写可视化代码(以避免使用昂贵的专有程序,如IDL).到目前为止,我已经使用过IDL和gnuplot.我想要做的是:
我使用fortran将二维数组写入未格式化的直接访问文件,我希望能够在python中读取.确切的测试代码如下.实际代码是一个巨大的并行代码,但数据输出几乎完全相同的格式.
program binary_out
implicit none
integer :: i,j,t,rec_array
double precision, dimension(100,100) :: fn
double precision, parameter :: p=2*3.1415929
INQUIRE(IOLENGTH=rec_array) fn
open(unit=10,file='test',status='new',form='unformatted',access='direct',recl=rec_array)
fn=0
write(10,rec=1) fn
do t=1,3
do i=1,100
do j=1,100
fn(i,j)=sin(i*p*t/100)*cos(j*p*t/100)
enddo
enddo
write(10,rec=t+1) fn
enddo
close(10)
end program binary_out
该程序应该给我t = 1的零和增加"岛"的数量以增加t的值.但是当我使用下面给出的python代码读取它时,我只是得到了零.如果我删除第一个零的写语句,我只是得到第一次切片,而不管我在python代码中使用的"timeslice"的值.我到目前为止的代码是:
#!/usr/bin/env python
import scipy
import glob
import numpy as np
import matplotlib.pyplot as plt
import os, sys
from pylab import *
def readslice(inputfilename,field,nx,ny,timeslice):
f=open(inputfilename,'r')
f.seek(timeslice*nx*ny)
field=np.fromfile(inputfilename,dtype='d',count=nx*ny)
field=np.reshape(field,(nx,ny))
return field
a=np.dtype('d')
a=readslice('test',a,100,100,2)
im=plt.imshow(a)
plt.show()
如果时间片等于i,我希望def readlice能够在第i个位置读取记录.为此,我试图使用f.seek,但它似乎不起作用.numpy.fromfile似乎开始从第一个记录本身读取.如何让numpy.fromfile从文件中的特定点读取？
我仍然在尝试习惯Python风格并深入研究文档.任何帮助和指示将不胜感激.
我有FORTRAN代码的二进制输出文件.想在python中阅读它.(使用FORTRAN读取并输出文本以读取python不是一个选项.长篇故事.)我可以以简单的方式阅读第一条记录:
>>> binfile=open('myfile','rb')
>>> pad1=struct.unpack('i',binfile.read(4))[0]
>>> ver=struct.unpack('d',binfile.read(8))[0]
>>> pad2=struct.unpack('i',binfile.read(4))[0]
>>> pad1,ver,pad2
(8,3.13,8)
正好.但这是一个大文件,我需要更有效地做到这一点.所以我尝试:
>>> (pad1,ver,pad2)=struct.unpack('idi',binfile.read(16))
这不会运行.给我一个错误并告诉我解压缩需要一个长度为20的参数.这对我来说没有任何意义,因为我上次检查时,4 + 8 + 4 = 16.当我输入并用16替换16时,它会运行,但这三个数字都填充了数字垃圾.有谁看到我做错了什么？谢谢!
以下是数据的编写方式(浮点数的2-D矩阵.我不确定大小).
open(unit=51,file='rmsd/'//nn_output,form='unformatted',access='direct',status='replace',&
recl=Npoints*sizeofreal)
!a bunch of code omitted
write(51,rec=idx-nstart+1) (real(dist(jdx)),jdx=1,Npoints)
f = open(inputfilename,'rb')
field = np.fromfile(f,dtype='float64')
但结果不正确.矩阵的对角线应为0(或非常接近),因为它是相似矩阵.我尝试了不同dtype的但我仍然无法得到正确的结果.
编辑:这是我得到的输出
array([ 1.17610188e+01, 2.45736970e+02, 7.79741823e+02, ...,
9.52930627e+05, 8.93743127e+05, 7.64186127e+05])
我还没有为重塑而烦恼,因为值应该在0到20之间.
编辑2:这是指向文件前100行的链接:https: //drive.google.com/file/d/0B2Mz7CoRS5g5SmRxTUg5X19saGs/view？usp =sharing
编辑3:文件顶部的声明
integer,parameter :: real_kind=8
!valuables for MPI
integer :: comm, nproc, myid, ierr
integer,allocatable :: idneigh(:)
real :: tmp
real,allocatable :: traj(:,:)
real(real_kind),allocatable :: dist(:)
real(real_kind),allocatable :: xx(:,:),yy(:,:),rot(:),weight(:)
character(200) :: nn_traj, nn_output, nn_neigh
integer,parameter :: sizeofreal=4,sizeofinteger=4