The latest version of IWFM uses hdf output file format. This makes the reading output files a bit more difficult since the outputs are not in ASCII format but way more efficient. Luckily Matlab, R and many other languages can handle hdf formats quite easily.
Groundwater Zone Budget file
The first step to read any hdf file in matlab is to read the information that is contained in the file
GB_bud_info = h5info(fullfile(c2vsim_path, 'Results','C2VSimFG_GW_ZBudget.hdf'));
The GB_bud_info is a structure with many fields that contain plenty of useful information. After examining the structure in my version the GB_bud_info contained the field Groups, which has the following structure
If we expand the dataset for the Layer1 we see the following structure
From this table we can find the row where each property is located . To read the data for a specific property we use the following. The Storage for example is the row 22. The following snippet loops through the layers and adds the entire storage from all layers
CVstorage = zeros(504, 32537); for ii = 1:4 CVstorage = CVstorage + h5read(GB_bud_info.Filename,... [GB_bud_info.Groups(1+ii).Name GB_bud_info.Name GB_bud_info.Groups(1+ii).Datasets(22).Name])'; end
The units of the hdf file are in feet. Since this is storage the units are cubic feet. To generate the famous storage depletion plot of Central valley for example after we have calculated the above we can do the following
% Create the time variable for the plot start_date = datetime(1973,10,1); tm = start_date + calmonths(0:503); % Calculate the storage from all elements CVcumulativeStorage = sum(CVstorage,2); % convert from feet to acre feet CVcumulativeStorage = CVcumulativeStorage * 2.29569e-5 % Convert to million acrefeet CVcumulativeStorage = CVcumulativeStorage /1000000 % and do the plot plot(tm,CVcumulativeStorage - CVcumulativeStorage(1),'linewidth',2) xlabel('Time') ylabel('Storage depletion [MAF]') grid on
In the case of storage the data have dimensions as we see in the ChunkSize field 32537 and we can somewhat safely assume that the 1st row correspond to 1st element etc. However the ChunkSize for the Streams (rows 23,24) is 9885.
StreamBud=h5read(GB_bud_info.Filename,... [GB_bud_info.Groups(2).Name GB_bud_info.Name GB_bud_info.Groups(2).Datasets(23).Name])';
The size of the StreamBud variable is 504×9885. Therefore from this matrix we cannot tell which element corresponds to the values. This information exists under Group(1).Datasets
The 6th row contains the correspondence between data and elements. To read it use
ids = h5read(GB_bud_info.Filename,[GB_bud_info.Groups(1).Name GB_bud_info.Name GB_bud_info.Groups(1).Datasets(6).Name])';
which returns a matrix 24×32537 (note that the above code transpose the output). Obviously the columns are related to the elements but how about the rows? This information exists in the 5th row of the above variable which we can read it as
Colnames = h5read(GB_bud_info.Filename,[GB_bud_info.Groups(1).Name GB_bud_info.Name GB_bud_info.Groups(1).Datasets(5).Name]);
The Colnames variable consist of 24 rows! The 3rd row is the Streams_Inflow (+). The elements that are associated with the stream budget are the nonzero ones in the 3rd row of the ids matrix which can be extracted easily
stream_elem_ids = find(ids(3,:) ~=0);
To be continued....