Xarray error when using nlmod.dims.get_modellayers_screens

[tomvansteijn]

Describe the bug

Xarray error when using nlmod.dims.get_modellayers_screens.

Expected behavior

Get layer numbers for well location and screen depth (z-coordinate).

Code to reproduce

ds = nlmod.get_ds([0., 1e3, 0., 1e3], model_ws="test")
ds = nlmod.grid.refine(ds)
nlmod.dims.get_modellayers_screens(ds, screen_top=0., screen_bottom=-10., icell2d=10)

Package and Python version

Python version     : 3.13.7
NumPy version      : 2.3.3
Xarray version     : 2025.12.0
Matplotlib version : 3.10.6
Flopy version      : 3.10.0

nlmod version      : 0.11.2

Additional context

Stack trace

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File d:\temp\conda-envs\bw_basis\Lib\site-packages\xarray\core\variable.py:[153](file:///D:/temp/conda-envs/bw_basis/Lib/site-packages/xarray/core/variable.py:153), in as_variable(obj, name, auto_convert)
    152 try:
--> 153     obj = Variable(dims_, data_, *attrs)
    154 except (TypeError, ValueError) as error:

File d:\temp\conda-envs\bw_basis\Lib\site-packages\xarray\core\variable.py:399, in Variable.__init__(self, dims, data, attrs, encoding, fastpath)
    379 """
    380 Parameters
    381 ----------
   (...)    397     unrecognized encoding items.
    398 """
--> [399](file:///D:/temp/conda-envs/bw_basis/Lib/site-packages/xarray/core/variable.py:399) super().__init__(
    400     dims=dims, data=as_compatible_data(data, fastpath=fastpath), attrs=attrs
    401 )
    403 self._encoding: dict[Any, Any] | None = None

File d:\temp\conda-envs\bw_basis\Lib\site-packages\xarray\namedarray\core.py:261, in NamedArray.__init__(self, dims, data, attrs)
    260 self._data = data
--> [261](file:///D:/temp/conda-envs/bw_basis/Lib/site-packages/xarray/namedarray/core.py:261) self._dims = self._parse_dimensions(dims)
    262 self._attrs = dict(attrs) if attrs else None

File d:\temp\conda-envs\bw_basis\Lib\site-packages\xarray\namedarray\core.py:505, in NamedArray._parse_dimensions(self, dims)
    504 if len(dims) != self.ndim:
--> [505](file:///D:/temp/conda-envs/bw_basis/Lib/site-packages/xarray/namedarray/core.py:505)     raise ValueError(
    506         f"dimensions {dims} must have the same length as the "
    507         f"number of data dimensions, ndim={self.ndim}"
    508     )
    509 if len(set(dims)) < len(dims):

ValueError: dimensions ('icell2d',) must have the same length as the number of data dimensions, ndim=0

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
Cell In[129], [line 3](vscode-notebook-cell:?execution_count=129&line=3)
      1 ds = nlmod.get_ds([0., 1e3, 0., 1e3], model_ws="test")
      2 ds = nlmod.grid.refine(ds)
----> [3](vscode-notebook-cell:?execution_count=129&line=3) nlmod.dims.get_modellayers_screens(ds, screen_top=0., screen_bottom=-10., icell2d=10)

File d:\temp\conda-envs\bw_basis\Lib\site-packages\nlmod\dims\layers.py:2158, in get_modellayers_screens(ds, screen_top, screen_bottom, xy, icell2d)
   2156 # make dataset of observations
   2157 ds_obs = ds.sel(icell2d=icell2d)
-> [2158](file:///D:/temp/conda-envs/bw_basis/Lib/site-packages/nlmod/dims/layers.py:2158) ds_obs["screen_top"] = (("icell2d"), screen_top)
   2159 ds_obs["screen_bot"] = (("icell2d"), screen_bottom)
   2160 dimname = "icell2d"

File d:\temp\conda-envs\bw_basis\Lib\site-packages\xarray\core\dataset.py:1405, in Dataset.__setitem__(self, key, value)
   1400     if isinstance(value, Dataset):
   1401         raise TypeError(
   1402             "Cannot assign a Dataset to a single key - only a DataArray or Variable "
   1403             "object can be stored under a single key."
   1404         )
-> [1405](file:///D:/temp/conda-envs/bw_basis/Lib/site-packages/xarray/core/dataset.py:1405)     self.update({key: value})
   1407 elif utils.iterable_of_hashable(key):
   1408     keylist = list(key)

File d:\temp\conda-envs\bw_basis\Lib\site-packages\xarray\core\dataset.py:5637, in Dataset.update(self, other)
   5609 def update(self, other: CoercibleMapping) -> None:
   5610     """Update this dataset's variables with those from another dataset.
   5611 
   5612     Just like :py:meth:`dict.update` this is a in-place operation.
   (...)   5635     Dataset.merge
   5636     """
-> [5637](file:///D:/temp/conda-envs/bw_basis/Lib/site-packages/xarray/core/dataset.py:5637)     merge_result = dataset_update_method(self, other)
   5638     self._replace(inplace=True, **merge_result._asdict())

File d:\temp\conda-envs\bw_basis\Lib\site-packages\xarray\structure\merge.py:1220, in dataset_update_method(dataset, other)
   1217                 value = value._replace(variable=variable)
   1218             other[key] = value
-> [1220](file:///D:/temp/conda-envs/bw_basis/Lib/site-packages/xarray/structure/merge.py:1220) return merge_core(
   1221     [dataset, other],
   1222     compat="broadcast_equals",
   1223     join="outer",
   1224     priority_arg=1,
   1225     indexes=dataset.xindexes,
   1226     combine_attrs="override",
   1227 )

File d:\temp\conda-envs\bw_basis\Lib\site-packages\xarray\structure\merge.py:767, in merge_core(objects, compat, join, combine_attrs, priority_arg, explicit_coords, indexes, fill_value, skip_align_args)
    764 for pos, obj in skip_align_objs:
    765     aligned.insert(pos, obj)
--> [767](file:///D:/temp/conda-envs/bw_basis/Lib/site-packages/xarray/structure/merge.py:767) collected = collect_variables_and_indexes(aligned, indexes=indexes)
    768 prioritized = _get_priority_vars_and_indexes(aligned, priority_arg, compat=compat)
    769 variables, out_indexes = merge_collected(
    770     collected,
    771     prioritized,
    772     compat=compat,
    773     combine_attrs=combine_attrs,
    774 )

File d:\temp\conda-envs\bw_basis\Lib\site-packages\xarray\structure\merge.py:409, in collect_variables_and_indexes(list_of_mappings, indexes)
    406     indexes_.pop(name, None)
    407     append_all(coords_, indexes_)
--> [409](file:///D:/temp/conda-envs/bw_basis/Lib/site-packages/xarray/structure/merge.py:409) variable = as_variable(variable, name=name, auto_convert=False)
    410 if name in indexes:
    411     append(name, variable, indexes[name])

File d:\temp\conda-envs\bw_basis\Lib\site-packages\xarray\core\variable.py:155, in as_variable(obj, name, auto_convert)
    153         obj = Variable(dims_, data_, *attrs)
    154     except (TypeError, ValueError) as error:
--> [155](file:///D:/temp/conda-envs/bw_basis/Lib/site-packages/xarray/core/variable.py:155)         raise error.__class__(
    156             f"Variable {name!r}: Could not convert tuple of form "
    157             f"(dims, data[, attrs, encoding]): {obj} to Variable."
    158         ) from error
    159 elif utils.is_scalar(obj):
    160     obj = Variable([], obj)

ValueError: Variable 'screen_top': Could not convert tuple of form (dims, data[, attrs, encoding]): ('icell2d', 0.0) to Variable.

[OnnoEbbens]

Hi Tom, get_modellayers_screens expects a collection of screens so you won't get the error if you use:

nlmod.dims.get_modellayers_screens(ds, screen_top=[0.], screen_bottom=[-10.], icell2d=[10])

We could modify get_modellayers_screens function to deal with a single well screen. I think it won't be too hard.

@tomvansteijn: Maybe if you are up for it you can draft a PR?

[tomvansteijn]

Ah I see. I didn't realize that the function requires lists as input. I can do this. Is it a good idea to extend the function to optionally return all layers that intersect with the screen?

[tomvansteijn]

So this is not a bug, I should change the label but I don't know how..

[OnnoEbbens]

I am not sure if you are allowed to remove a label, so I removed it.

I don't think it is easy to adapt get_modellayers_screens to include all modellayers because it will add an extra dimension to the output array. However you can get the top and bottom modellayer of each well screen using get_modellayers_indexer with full_output=True, see the example below:

ds = nlmod.get_ds([0., 1e3, 0., 1e3], botm=[-1,-7,-11,-15], model_ws="test")
ds = nlmod.grid.refine(ds)
df = pd.DataFrame({'x':[100], 'y':[100], 'screen_top':[0.], 'screen_bottom':[-10.]})
ds_obs = nlmod.dims.get_modellayers_indexer(ds, df, full_output=True)
ds_obs[['modellayer_top', 'modellayer_bot']]

I get that this is a bit confusing. The functions are part of a workflow for comparing model results with observations and not very easy to call separately. This notebook has an example on how to use the full workflow.