ml.aquifer_summary() doesn't work for "Model"?

[Vera-geo]

Hi! I am a beginner and just trying to figure it all out.
Why does the code below work for a "ModelMaq" model but not for a "Model"?
ml.aquifer_summary()

Looking at the references in the help, it should be the same line of code. I build 2 different scripts one for Model and one for ModelMaq, both are called 'ml'.
I hope someone can explain.

[MattBrst]

Hi Vera, if you could post the code snippets, it is often much easier to help you.

[Vera-geo]

Sure!

import matplotlib.pyplot as plt
import numpy as np
import timflow.steady as tfs

ml = tfs.Model(
z=[-1.2, -4.2, -7.2, -15, -25, -40], #maaiveld, depth of layers
kaq=[4, 3, 15, 50], #k of aquifers
c=[440, 6.2, 0.66, 0.3], # c between aquifers
npor=[0.3, 0.3, 0.3, 0.3, 0.3],
ltype=["l", "a", "a", "a", "a"],
)

ax = ml.plots.xsection(params=True, units={"k": "m/d", "c":"d"})
ax.set_frame_on(False)

print(ml.aq.c)

ml.aquifer_summary()

While it does work for:

import timflow.steady as tfs

ml = tfs.ModelMaq(z=[0, -5, -10, -20, -25, -35], kaq=[10, 30, 20], c=[2000, 5000])

ml.aquifer_summary()

[dbrakenhoff]

Hi @Vera-geo, the two main methods of building a Model in timflow are ModelMaq (a stack of aquifers and leaky layers) and Model3D (a stack of aquifers, where resistance to vertical flow is determined by an anisotropy factor). The Model class is a generic entry point that you can use to build more complex layer models (switching between these two options, I believe), but it means you have to provide all the details yourself. I rarely use it, so my guess is you can probably focus on using one of ModelMaq or Model3D?

ml.aquifer_summary() works for both of these Model classes.

You are of course right to point out that it is somewhat confusing that aquifer_summary() doesn't work for the Model class. The reason I didn't implement that is because of the reason above, and its easy to write simple code to to create an aquifer summary for the Maq and 3D class, but slightly more complicated for the more generic Model class. But I will make an issue and add it to the todo list, and for now implement a slightly more informative error that Model is not supported yet.

[dbrakenhoff]

And thanks for sharing your code. Your model can be built with Model3D with topboundary="semi".

ml3d = tfs.Model3D(
    z=[-4.2, -7.2, -15, -25, -40],  # maaiveld, depth of layers
    kaq=[4, 3, 15, 50],  # k of aquifers
    topboundary="semi",
    topres=440,
    topthick=3,
)
ml3d.plots.xsection()

You can play around with kzoverkh to set the anisotropy factor to get the resistances you want. You can see what they are through ml3d.aq.c.

[Vera-geo]

@dbrakenhoff thanks!! So I wasn't going mad. Your reply is so fast, thank you.
@MattBrst: Thanks for trying to help, so fast.

I can focus on using ModelMaq or Model3D for sure, but since I saw it I tried it. Could you explain what you mean by 'switching' between 2 options?

Will you make that adjustion also in the help? For me this text in the help made me believe it should be possible?

help(tfs.model.Model.aquifer_summary)
Help on function aquifer_summary in module timflow.steady.model:

aquifer_summary(self)
Return DataFrame with summary of aquifer(s) parameters in model.

Returns
-------
pandas.DataFrame
    dataframe with summary of aquifer(s) parameters

[dbrakenhoff]

Will you make that adjustion also in the help? For me this text in the help made me believe it should be possible

Yes for now we will add an informative error and improve the docstring. The issue of writing a generic aquifer_summary that also supports the Model class will be put on the issue list. I'm afraid this won't have our priority at the moment (too many other interesting things to do 😇), so we'll label it with "help wanted" and if you feel like having a go, feel free to do so :). We'd be happy to help and take a look at your code if you run into any challenges of course.

[Vera-geo]

@dbrakenhoff Funnily enough I can't seem to make it work for the 3D either, ModelMaq does still work.

ml3d = tfs.Model3D(
z=[-4.2, -7.2, -15, -25, -40], # maaiveld, depth of layers
kaq=[4, 3, 15, 50], # k of aquifers
topboundary="semi",
topres=440,
topthick=3,
)

ml3d.plots.xsection() --> works fine

ml3d.aquifer_summary() --> doesn't work

---

ValueError Traceback (most recent call last)
Cell In[12], line 1
----> 1 ml3d.aquifer_summary()

File ~\AppData\Local\anaconda3\Lib\site-packages\timflow\steady\model.py:592, in Model.aquifer_summary(self)
590 aqs = {}
591 if not isinstance(self.aq, SimpleAquifer):
--> 592 aqs["background"] = self.aq.summary()
593 for i, iaq in enumerate(self.aq.inhomlist):
594 aqs[f"inhom{i}"] = iaq.summary()

File ~\AppData\Local\anaconda3\Lib\site-packages\timflow\steady\aquifer.py:183, in AquiferData.summary(self)
181 else:
182 summary.iloc[~maskaq, 2] = self.Hll
--> 183 summary.iloc[~maskaq, 4] = self.c
184 summary.loc[:, "layer"] = self.layernumber
185 return summary

File ~\AppData\Local\anaconda3\Lib\site-packages\pandas\core\indexing.py:912, in _LocationIndexer.setitem(self, key, value)
909 self._has_valid_setitem_indexer(key)
911 iloc = self if self.name == "iloc" else self.obj.iloc
--> 912 iloc._setitem_with_indexer(indexer, value, self.name)

File ~\AppData\Local\anaconda3\Lib\site-packages\pandas\core\indexing.py:1943, in _iLocIndexer._setitem_with_indexer(self, indexer, value, name)
1940 # align and set the values
1941 if take_split_path:
1942 # We have to operate column-wise
-> 1943 self._setitem_with_indexer_split_path(indexer, value, name)
1944 else:
1945 self._setitem_single_block(indexer, value, name)

File ~\AppData\Local\anaconda3\Lib\site-packages\pandas\core\indexing.py:1999, in _iLocIndexer._setitem_with_indexer_split_path(self, indexer, value, name)
1994 if len(value) == 1 and not is_integer(info_axis):
1995 # This is a case like df.iloc[:3, [1]] = [0]
1996 # where we treat as df.iloc[:3, 1] = 0
1997 return self._setitem_with_indexer((pi, info_axis[0]), value[0])
-> 1999 raise ValueError(
2000 "Must have equal len keys and value "
2001 "when setting with an iterable"
2002 )
2004 elif lplane_indexer == 0 and len(value) == len(self.obj.index):
2005 # We get here in one case via .loc with a all-False mask
2006 pass

ValueError: Must have equal len keys and value when setting with an iterable

[dbrakenhoff]

Intriguiing, it's working for me:

The dev branch probably includes a fix. If you're comfortable installing the develop version of timflow:

pip install git+https://github.com/timflow-org/timflow.git@dev

Otherwise download the repository when you've selected the dev branch from github as a zip, unpack it and use pip install . from inside the unpacked directory.

And otherwise, we are preparing a new release early next week, so after that a simple pip install -U timflow will do the trick.

[Vera-geo]

Yes installed the dev, that did the trick!