[FEATURE] - Add reader for LAK STAGE FILEOUT binary file

[FTW12345]

Title
Feature request: imod.mf6 utility to read LAK stagefile.bin into an xarray object
Problem
To visualise seasonal lake dynamics, we need time series of lake stage per stress period.
When lake.lak, options contains 'STAGE FILEOUT C:/.../stagefile.bin', then MF6 writes the binary file stagefile.bin. Alas, iMOD Python lacks a function to read the written binary file.
Solution
Proposed API
• imod.mf6.read_lak_stage(path, nlakes=None) -> xr.DataArray
o dims: ("time", "lake")
o coords: time from tdis, lake as 1..NLAKES (or names if available)
o dtype: float64
Reference structure (tested)
Each record contains:
• int32 kstp, int32 kper, float64 pertim, float64 totim
• 16-char text label ("STAGE")
• int32 n1, int32 n2, int32 n3 (counters)
• then float64 stage values (NLAKES entries)
Minimal reader we used (works for 1 lake)

import struct
from pathlib import Path

def read_lak_stagebin(bin_path: Path):
    recs = []
    with open(bin_path, "rb") as f:
        while True:
            hdr = f.read(24)
            if len(hdr) < 24:
                break
            kstp, kper, pertim, totim = struct.unpack("<iid d".replace(" ",""), hdr)
            label = f.read(16).decode("ascii", "ignore").strip()
            if label.upper() != "STAGE":
                raise ValueError(f"Unexpected label: {label!r}")
            f.read(12)  # three int32
            stage = struct.unpack("<d", f.read(8))[0]
            recs.append((kper, kstp, pertim, totim, stage))
    return recs

[JoerivanEngelen]

Hi @FTW12345, we are indeed lacking this functionality. And thanks for writing a concise issue. The code example provides a good starting point. The Lake package hasn't seen much use yet by us, so that's probably why you're running in these problems.

I hope I can carve out some time to do this on the side, the Lake package isn't the main focus for us at the moment, as we currently do not have any projects applying it. You can see what we have broadly planned for iMOD on our roadmap. If you are willing to create a PR for this yourself with some guidance, let me know in this issue and I'll point you to the locations where logic should be put.

