Reading a binary file in Python using a control file with multiple variables
I am trying to read a binary file in Python using a Fortran code as an example.
The binary file is called data.grads and I have a control file called data.ctl that should allow me to understand how to read the binary file. As I said, I have a Fortran code which my advisor wrote to explain me the process of reading the binary file and constructing arrays corresponding to the different variables inside the control file (speed, temperature, pressure, etc).
I read multiple posts on the matter and I am having some trouble understanding how to use the data inside the binary file.
Here are some posts that I checked:
- Reading a direct access fortran unformatted file in Python
- Reading direct access binary file format in Python
- Reading records from a Fortran binary file in Python
- Reading Fortran binary file in Python
- Reading fortran direct access data and writing formatted data - faster with python than with fortran?
- Python reading unformatted direct access Fortran 90 gives incorrect output
The binary file stores simulations results of different properties of the atmosphere that scientists use to plot, for instance, temperature and pressure sections such as these:
I was told that the control file is key to understanding how to get anything out of the binary file.
I am able to read the file, but I don’t know how to access results for a specific variable. Here is a piece of code that I was able to derive from some of the posts that I quoted before:
filename = "/path/data.grads"
nlat = 32
nlon = 67
f = open(filename, 'rb')
recl = np.fromfile(f, dtype='int32', count=4*nlat*nlon)
f.seek(4)
field = np.fromfile(f, dtype='float32',count=-1)
print('Record length=',recl)
print(field, len(field))
It returns the following:
Record length= [-134855229 -134855229 -134855229 ... -134855229 -134855229 -134855229]
[-9.99e+33 -9.99e+33 -9.99e+33 ... -9.99e+33 -9.99e+33 -9.99e+33] 10319462399
Can someone who has experience in the subject help me figure out how to access the different variables using the control file ?
If you need more explanations, please tell me, I will edit my posts and add as much information as you want.
Unfortunately, I cannot share the binary file because it is on a server and weighs about 40 GB…
I share:
- the content of the control file,
- and the Fortran code (it lacks some code because it was written as an explanation).
Control file (data.ctl)
DSET ^data.grads
UNDEF -9.99e33
XDEF 64 LINEAR 0.0 5.6250
YDEF 32 LEVELS
-85.761 -80.269 -74.745 -69.213 -63.679 -58.143 -52.607 -47.070 -41.532
-35.995 -30.458 -24.920 -19.382 -13.844 -8.307 -2.769 2.769 8.307
13.844 19.382 24.920 30.458 35.995 41.532 47.070 52.607 58.143
63.679 69.213 74.745 80.269 85.761
ZDEF 67 LEVELS
.000696 .08558 .1705 .2554 .3402 .4544 .6274 .8599 1.152 1.505 1.918 2.392
2.928 3.495 4.063 4.781 5.816 7.181 8.889 10.95 13.39 16.07 18.81 21.61
24.47 27.39 30.37 33.40 36.50 39.64 42.77 45.90 49.04 52.17 55.31 58.44
61.57 64.71 67.84 70.98 74.11 77.24 80.38 83.51 86.65 89.78 92.92 96.05
99.18 102.3 105.5 108.6 111.7 114.9 118.0 121.1 124.3 127.4 130.5 133.7
136.8 139.9 143.1 146.2 149.3 152.5 160.0
TDEF 3120 LINEAR 01JAN2000 1HR
VARS 31
u 67 99 u (m/s)
v 67 99 v (m/s)
w 67 99 w (m/s)
T 67 99 T (K)
dia 67 99 diagnostics (see table)
ps 0 99 ps
Ts 0 99 Ts (K)
h2og 67 99 Water vapor [kg/m^2]
h2oim1 67 99 Water ice mass for dust 0.30E-07 [kg/m^2]
h2oim2 67 99 Water ice mass for dust 0.10E-06 [kg/m^2]
h2oim3 67 99 Water ice mass for dust 0.30E-06 [kg/m^2]
h2oim4 67 99 Water ice mass for dust 0.10E-05 [kg/m^2]
h2oin1 67 99 Water ice number for dust 0.30E-07 [number/m^2]
h2oin2 67 99 Water ice number for dust 0.10E-06 [number/m^2]
h2oin3 67 99 Water ice number for dust 0.30E-06 [number/m^2]
h2oin4 67 99 Water ice number for dust 0.10E-05 [number/m^2]
h2ois 0 99 surface h2o ice [kg/m^2]
p 67 99 Pressure [Pa]
h 67 99 Height above the surface [m]
dipre 67 99 Delta pressure [Pa]
surf 67 99 space of cell factor [sin*cos]
dm 67 99 CO2 ice cloud mass concentration [kg/m^3]
cap 0 99 Surface CO2 ice mass [kg]
hcap 0 99 Surface CO2 ice [kg/m^2]
gdq_op 0 99 dust from rad.mod [opacity]
gdq_mix 67 99 dust from rad.mod [mix.r]
dust_op 0 99 dust from tracers [opacity]
dust_n1 67 99 Dust num dens from tracers for R=0.30E-07 [m^-3]
dust_n2 67 99 Dust num dens from tracers for R=0.10E-06 [m^-3]
dust_n3 67 99 Dust num dens from tracers for R=0.30E-06 [m^-3]
dust_n4 67 99 Dust num dens from tracers for R=0.10E-05 [m^-3]
ENDVARS
Fortran file
open(12,file='data.grads',status='unknown',
& form='unformatted',access='direct',
& recl = 4*nlat*nlon )
krec=0
do l = 1,nlat
krec = krec+1
write(12,rec=krec,ERR=900)real(u3d(1:nlon,1:nlat,l))
enddo
do l = 1,nlat
krec = krec+1
write(12,rec=krec,ERR=900)real(v3d(:,:,l))
enddo
do l = 1,nlat
krec = krec+1
write(12,rec=krec,ERR=900)real(w3d(:,:,l))
enddo
do l = 1,nlat
krec = krec+1
write(12,rec=krec,ERR=900)real(T3d(:,:,l))
enddo
do l = 1,nlat
krec = krec+1
write(12,rec=krec,ERR=900)real(D3(:,:,l))
enddo
krec = krec+1
write(12,rec=krec,ERR=900)real(ps3d(:,:))
krec = krec+1
write(12,rec=krec,ERR=900)real(Ts3d(:,:))
do n = 1,4
do l = 1,nlat
krec = krec+1
write(12,rec=krec,ERR=900)real(trace4D(:,:,l,n))
enddo
enddo
do n = 1,4
krec = krec+1
write(12,rec=krec,ERR=900)real(trace2D(:,:,n))
enddo
do l = nlat,1,-1
krec = krec+1
write(12,rec=krec,ERR=900)real(pre3d(:,:,l))
enddo
do l = 1,nlat
krec = krec+1
write(12,rec=krec,ERR=900)real(alth3d(:,:,l))
enddo
do l = 1,nlat
krec = krec+1
write(12,rec=krec,ERR=900)real(dipre3d(:,:,l))
enddo
do l = 1,nlat
krec = krec+1
write(12,rec=krec,ERR=900)real(surf3d(:,:,l))
enddo
do l = 1,nlat
krec = krec+1
! Convert DM to DM/grid volume
write(12,rec=krec,ERR=900) real(DM4(:,latmask(:),l)
& *GRAV*pre3d(:,:,nlat-l+1)
& /(T3d(:,:,nlat-l+1)*dipre3d(:,:,nlat-l+1)*RGAS)
& /surf3d(:,:,nlat-l+1) )
enddo
krec = krec+1
write(12,rec=krec,ERR=900)real(CAP4(:,latmask(:)))
krec = krec+1
write(12,rec=krec,ERR=900)real(HCAP4(:,latmask(:)))
krec = krec+1
write(12,rec=krec,ERR=900)real(DDUSTA3d(:,:))
do l = 1,nlat
krec = krec+1
write(12,rec=krec,ERR=900)real(GDQ3d(:,:,l))
enddo
krec = krec+1
write(12,rec=krec,ERR=900)real(tau_dust3d(:,:))
do n = 1,4
do l = 1,nlat
krec = krec+1
write(12,rec=krec,ERR=900)real(dust_n3d(:,:,l,n))
enddo
enddo
EDIT: First lines of the binary file (xxd -l 100 data.grads)
0000000: c345 f6f7 c345 f6f7 c345 f6f7 c345 f6f7 .E...E...E...E..
0000010: c345 f6f7 c345 f6f7 c345 f6f7 c345 f6f7 .E...E...E...E..
0000020: c345 f6f7 c345 f6f7 c345 f6f7 c345 f6f7 .E...E...E...E..
0000030: c345 f6f7 c345 f6f7 c345 f6f7 c345 f6f7 .E...E...E...E..
0000040: c345 f6f7 c345 f6f7 c345 f6f7 c345 f6f7 .E...E...E...E..
0000050: c345 f6f7 c345 f6f7 c345 f6f7 c345 f6f7 .E...E...E...E..
0000060: c345 f6f7 .E..
The numpy equivalent (or rather counterpart, since you showed the writer) of that Fortran begins
import numpy
nlat=32 # from data.ctl
nlon=64
nz=64 # or 67?
dt=numpy.dtype((numpy.float32,(nlat,nlon)))
def read(n=nz): return numpy.fromfile(f,dt,n)
def read1(): return read(1)[0]
with open("data.grads",'rb') as f:
u3d=read()
v3d=read()
# …
ps3d=read1()
# …
trace4D=numpy.fromfile(f,numpy.dtype((dt,nlat)),4)
trace2D=read(4)
pre3d=read()[::-1]
# …
There is some confusion between the nlat and nz variables; since you said the code was altered, the data.ctl might be the better reference.