How to convert a unformatted fortran file(modflow output ) to numpy array

Question

I have a modflow output file with extension hds. Google drive link fora file. Its a unformatted fortran file. I need to convert it to numpy array, I have tried :

floattype = 'f4'
a = np.fromfile("lake_example.hds", np.dtype([('kstp','i4'),('kper','i4'),('pertim',floattype),('totim',floattype),('text','a16'),('ncol','i4'),('nrow','i4'),('ilay','i4')]))
print a
print a.shape

github link for code: https://github.com/Kirubaharan/hydrology/blob/master/gw_tut.py

I am a trying a tutorial from this link. Since I am on linux I cannot use flopy's method to get the output array from the file. So I am trying to use np.fromfile, But I am having problem in getting the output.

My output is like this now:

[ (44, 1, 1.401298464324817e-45, 1.0, '\x00\x00\x80?            ', 1145128264, 11, 11)
 (1, 44, 6.782284567332115e-43, 100.0, '\x00\x00\xc8B\x00\x00\xc8B\x00\x00\xc8B\x00\x00\xc8B', 1120403456, 1120403456, 1120403456)
 (1120403456, 1120403456, 100.0, 100.0, '\x0c\xbf\xc7B\x18~\xc7B=@\xc7B\xce\x0e\xc7B', 1120336356, 1120341710, 1120354365)
 (1120370200, 1120386828, 100.0, 100.0, '\x18~\xc7B\x0e\xf9\xc6B\xf0s\xc6B\xaa\x00\xc6B', 1120258308, 1120272554, 1120302064)
 (1120336142, 1120370200, 100.0, 100.0, '=@\xc7B\xf0s\xc6B\xf8\x94\xc5B\x91\xb3\xc4B', 1120149448, 1120187281, 1120244984)
 (1120302064, 1120354365, 100.0, 100.0, '\xce\x0e\xc7B\xaa\x00\xc6B\x91\xb3\xc4B\xac\xff\xc2B', 1119940155, 1120075692, 1120187281)
 (1120272554, 1120341710, 100.0, 100.0, '\xe4\xf9\xc6B\x04\xc9\xc5B\xc8\x1f\xc4B;\xee\xc0B', 1119092736, 1119940155, 1120149448)
 (1120258308, 1120336356, 100.0, 100.0, '\xce\x0e\xc7B\xaa\x00\xc6B\x91\xb3\xc4B\xac\xff\xc2B', 1119940155, 1120075692, 1120187281)
 (1120272554, 1120341710, 100.0, 100.0, '=@\xc7B\xf0s\xc6B\xf8\x94\xc5B\x91\xb3\xc4B', 1120149448, 1120187281, 1120244984)
 (1120302064, 1120354365, 100.0, 100.0, '\x18~\xc7B\x0e\xf9\xc6B\xf0s\xc6B\xaa\x00\xc6B', 1120258308, 1120272554, 1120302064)
 (1120336142, 1120370200, 100.0, 100.0, '\x0c\xbf\xc7B\x18~\xc7B=@\xc7B\xce\x0e\xc7B', 1120336356, 1120341710, 1120354365)

I have included only few lines of output.

For header information you can refer their source code:https://github.com/modflowpy/flopy/blob/master/flopy/utils/binaryfile.py#L30g

Thomas Ritschel · Accepted Answer · 2015-02-27T14:59:25.497

Your code doesn't match the structure of your data file:

00000000  2c 00 00 00 01 00 00 00  01 00 00 00 00 00 80 3f  |,..............?|
00000010  00 00 80 3f 20 20 20 20  20 20 20 20 20 20 20 20  |...?            |
00000020  48 45 41 44 0b 00 00 00  0b 00 00 00 01 00 00 00  |HEAD............|
00000030  2c 00 00 00 e4 01 00 00  00 00 c8 42 00 00 c8 42  |,..........B...B|
00000040  00 00 c8 42 00 00 c8 42  00 00 c8 42 00 00 c8 42  |...B...B...B...B|
00000050  00 00 c8 42 00 00 c8 42  00 00 c8 42 00 00 c8 42  |...B...B...B...B|

Each data block has it's own 56 bytes header consisting of: 3 integers (i4), 2 floating point values (f4), 16 characters, and again 5 integers (i4):

44 1 1
1.0 1.0
            HEAD
11 11 1 44 484

Then the data block follows (11x11 floating point values):

100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0
100.0 99.87313842773438 99.74627685546875 99.6254653930664 ...

I'm not sure whether this could be imported directly into a numpy array.

The following sample code will loop over the whole file and extracts header and data for each block:

#!/usr/bin/python

import struct
import numpy as np

infile = open("lake_example.hds","rb")

blockdata = []

while infile.read(1):
    infile.seek(-1,1)
    data = infile.read(56)
    n = struct.unpack('<3i4', data[0:12])
#    print n[0], n[1], n[2]
    n = struct.unpack('<2f4', data[12:20])
#    print n[0], n[1]
#    print data[20:36]
    n = struct.unpack('<5i4', data[36:56])
#    print n[0], n[1], n[2], n[3], n[4]
    ncol = n[0]
    nrow = n[1]
    a = np.fromfile(infile,dtype='f4',count=ncol*nrow).reshape((ncol,nrow))
    blockdata.append(a)
    data = infile.read(4)
    n = struct.unpack('<i4', data)
#    print n[0]

for block in blockdata:
    print block

Most likely you also will need some of the information from the block headers (see print statements).

I need to convert the second data block(11x11 floating point values) — Kirubaharan J, Feb 27 '15 at 13:11
Can you please tell how to extract different data blocks in python — Kirubaharan J, Feb 27 '15 at 13:13
I have changed my example code so that it now reads all data blocks from the file. — Thomas Ritschel, Feb 27 '15 at 15:01

score 2 · Answer 2 · answered May 11 '15 at 01:53

2

See also the 'flopy.utils.binaryfile Module': http://modflowpy.github.io/flopydoc/binaryfile.html

Look at Flopy-3 Tutorial 2 (Unconfined Transient Flow Model), in the plotting section:

From: http://modflowpy.github.io/flopydoc/tutorial2.html

A 'headobj' is first defined:

headobj = bf.HeadFile(modelname+'.hds')

and the heads are extracted as follows:

head = headobj.get_data(totim=time)

Runs on Debian

answered May 11 '15 at 01:53

klobrille

21
2

1

Thanks klobrille , I have figured it out . You have to use the following options while compiling the Modflow. In makefile, use F90= gfortran . In openspec.inc change DATA ACCESS/'SEQUENTIAL'/ toDATA ACCESS/'STREAM'/ change DATA FORM/'BINARY'/ to DATA FORM/'UNFORMATTED'/ – Kirubaharan J May 18 '15 at 14:19

How to convert a unformatted fortran file(modflow output ) to numpy array

2 Answers2