NB: I modified your post slightly by putting ```python around your example code, so that it is properly formatted in the issue.

[FTW12345]

Hello Joeri,
Thanks for your reply.
I'm happy to help and create a PR.
However, I have no experience with Pull Requests.
I am thus unfamiliar with what my contribution can be beyond what I already send in with the issue (meaning I am not a Python programmer).
Martin

[FTW12345]

Hello Joeri,
Reading on our roadmap I see your Wishlist ... "iMOD Python | Export data to VTK".
FYI: I created c. 300 script lines to transform the MODFLOW 6 output files .hds, .cbc, .ucn, lakestage.bin and lake_conc.bin into .pvd and belonging *.vts files, ready to open in Paraview.
Would they be of any help?
I presume that using them for a new imod python function, requires multiple specific imod python checks to add (with which I am unfamiliar).
Martin

[JoerivanEngelen]

Hi Martin,

I just spent a bit more time investigating this, it turns out that we DO have a function for this, added by my colleague @HendrikKok!
It is a bit obscure, but it is named imod.mf6.open_dvs. Could you test if this works for you?

We could do with some clearer naming for the function names so that it is more findable for users. For example by creating functions with names like open_lak_stage, open_uzf_water_content, open_sfr_stage, which just call open_dvs behind-the-scenes. What do you think?

Considering your Python coding experience and the PR, I started my career as a PhD student at the UU in geohydrology as well (inferring you too from your profile), which for me was a good starting point to gain experience in coding. Contributing to packages is a fast way to learn new coding concepts and tricks, it however comes with the frustration that package maintainers might ask you to rewrite certain parts of code or reject some of your ideas (but still less frustrating as getting a paper rejected, in my experience, because time scales are smaller).
I get though that time flies in academia, there are many people (university, professors, teachers, students) who want things from you already, and your main focus is finishing your project successfully.

Considering the VTK exported code, would be great if you could share it in this issue that I just created: #1759 !

[FTW12345]

Hello Joeri,
I looked at the help that is available for open_dvs (https://deltares.github.io/imod-python/api/generated/mf6/imod.mf6.open_dvs.html#imod.mf6.open_dvs): open_dvs uses the 3D model grid (dis.dis.grb).
The LAK results (stagefile.bin; stage of lake) and the LKT results (conc_lake.bin; concentration in the lake) are 1D.
Therefore, the function open_dvs is alas not applicable as reader for LAK and LKT results.
This means that we can start a Pull Request for the imod python functions:

• imod.mf6.open_lak_stage;
• imod.mf6.open_lkt_conc.

Martin

[JoerivanEngelen]

You are right, I've investigated a bit further with @HendrikKok and here are our findings:

TL;DR: The proposed changes would entail quite a refactoring, so I think doing this as your first PR to iMOD Python is probably too ambitious.

Ugly hack with current code base

imod.mf6.open_dvs has an indices argument, which is flat single array used to map dvs data to the MODFLOW grid. This works fine for the UZF, as a single UZF cell can only be coupled to one MODFLOW cell (or at least that is the most common usecase). For lakes however, this is not the case. One lake can be assigned to multiple modflow cells, and multiple lakes can be assigned to one modflow cell. As an ugly workaround to get what you want, a 2d array of size (ntime, nlake), we can abuse the indices argument to temporarily assign lake stages to a modflow cellid that equals the lake index, and then selecting this again. This is less performant than your solution, but might be good enough for your usecase. I tested this out on this example here and it seemed to work:

# Work around the mapping of lakes to modflow grid cells, by assigning lake
# stages to modflow cellids that are equal to lake numbers, then selecting those
# lake numbers again.
lake_numbers = np.array([0])
stages_grid_scratch = imod.mf6.open_dvs(stagefile, grbfile, lake_numbers)
cellid_dims = ["layer", "y", "x"]
stages_lake = stages_grid_scratch.stack(cellid=cellid_dims).isel(cellid=lake_numbers)

Flaws of the current approach

The idea behind how open_dvs now works is that the variables are mapped from a 1D array to the data in the form users have provided to construct the model with the indices argument, in the case of the UZF package this is a 3D grid. However, I've noticed the following flaws of the current approach:

1. It's now easy to make a mistake in determining these indices, and there is no input validation. Providing indices with the wrong dimensions result in garbage data with the wrong time dimension. So it would be better if users need to provide a grid with UZF numbers directly, which is used to assign uzf data to the grid. This can be generated with UnsaturatedZoneFlow._create_uzf_numbers(), which we could expose as a public classmethod. Furthermore, for extra safety, we could let users provide times themselves, so that ntime can be provided from that, instead of start_simulation_data.
2. Furthermore, how lakes should be mapped with the present approach is less evident: You can couple multiple modflow cells to 1 lake and multiple lakes to 1 modflow cell with the current API we have in iMOD Python for LakeData and Lake. I think easiest is to map lake stages to the lake_dim of in Lake(...).dataset, so as a 1D array. Alternatively, we could map these to the model grid, but since multiple lakes can be connected to a single modflow cell, this would require either accepting a grid of shape (lake_dim, layer, y, x) and then returning (time, lake_dim, layer, y, x), which can become quite huge with multiple lakes. Alternatively, we could aggregate this data, as we do for the cbc data with a summation, but how lake stages should be aggregated I find less evident.
3. In addition, for the SFR package, which we plan to support in the future at some point, we plan to use a UGRID1D network to specify the streams. Ideally, we would like to map SFR data back to this UGRID network. This has an entirely different shape than the model grid.

Proposed approach

We think the way forward is to change open_dvs, so that that only returns 2d arrays of shape (ntime, n_dvs), based on an ntime. This can then be mapped to a shape, based on the package type:

• open_uzf_water_content will map uzf data to the model grid, so have dims (time, layer, y, x)
• open_lake_stages will not map lake stages to a grid, and will have dims (time, lake_dim)
• open_sfr_stage would map to a sfr stages to the 1D ugrid network, and will have (time, sfr_edge)

In general, we set the following requirements to ourselves when changing/adding new features to iMOD Python:

• Avoid additional code duplication as much as possible
• Present functionality should be preserved in some way as much as possible, so in this case the functionality to map UZF nodes to the model grid should still be preserved.

Broadly speaking, we would therefore have to change the following:

• Either generalize reshaping logic from read_dvs_timestep functions or move this logic to a later stage, somewhere in the package specifc functions. If the latter is chosen, we should test if this doesn't hit performance too badly for large grids, given that the current approach nicely reshapes for each chunk very early in the reading routine.
• Think of clever way to determine ntime based on user input
• Implement separate open_lake_stages and open_uzf_water_content functions.

This goes beyond just adding two extra functions, unfortunately, so I reckon this might be a bit too much for a first contribution.

[Huite]

You could consider adding the output reading method as a function to the relevant classes.

I think imod.mf6.Lake.open_stages() reads naturally. Best option seems to make it a static method which takes I guess takes a DIS grid (or something) to extract the topology?
Same for UZF, etc.