In many cases its very useful to be able to read and process in batch mode the output file of Modflow (Optimization, inverse modeling, Monte Carlo simulation etc...)
The USGS has provided a utility program to convert the binary output to ascii, which is easier to read, yet as you will see below it is rather straightforward to read the binary files without using an intermediate step.
In general, reading binary files requires precise knowledge of how the data are written on the file. Finding this information was actually more difficult than writing the program itself. Go to http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/ and then and then under the FAQ look for "How can I read data from a binary file generated by MODFLOW?". At the bottom of the description click the plus next to Array Data and there they are!!.
Let's suppose we want to read the "mymodel.hds" binary file that contains the head field of a multi-layer aquifer with Nlay layers.
To read the file we will use few Matlab/Octave commands, fopen/fclose, fread.
>>fid=fopen("mymodel.hds",'r'); % open the file for reading
Next we will start reading the data but as the description suggest we should know whether the data are written in single-precision or double-precision format. As we would never know beforehand the correct format, we just make an assumption and check this later.
First we need to read an integer of 4 bytes, which is the the time step number.
>>KPER = fread(fid, 1, 'uint');% read the stress period number
>>PERTIM = fread(fid, 1, 'float'); % the time in the current stress period, which is a real number and that's why we use 'float' format.
>> TOTIM = fread(fid, 1, 'float'); %the total elapsed time, which again is a real number.
Next we have to read a word, which is actually the moment of truth. If the assumptions about the single/double precision format or the number of bytes (4 or 8) for the real numbers are correct, we should read something that makes sense. When we read a head field, for example, we should simply read : ' HEAD'.
According to the description we have to read 16 ANSI characters
[32 32 32 32 32 72 69 65 68], which are the ANSI characters of the word we just read.
To convert it to a readable format use:
If DESC matches one of the values of the table, on the FAQ then our assumptions are correct.
Next we read 3 integer numbers
>>NCOL = fread(fid, 1, 'uint'); % number of columns
>>NROW = fread(fid, 1, 'uint'); % number of rows
>>ILAY = fread(fid, 1, 'uint'); % layer number
Finally we have to read the head field itself using a loop
In case of multi layer we need to put all the above commands into a loop.
Last we need to close the file (which is something I always forget)
>>fclose(fid);
(I have used many times the above implementation without any problems. Therefore I cannot give recommendations if the file uses a different format).
The USGS has provided a utility program to convert the binary output to ascii, which is easier to read, yet as you will see below it is rather straightforward to read the binary files without using an intermediate step.
In general, reading binary files requires precise knowledge of how the data are written on the file. Finding this information was actually more difficult than writing the program itself. Go to http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/ and then and then under the FAQ look for "How can I read data from a binary file generated by MODFLOW?". At the bottom of the description click the plus next to Array Data and there they are!!.
Let's suppose we want to read the "mymodel.hds" binary file that contains the head field of a multi-layer aquifer with Nlay layers.
To read the file we will use few Matlab/Octave commands, fopen/fclose, fread.
>>fid=fopen("mymodel.hds",'r'); % open the file for reading
Next we will start reading the data but as the description suggest we should know whether the data are written in single-precision or double-precision format. As we would never know beforehand the correct format, we just make an assumption and check this later.
First we need to read an integer of 4 bytes, which is the the time step number.
>>KSPT = fread(fid, 1, 'uint');
similarly:>>KPER = fread(fid, 1, 'uint');% read the stress period number
>>PERTIM = fread(fid, 1, 'float'); % the time in the current stress period, which is a real number and that's why we use 'float' format.
>> TOTIM = fread(fid, 1, 'float'); %the total elapsed time, which again is a real number.
Next we have to read a word, which is actually the moment of truth. If the assumptions about the single/double precision format or the number of bytes (4 or 8) for the real numbers are correct, we should read something that makes sense. When we read a head field, for example, we should simply read : ' HEAD'.
According to the description we have to read 16 ANSI characters
>> DESC = fread(fid, 16, 'char');
The DESC now is an array that looks like this:[32 32 32 32 32 72 69 65 68], which are the ANSI characters of the word we just read.
To convert it to a readable format use:
>>DESC=char(DESC');
DESC is now the word ' HEAD' in ascii characters.If DESC matches one of the values of the table, on the FAQ then our assumptions are correct.
Next we read 3 integer numbers
>>NCOL = fread(fid, 1, 'uint'); % number of columns
>>NROW = fread(fid, 1, 'uint'); % number of rows
>>ILAY = fread(fid, 1, 'uint'); % layer number
Finally we have to read the head field itself using a loop
>> for ir = 1 : NROW
for ic = 1 : NCOL
H(ir,ic,ILAY)=fread(fid, 1, float);
end
end
In case of multi layer we need to put all the above commands into a loop.
Last we need to close the file (which is something I always forget)
>>fclose(fid);
(I have used many times the above implementation without any problems. Therefore I cannot give recommendations if the file uses a different format).
10 comments:
Thanks. This information is very helpful.
How would you read .ddn files?
Do you have a matlab script for that too? can you upload it?
Thanks
I have used a very similar code to read the cell by cell flow file, which is more complex than the head and drawdown.
If you can wait few days I'll write the script.
I can provide a "R" function reading the heads of each layer from the binary file. Interested?
Hi George, I am looking for a Matlab function to read a binary cell by cell flow file. Do you have one you could share? Thanks!
I have a code but its not general purpose. Its just a script that reads few types of properties of the cell-by-cell file. You may need to modify it a bit.
Send me an email at gkourakos at ucdavis dot edu
For those who do not visit this post accidentally you might be also interested in this:
http://www.subsurface.gr/index.php/software/modflow-utilities
Is there a way to do this for those of us without MATLAB?
If you have linux or can install it as second OS, the script should work with Octave(Octave is a free alternative to matlab). It is possible to install Octave under windows but its not going to be easy.
If you dont use matlab, then choosing how to read the Modflow output depends on how you intent to use the results of modflow. In my case I'm always post-processing the Modflow results with Matlab. If you would like to use ARCGIS for example, then a good choise is to write a python script to read the files so that you can take advantage of the arcpy extension
In general I think you can convert the code to any modern language. Many years ago I had done this using Visual Basic.
C/C++, python, R, all these languages provide the methods to read all type of numbers
Hello,
I know this is old but if the person who volunteered the R script sees this, I would greatly appreciate it. I can be contacted at bbuzz31@yahoo.com
Thanks,
Brett
Post a Comment