From fbd391e502d3615d9871fe723b39577e9487abc3 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Thu, 1 Oct 2026 17:08:09 +0200
Subject: [PATCH 01/27] Start working on opening idfs by time

---
 imod/formats/array_io/reading.py |  12 +-
 imod/formats/idf.py              | 189 +++++++++++++++++++++++++------
 2 files changed, 163 insertions(+), 38 deletions(-)

diff --git a/imod/formats/array_io/reading.py b/imod/formats/array_io/reading.py
index 412513214..ddde13c4e 100644
--- a/imod/formats/array_io/reading.py
+++ b/imod/formats/array_io/reading.py
@@ -280,8 +280,8 @@ def _dask(path, attrs=None, pattern=None, _read=None, header=None):
     return x, attrs
 
 
-def _load(paths, use_cftime, _read, headers):
-    """Combine a list of paths to IDFs to a single xarray.DataArray"""
+def _load_as_components(paths, use_cftime, _read, headers):
+    """Combine a list of paths to IDFs to a single xarray.DataArray and return as components"""
     # this function also works for single IDFs
     names = [h["name"] for h in headers]
     _all_equal(names, "names")
@@ -320,8 +320,14 @@ def _load(paths, use_cftime, _read, headers):
             nested_dict.set_nested(groupby, groupbykeys, da)
         dask_arrays = nested_dict.sorted_nested_dict(groupby)
         dask_array = _ndconcat(dask_arrays, ndim)
+    return dask_array, coords, dims, names[0]
+
 
-    out = xr.DataArray(dask_array, coords, dims, name=names[0])
+def _load(paths, use_cftime, _read, headers):
+    dask_array, coords, dims, name = _load_as_components(
+        paths, use_cftime, _read, headers
+    )
+    out = xr.DataArray(dask_array, coords, dims, name=name)
 
     first_attrs = headers[0]
 
diff --git a/imod/formats/idf.py b/imod/formats/idf.py
index 68eb1e1c8..c0443e8b6 100644
--- a/imod/formats/idf.py
+++ b/imod/formats/idf.py
@@ -12,7 +12,7 @@
 from collections.abc import Iterable
 from pathlib import Path
 from re import Pattern
-from typing import Any, DefaultDict
+from typing import Any, Callable, DefaultDict
 
 import dask
 import dask.array
@@ -20,7 +20,7 @@
 import xarray as xr
 
 import imod
-from imod.formats import array_io
+from imod.formats import _load_as_components, array_io
 from imod.typing.structured import merge_partitions_as_da_components
 
 # Make sure we can still use the built-in function...
@@ -291,6 +291,150 @@ def check_subdomain_consistency(
         )
 
 
+def _open_idf_data(
+    paths: list[str],
+    use_cftime: bool,
+    pattern: str | Pattern,
+):
+    """Open IDF files and return as a DataArray."""
+    headers = [header(p, pattern) for p in paths]
+    data, _, _, _ = _load_as_components(
+        paths, use_cftime=use_cftime, _read=array_io.read, headers=headers
+    )
+    return data
+
+
+def _open_idf_chunked_by_time(
+    open_func: Callable,
+    grouped_by_time: dict[Any, Any],
+    template: xr.DataArray,
+    has_time: bool,
+    use_cftime: bool,
+    pattern: str | Pattern,
+):
+    """
+    Open IDF files chunked by time using a template DataArray.
+
+    Parameters
+    ----------
+    open_func : Callable
+        Function to open individual IDF files. Needs to have group, use_cftime,
+        pattern as arguments. Needs to return dask array as data.
+    grouped_by_time : dict[Any, Any]
+        Dictionary grouping file paths by time.
+    template : xr.DataArray
+        Template DataArray to define the shape, dims, and coords.
+    has_time : bool
+        Whether the data has a time dimension.
+    use_cftime : bool
+        Whether to use cftime for time coordinates.
+    pattern : str | Pattern
+        Pattern to match file paths.
+
+    Returns
+    -------
+    xr.DataArray
+        Combined DataArray with data from all subdomains, chunked by time.
+
+
+    """
+    shape = template.shape  # e.g. (1, nlayer, nrow, ncol)
+    dims = template.dims  # e.g. ("time", "layer", "y", "x")
+    dtype = template.dtype
+    if has_time:
+        time_axis = list(dims).index("time")
+    else:
+        time_axis = -1  # steady-state, no time dimension
+
+    # Sort and convert times before calling _merge_subdomains so that
+    # use_cftime is already correct when the template is built.
+    raw_times_sorted = sorted(grouped_by_time.keys())
+    converted_times, use_cftime = imod.util.time._convert_datetimes(
+        raw_times_sorted, use_cftime
+    )
+
+    # Delayed tasks for each timestep, which will be concatenated into a single
+    # dask array. One delayed task per timestep → outer graph depth 3, O(n_time)
+    # tasks
+    dask_arrays = []
+    for time_key in raw_times_sorted:
+        group = grouped_by_time[time_key]
+        timestep_data = dask.delayed(open_func)(group, use_cftime, pattern)
+        dask_arrays.append(
+            dask.array.from_delayed(timestep_data, shape=shape, dtype=dtype)
+        )
+    data = dask.array.concatenate(dask_arrays, axis=time_axis)
+
+    # Build the full time coordinate
+    coords = dict(template.coords)
+    if has_time:
+        if use_cftime:
+            time_coord = xr.CFTimeIndex(converted_times)
+        else:
+            time_coord = np.array(converted_times, dtype="datetime64[ns]")
+        coords["time"] = time_coord
+
+    return xr.DataArray(data, coords, dims, name=template.name, attrs=template.attrs)
+
+
+def open_idf_by_time(path, use_cftime=False, pattern=None):
+    """
+    Open IDF files grouped by time.
+
+    Parameters
+    ----------
+    path : str or Path
+        Global path.
+    use_cftime : bool, optional
+        Whether to use cftime for time coordinates.
+    pattern : str, regex pattern, optional
+        If no pattern is provided, the function will first try:
+        "{name}_c{species}_{time}_l{layer}_p{subdomain}"
+        and if that fails:
+        "{name}_{time}_l{layer}_p{subdomain}"
+        Following the iMOD5/iMOD-WQ filename conventions.
+
+    Returns
+    -------
+    xarray.DataArray
+
+    """
+
+    paths = sorted(glob.glob(str(path)))
+    parsed = [imod.util.path.decompose(path, pattern) for path in paths]
+
+    has_time = "time" in parsed[0]
+
+    # Group by time (datetime.datetime from decompose), then by subdomain.
+    # Each delayed task processes one timestep, keeping the outer graph at O(n_time).
+    grouped_by_time: DefaultDict[Any, list] = defaultdict(list)
+
+    for match, p in zip(parsed, paths):
+        if has_time:
+            time_key = match["time"]
+        else:
+            # Work around for files without time dimension
+            # (imod.util.time._convert_datetimes special-cases this string)
+            time_key = "steady-state"
+        grouped_by_time[time_key].append(p)
+
+    # Sort and convert times before calling idf.open so that
+    # use_cftime is already correct when the template is built.
+    raw_times_sorted = sorted(grouped_by_time.keys())
+
+    # Call idf.open eagerly for the first timestep to obtain a
+    # coordinate template. No data is computed — only the coordinate arrays
+    # (which are numpy) are used; the dask data array is discarded.
+    first_time_key = raw_times_sorted[0]
+    template = open(
+        grouped_by_time[first_time_key], use_cftime=use_cftime, pattern=pattern
+    )
+
+    return _open_idf_chunked_by_time(
+        _open_idf_data, grouped_by_time, template, has_time, use_cftime, pattern
+    )
+
+
 def open_subdomains(
     path: str | Path, use_cftime: bool = False, pattern: str | Pattern = None
 ) -> xr.DataArray:
@@ -352,9 +496,6 @@ def open_subdomains(
     # Sort and convert times before calling _merge_subdomains so that
     # use_cftime is already correct when the template is built.
     raw_times_sorted = sorted(grouped_by_time.keys())
-    converted_times, use_cftime = imod.util.time._convert_datetimes(
-        raw_times_sorted, use_cftime
-    )
 
     # Call _merge_subdomains eagerly for the first timestep to obtain a
     # coordinate template. No data is computed — only the coordinate arrays
@@ -364,36 +505,14 @@ def open_subdomains(
         grouped_by_time[first_time_key], use_cftime, pattern
     )
 
-    shape = template.shape  # e.g. (1, nlayer, nrow, ncol)
-    dims = template.dims  # e.g. ("time", "layer", "y", "x")
-    dtype = template.dtype
-    if has_time:
-        time_axis = list(dims).index("time")
-    else:
-        time_axis = -1  # steady-state, no time dimension
-
-    # Delayed tasks for each timestep, which will be concatenated into a single
-    # dask array. One delayed task per timestep → outer graph depth 3, O(n_time)
-    # tasks
-    merged = []
-    for time_key in raw_times_sorted:
-        group = grouped_by_time[time_key]
-        timestep_data = dask.delayed(_merge_subdomains_values)(
-            group, use_cftime, pattern
-        )
-        merged.append(dask.array.from_delayed(timestep_data, shape=shape, dtype=dtype))
-    data = dask.array.concatenate(merged, axis=time_axis)
-
-    # Build the full time coordinate
-    coords = dict(template.coords)
-    if has_time:
-        if use_cftime:
-            time_coord = xr.CFTimeIndex(converted_times)
-        else:
-            time_coord = np.array(converted_times, dtype="datetime64[ns]")
-        coords["time"] = time_coord
-
-    return xr.DataArray(data, coords, dims, name=template.name, attrs=template.attrs)
+    return _open_idf_chunked_by_time(
+        _merge_subdomains_values,
+        grouped_by_time,
+        template,
+        has_time,
+        use_cftime,
+        pattern,
+    )
 
 
 def open_dataset(globpath, use_cftime=False, pattern=None):

From 0bc5815b4b7607d6a426e3757c9f2d191ab4baa5 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Thu, 1 Oct 2026 17:27:37 +0200
Subject: [PATCH 02/27] Start making test and make sure it works

---
 imod/formats/array_io/reading.py    |  4 ++--
 imod/formats/idf.py                 |  8 ++++----
 imod/tests/test_formats/test_idf.py | 15 +++++++++++++++
 3 files changed, 21 insertions(+), 6 deletions(-)

diff --git a/imod/formats/array_io/reading.py b/imod/formats/array_io/reading.py
index ddde13c4e..aa23ccd54 100644
--- a/imod/formats/array_io/reading.py
+++ b/imod/formats/array_io/reading.py
@@ -280,7 +280,7 @@ def _dask(path, attrs=None, pattern=None, _read=None, header=None):
     return x, attrs
 
 
-def _load_as_components(paths, use_cftime, _read, headers):
+def load_as_components(paths, use_cftime, _read, headers):
     """Combine a list of paths to IDFs to a single xarray.DataArray and return as components"""
     # this function also works for single IDFs
     names = [h["name"] for h in headers]
@@ -324,7 +324,7 @@ def _load_as_components(paths, use_cftime, _read, headers):
 
 
 def _load(paths, use_cftime, _read, headers):
-    dask_array, coords, dims, name = _load_as_components(
+    dask_array, coords, dims, name = load_as_components(
         paths, use_cftime, _read, headers
     )
     out = xr.DataArray(dask_array, coords, dims, name=name)
diff --git a/imod/formats/idf.py b/imod/formats/idf.py
index c0443e8b6..1cef59eec 100644
--- a/imod/formats/idf.py
+++ b/imod/formats/idf.py
@@ -20,7 +20,7 @@
 import xarray as xr
 
 import imod
-from imod.formats import _load_as_components, array_io
+from imod.formats import array_io
 from imod.typing.structured import merge_partitions_as_da_components
 
 # Make sure we can still use the built-in function...
@@ -298,8 +298,8 @@ def _open_idf_data(
 ):
     """Open IDF files and return as a DataArray."""
     headers = [header(p, pattern) for p in paths]
-    data, _, _, _ = _load_as_components(
-        paths, use_cftime=use_cftime, _read=array_io.read, headers=headers
+    data, _, _, _ = array_io.reading.load_as_components(
+        paths, use_cftime=use_cftime, _read=_read, headers=headers
     )
     return data
 
@@ -377,7 +377,7 @@ def _open_idf_chunked_by_time(
     return xr.DataArray(data, coords, dims, name=template.name, attrs=template.attrs)
 
 
-def open_idf_by_time(path, use_cftime=False, pattern=None):
+def open_by_time(path, use_cftime=False, pattern=None):
     """
     Open IDF files grouped by time.
 
diff --git a/imod/tests/test_formats/test_idf.py b/imod/tests/test_formats/test_idf.py
index af19ad894..8cada085c 100644
--- a/imod/tests/test_formats/test_idf.py
+++ b/imod/tests/test_formats/test_idf.py
@@ -299,6 +299,21 @@ def test_open_subdomains_error(subdomains, expected, equidistant, tmp_path):
         idf.open_subdomains(tmp_path / "subdomains_*.idf")
 
 
+@parametrize_with_cases(
+    "subdomains,expected,equidistant", cases=SubdomainCases, has_tag="no_species"
+)
+def test_open_by_time(subdomains, expected, equidistant, tmp_path):
+    _save_subdomains_no_species(subdomains, tmp_path)
+
+    # Test with pattern
+    pattern = r"{name}_{time}_l{layer}_p000"
+
+    da = idf.open_by_time(tmp_path / "subdomains_*_p000.idf", pattern=pattern).load()
+
+    assert da.dims == ("time", "layer", "y", "x")
+    assert da.name == "subdomains"
+
+
 def test_xycoords_equidistant():
     dx, dy = 1.0, -1.0
     xmin, xmax = 0.0, 4.0

From 03207043821e28efe5e2c43d3e702c86a9b5ef34 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Thu, 1 Oct 2026 17:33:56 +0200
Subject: [PATCH 03/27] Add chunks assertion

---
 imod/tests/test_formats/test_idf.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/imod/tests/test_formats/test_idf.py b/imod/tests/test_formats/test_idf.py
index 8cada085c..31586ba9e 100644
--- a/imod/tests/test_formats/test_idf.py
+++ b/imod/tests/test_formats/test_idf.py
@@ -308,10 +308,11 @@ def test_open_by_time(subdomains, expected, equidistant, tmp_path):
     # Test with pattern
     pattern = r"{name}_{time}_l{layer}_p000"
 
-    da = idf.open_by_time(tmp_path / "subdomains_*_p000.idf", pattern=pattern).load()
+    da = idf.open_by_time(tmp_path / "subdomains_*_p000.idf", pattern=pattern)
 
     assert da.dims == ("time", "layer", "y", "x")
     assert da.name == "subdomains"
+    assert da.chunks == ((1,), (da.sizes["layer"],), (da.sizes["y"],), (da.sizes["x"],))
 
 
 def test_xycoords_equidistant():

From 27d3c7d32b73aa6a76e53bcb61da39a855a88a70 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Mon, 5 Oct 2026 16:13:03 +0200
Subject: [PATCH 04/27] Improve docstring load_as_components

---
 imod/formats/array_io/reading.py | 27 ++++++++++++++++++++++++++-
 1 file changed, 26 insertions(+), 1 deletion(-)

diff --git a/imod/formats/array_io/reading.py b/imod/formats/array_io/reading.py
index aa23ccd54..8f4ef3bc5 100644
--- a/imod/formats/array_io/reading.py
+++ b/imod/formats/array_io/reading.py
@@ -281,7 +281,32 @@ def _dask(path, attrs=None, pattern=None, _read=None, header=None):
 
 
 def load_as_components(paths, use_cftime, _read, headers):
-    """Combine a list of paths to IDFs to a single xarray.DataArray and return as components"""
+    """
+    Combine a list of paths to IDFs to components that can be used to construct
+    a single xarray.DataArray.
+
+    Parameters
+    ----------
+    paths : list[str]
+        List of file paths to IDF files.
+    use_cftime : bool
+        Whether to use cftime for time coordinates.
+    _read : Callable
+        Function to read individual IDF files.
+    headers : list[dict]
+        List of headers corresponding to each IDF file.
+
+    Returns
+    -------
+    dask.array
+        Dask array containing the combined data from all IDF files.
+    dict
+        Dictionary of coordinates for the DataArray.
+    list[str]
+        List of dimension names for the DataArray.
+    str
+        Name of the DataArray.
+    """
     # this function also works for single IDFs
     names = [h["name"] for h in headers]
     _all_equal(names, "names")

From 14d4188a4d4051358870f87726b431db09a5906e Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Mon, 5 Oct 2026 16:38:22 +0200
Subject: [PATCH 05/27] Update faulty docstrings

---
 imod/formats/idf.py | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/imod/formats/idf.py b/imod/formats/idf.py
index 1cef59eec..78bd31578 100644
--- a/imod/formats/idf.py
+++ b/imod/formats/idf.py
@@ -296,7 +296,7 @@ def _open_idf_data(
     use_cftime: bool,
     pattern: str | Pattern,
 ):
-    """Open IDF files and return as a DataArray."""
+    """Open IDF files and return as numpy array."""
     headers = [header(p, pattern) for p in paths]
     data, _, _, _ = array_io.reading.load_as_components(
         paths, use_cftime=use_cftime, _read=_read, headers=headers
@@ -388,11 +388,10 @@ def open_by_time(path, use_cftime=False, pattern=None):
     use_cftime : bool, optional
         Whether to use cftime for time coordinates.
     pattern : str, regex pattern, optional
-        If no pattern is provided, the function will first try:
-        "{name}_c{species}_{time}_l{layer}_p{subdomain}"
-        and if that fails:
-        "{name}_{time}_l{layer}_p{subdomain}"
-        Following the iMOD5/iMOD-WQ filename conventions.
+        If the filenames do match default naming conventions of
+        {name}_{time}_l{layer}, a custom pattern can be defined here either
+        as a string, or as a compiled regular expression pattern. See the
+        examples below.
 
     Returns
     -------

From ac781253ec37246087eee5751f194bf373e7cb30 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Mon, 5 Oct 2026 16:38:37 +0200
Subject: [PATCH 06/27] Finish unittests for open_by_time

---
 imod/tests/test_formats/test_idf.py | 76 ++++++++++++++++++++++++++---
 1 file changed, 69 insertions(+), 7 deletions(-)

diff --git a/imod/tests/test_formats/test_idf.py b/imod/tests/test_formats/test_idf.py
index 31586ba9e..db254cf0e 100644
--- a/imod/tests/test_formats/test_idf.py
+++ b/imod/tests/test_formats/test_idf.py
@@ -1,6 +1,7 @@
 import datetime
 
 import numpy as np
+import pandas as pd
 import pytest
 import xarray as xr
 from pytest import approx
@@ -299,20 +300,81 @@ def test_open_subdomains_error(subdomains, expected, equidistant, tmp_path):
         idf.open_subdomains(tmp_path / "subdomains_*.idf")
 
 
+class TemporalCases:
+    def create_da(self, ntime):
+        nlayer, nrow, ncol = 3, 6, 8
+        dx = 1.0
+        dy = -1.0
+        xmin, xmax = 0.0, 8.0
+        ymin, ymax = 0.0, 6.0
+        layer = [1, 2, 3]
+
+        time = pd.date_range("2000-01-01", periods=ntime, freq="D")
+
+        kwargs = {"name": "temporal_data", "dims": ("time", "layer", "y", "x")}
+        kwargs["coords"] = util.spatial._xycoords((xmin, xmax, ymin, ymax), (dx, dy))
+        kwargs["coords"]["layer"] = layer
+        kwargs["coords"]["time"] = time
+        kwargs["data"] = np.ones((ntime, nlayer, nrow, ncol), dtype=np.float64)
+        return xr.DataArray(**kwargs).cumsum(dim="time")
+
+    def case_single_time(self):
+        ntime = 1
+        return self.create_da(ntime=ntime), ntime
+
+    def case_multiple_times(self):
+        ntime = 5
+        return self.create_da(ntime=ntime), ntime
+
+
 @parametrize_with_cases(
-    "subdomains,expected,equidistant", cases=SubdomainCases, has_tag="no_species"
+    "temporal_data,ntime", cases=TemporalCases
 )
-def test_open_by_time(subdomains, expected, equidistant, tmp_path):
-    _save_subdomains_no_species(subdomains, tmp_path)
+def test_open_by_time__with_pattern(temporal_data, ntime, tmp_path):
+    idf.save(tmp_path / "temporal_data", temporal_data)
 
     # Test with pattern
-    pattern = r"{name}_{time}_l{layer}_p000"
+    pattern = r"{name}_{time}_l{layer}"
 
-    da = idf.open_by_time(tmp_path / "subdomains_*_p000.idf", pattern=pattern)
+    da = idf.open_by_time(tmp_path / "temporal_data_*.idf", pattern=pattern)
 
     assert da.dims == ("time", "layer", "y", "x")
-    assert da.name == "subdomains"
-    assert da.chunks == ((1,), (da.sizes["layer"],), (da.sizes["y"],), (da.sizes["x"],))
+    assert da.name == "temporal_data"
+    assert da.sizes["layer"] == 3
+    assert da.sizes["y"] == 6
+    assert da.sizes["x"] == 8
+    assert da.sizes["time"] == ntime
+    time_chunk_shape = (1,) * ntime
+    assert da.chunks == (time_chunk_shape, (da.sizes["layer"],), (da.sizes["y"],), (da.sizes["x"],))
+
+    # Compute and see if no error is thrown
+    da = da.compute()
+
+    np.testing.assert_allclose(da, temporal_data)
+
+
+@parametrize_with_cases(
+    "temporal_data,ntime", cases=TemporalCases
+)
+def test_open_by_time__without_pattern(temporal_data, ntime, tmp_path):
+    idf.save(tmp_path / "temporal_data", temporal_data)
+
+    # Test without pattern
+    da = idf.open_by_time(tmp_path / "temporal_data_*.idf")
+
+    assert da.dims == ("time", "layer", "y", "x")
+    assert da.name == "temporal_data"
+    assert da.sizes["layer"] == 3
+    assert da.sizes["y"] == 6
+    assert da.sizes["x"] == 8
+    assert da.sizes["time"] == ntime
+    time_chunk_shape = (1,) * ntime
+    assert da.chunks == (time_chunk_shape, (da.sizes["layer"],), (da.sizes["y"],), (da.sizes["x"],))
+
+    # Compute and see if no error is thrown
+    da = da.compute()
+
+    np.testing.assert_allclose(da, temporal_data)
 
 
 def test_xycoords_equidistant():

From b001c4520d7a5bf5b60b118afb74623b1c431cf1 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Tue, 6 Oct 2026 07:50:36 +0200
Subject: [PATCH 07/27] Put path handling logic in separate helper function and
 use in more places in idf.py

---
 imod/formats/array_io/reading.py | 12 ++++++++++--
 imod/formats/idf.py              |  9 ++++-----
 2 files changed, 14 insertions(+), 7 deletions(-)

diff --git a/imod/formats/array_io/reading.py b/imod/formats/array_io/reading.py
index 8f4ef3bc5..893ed4afe 100644
--- a/imod/formats/array_io/reading.py
+++ b/imod/formats/array_io/reading.py
@@ -364,7 +364,7 @@ def _load(paths, use_cftime, _read, headers):
     return out
 
 
-def _open(path, use_cftime, pattern, header, _read):
+def handle_path(path: pathlib.Path | str | list[pathlib.Path]) -> list[pathlib.Path]:
     if isinstance(path, pathlib.Path):
         path = str(path)
 
@@ -373,8 +373,16 @@ def _open(path, use_cftime, pattern, header, _read):
     else:
         paths = [pathlib.Path(p) for p in glob.glob(path)]
 
-    headers = [header(p, pattern) for p in paths]
     n = len(paths)
     if n == 0:
         raise FileNotFoundError(f"Could not find any files matching {path}")
+
+    return paths
+
+
+def _open(path, use_cftime, pattern, header, _read):
+    paths = handle_path(path)
+
+    headers = [header(p, pattern) for p in paths]
+
     return _load(paths, use_cftime, _read, headers)
diff --git a/imod/formats/idf.py b/imod/formats/idf.py
index 78bd31578..d771e89de 100644
--- a/imod/formats/idf.py
+++ b/imod/formats/idf.py
@@ -399,7 +399,7 @@ def open_by_time(path, use_cftime=False, pattern=None):
 
     """
 
-    paths = sorted(glob.glob(str(path)))
+    paths = array_io.reading.handle_path(path)
     parsed = [imod.util.path.decompose(path, pattern) for path in paths]
 
     has_time = "time" in parsed[0]
@@ -442,8 +442,8 @@ def open_subdomains(
 
     Parameters
     ----------
-    path : str or Path
-        Global path.
+    path : str, Path or list
+        Global path or list of paths to open.
     use_cftime : bool, optional
     pattern : str, regex pattern, optional
         If no pattern is provided, the function will first try:
@@ -461,8 +461,7 @@ def open_subdomains(
     # merges the subdomains into one DataArray, in a delayed manner, chunked per
     # timestep. A lot of logic in this function is about grouping the files by
     # subdomain and time, and setting the right time coordinate again.
-
-    paths = sorted(glob.glob(str(path)))
+    paths = array_io.reading.handle_path(path)
 
     if pattern is None:
         # If no pattern provided test if

From 0a80190bad8c506aa3fac6c56038a2128b50f888 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Tue, 6 Oct 2026 14:29:26 +0200
Subject: [PATCH 08/27] Add docstring

---
 imod/formats/array_io/reading.py | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/imod/formats/array_io/reading.py b/imod/formats/array_io/reading.py
index 893ed4afe..fccb6e8ad 100644
--- a/imod/formats/array_io/reading.py
+++ b/imod/formats/array_io/reading.py
@@ -365,6 +365,24 @@ def _load(paths, use_cftime, _read, headers):
 
 
 def handle_path(path: pathlib.Path | str | list[pathlib.Path]) -> list[pathlib.Path]:
+    """
+    Handle a path input and return a list of pathlib.Path objects.
+
+    Parameters
+    ----------
+    path : pathlib.Path, str, or list of pathlib.Path
+        The input path(s) to handle.
+
+    Returns
+    -------
+    list[pathlib.Path]
+        List of pathlib.Path objects corresponding to the input path(s).
+
+    Raises
+    ------
+    FileNotFoundError
+        If no files matching the input path(s) are found.
+    """
     if isinstance(path, pathlib.Path):
         path = str(path)
 

From 025a73be5a9a506d37c10de8a75c5d0ef4b4b3e0 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Tue, 6 Oct 2026 14:30:43 +0200
Subject: [PATCH 09/27] Chunk idf data by time

---
 imod/formats/idf.py     | 93 ++++++++++++++++++++++++++---------------
 imod/formats/prj/prj.py |  4 +-
 2 files changed, 62 insertions(+), 35 deletions(-)

diff --git a/imod/formats/idf.py b/imod/formats/idf.py
index d771e89de..94fdba984 100644
--- a/imod/formats/idf.py
+++ b/imod/formats/idf.py
@@ -27,7 +27,7 @@
 f_open = open
 
 
-def header(path, pattern):
+def header(path: str | pathlib.Path, pattern: Pattern | str) -> dict[str, Any]:
     """Read the IDF header information into a dictionary"""
     attrs = imod.util.path.decompose(path, pattern)
     with f_open(path, "rb") as f:
@@ -290,27 +290,42 @@ def check_subdomain_consistency(
             f"Each subdomain must have the same number of IDF files, found: {n_idf_per_subdomain}"
         )
 
-
-def _open_idf_data(
-    paths: list[str],
+def _open_idf_data_from_group(
+    group: tuple[list[str], list[dict[str, Any]]],
     use_cftime: bool,
-    pattern: str | Pattern,
 ):
-    """Open IDF files and return as numpy array."""
-    headers = [header(p, pattern) for p in paths]
+    """Open IDF files from a group and return as numpy array."""
+    group_paths, group_headers = group
+    # load_as_components mutates header dictionaries while preparing lazy reads.
+    # Copy to keep repeated/template reads safe.
+    headers = [h.copy() for h in group_headers]
     data, _, _, _ = array_io.reading.load_as_components(
-        paths, use_cftime=use_cftime, _read=_read, headers=headers
+        group_paths, use_cftime=use_cftime, _read=_read, headers=headers
     )
     return data
 
 
+def _open_idf_from_group(
+    group: tuple[list[str], list[dict[str, Any]]],
+    use_cftime: bool,
+) -> xr.DataArray:
+    """Open IDF files and return a DataArray template with header-derived dims/coords."""
+    # load_as_components mutates header dictionaries while preparing lazy reads.
+    group_paths, group_headers = group
+    group_headers = [h.copy() for h in group_headers]
+    data, coords, dims, name = array_io.reading.load_as_components(
+        group_paths, use_cftime=use_cftime, _read=_read, headers=group_headers
+    )
+    return xr.DataArray(data, coords, dims, name=name)
+
+
 def _open_idf_chunked_by_time(
     open_func: Callable,
     grouped_by_time: dict[Any, Any],
     template: xr.DataArray,
     has_time: bool,
     use_cftime: bool,
-    pattern: str | Pattern,
+    *args,
 ):
     """
     Open IDF files chunked by time using a template DataArray.
@@ -328,8 +343,8 @@ def _open_idf_chunked_by_time(
         Whether the data has a time dimension.
     use_cftime : bool
         Whether to use cftime for time coordinates.
-    pattern : str | Pattern
-        Pattern to match file paths.
+    *args : tuple
+        Additional arguments to pass to the open_func.
 
     Returns
     -------
@@ -359,7 +374,7 @@ def _open_idf_chunked_by_time(
     dask_arrays = []
     for time_key in raw_times_sorted:
         group = grouped_by_time[time_key]
-        timestep_data = dask.delayed(open_func)(group, use_cftime, pattern)
+        timestep_data = dask.delayed(open_func)(group, use_cftime, *args)
         dask_arrays.append(
             dask.array.from_delayed(timestep_data, shape=shape, dtype=dtype)
         )
@@ -377,7 +392,7 @@ def _open_idf_chunked_by_time(
     return xr.DataArray(data, coords, dims, name=template.name, attrs=template.attrs)
 
 
-def open_by_time(path, use_cftime=False, pattern=None):
+def open_by_time(path, use_cftime=False, pattern=None, headers=None):
     """
     Open IDF files grouped by time.
 
@@ -392,45 +407,55 @@ def open_by_time(path, use_cftime=False, pattern=None):
         {name}_{time}_l{layer}, a custom pattern can be defined here either
         as a string, or as a compiled regular expression pattern. See the
         examples below.
+    headers:  list, optional
+        List of headers. Can be provided to avoid re-reading headers from the
+        files.
 
     Returns
     -------
     xarray.DataArray
 
     """
-
     paths = array_io.reading.handle_path(path)
-    parsed = [imod.util.path.decompose(path, pattern) for path in paths]
+    if headers is None:
+        headers = [header(p, pattern) for p in paths]
 
-    has_time = "time" in parsed[0]
-
-    # Group by time (datetime.datetime from decompose), then by subdomain.
-    # Each delayed task processes one timestep, keeping the outer graph at O(n_time).
-    grouped_by_time: DefaultDict[Any, list] = defaultdict(list)
+    # If headers were supplied by the caller, they may already contain
+    # dimensions like time/layer that are absent in the filename itself.
+    # Prefer header metadata for grouping and template construction.
+    has_time = "time" in headers[0].get("dims", [])
 
-    for match, p in zip(parsed, paths):
+    grouped_paths_by_time: DefaultDict[Any, list] = defaultdict(list)
+    grouped_headers_by_time: DefaultDict[Any, list] = defaultdict(list)
+    for p, h in zip(paths, headers):
         if has_time:
-            time_key = match["time"]
+            time_key = h["time"]
         else:
             # Work around for files without time dimension
             # (imod.util.time._convert_datetimes special-cases this string)
             time_key = "steady-state"
-        grouped_by_time[time_key].append(p)
+        grouped_paths_by_time[time_key].append(p)
+        grouped_headers_by_time[time_key].append(h)
 
-    # Sort and convert times before calling idf.open so that
-    # use_cftime is already correct when the template is built.
-    raw_times_sorted = sorted(grouped_by_time.keys())
+    raw_times_sorted = sorted(grouped_paths_by_time.keys())
+
+    grouped_by_time: dict[Any, tuple[list[str], list[dict[str, Any]]]] = {
+        key: (grouped_paths_by_time[key], grouped_headers_by_time[key])
+        for key in raw_times_sorted
+    }
 
-    # Call idf.open eagerly for the first timestep to obtain a
-    # coordinate template. No data is computed — only the coordinate arrays
-    # (which are numpy) are used; the dask data array is discarded.
     first_time_key = raw_times_sorted[0]
-    template = open(
-        grouped_by_time[first_time_key], use_cftime=use_cftime, pattern=pattern
+    template = _open_idf_from_group(
+        grouped_by_time[first_time_key],
+        use_cftime=use_cftime,
     )
 
     return _open_idf_chunked_by_time(
-        _open_idf_data, grouped_by_time, template, has_time, use_cftime, pattern
+        _open_idf_data_from_group,
+        grouped_by_time,
+        template,
+        has_time,
+        use_cftime,
     )
 
 
@@ -502,6 +527,8 @@ def open_subdomains(
     template = _merge_subdomains_to_dataarray(
         grouped_by_time[first_time_key], use_cftime, pattern
     )
+    # Collect additional arguments to pass to the open function.
+    args_for_func = (pattern,)
 
     return _open_idf_chunked_by_time(
         _merge_subdomains_values,
@@ -509,7 +536,7 @@ def open_subdomains(
         template,
         has_time,
         use_cftime,
-        pattern,
+        *args_for_func,
     )
 
 
diff --git a/imod/formats/prj/prj.py b/imod/formats/prj/prj.py
index dabaf2478..96c136c41 100644
--- a/imod/formats/prj/prj.py
+++ b/imod/formats/prj/prj.py
@@ -568,10 +568,10 @@ def _create_dataarray_from_paths(
     factor = _get_array_transformation_parameters(headers, "factor", dim)
     addition = _get_array_transformation_parameters(headers, "addition", dim)
     da = _try_read_with_func(
-        imod.formats.array_io.reading._load,
+        imod.formats.idf.open_by_time,
         paths,
         use_cftime=False,
-        _read=imod.idf._read,
+        pattern=None,
         headers=headers,
     )
 

From 0dfa390d62abbfdc1772650aa28b1c25c784ece2 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Tue, 6 Oct 2026 14:30:50 +0200
Subject: [PATCH 10/27] Format

---
 imod/tests/test_formats/test_idf.py | 22 ++++++++++++++--------
 1 file changed, 14 insertions(+), 8 deletions(-)

diff --git a/imod/tests/test_formats/test_idf.py b/imod/tests/test_formats/test_idf.py
index db254cf0e..86aff53ec 100644
--- a/imod/tests/test_formats/test_idf.py
+++ b/imod/tests/test_formats/test_idf.py
@@ -327,9 +327,7 @@ def case_multiple_times(self):
         return self.create_da(ntime=ntime), ntime
 
 
-@parametrize_with_cases(
-    "temporal_data,ntime", cases=TemporalCases
-)
+@parametrize_with_cases("temporal_data,ntime", cases=TemporalCases)
 def test_open_by_time__with_pattern(temporal_data, ntime, tmp_path):
     idf.save(tmp_path / "temporal_data", temporal_data)
 
@@ -345,7 +343,12 @@ def test_open_by_time__with_pattern(temporal_data, ntime, tmp_path):
     assert da.sizes["x"] == 8
     assert da.sizes["time"] == ntime
     time_chunk_shape = (1,) * ntime
-    assert da.chunks == (time_chunk_shape, (da.sizes["layer"],), (da.sizes["y"],), (da.sizes["x"],))
+    assert da.chunks == (
+        time_chunk_shape,
+        (da.sizes["layer"],),
+        (da.sizes["y"],),
+        (da.sizes["x"],),
+    )
 
     # Compute and see if no error is thrown
     da = da.compute()
@@ -353,9 +356,7 @@ def test_open_by_time__with_pattern(temporal_data, ntime, tmp_path):
     np.testing.assert_allclose(da, temporal_data)
 
 
-@parametrize_with_cases(
-    "temporal_data,ntime", cases=TemporalCases
-)
+@parametrize_with_cases("temporal_data,ntime", cases=TemporalCases)
 def test_open_by_time__without_pattern(temporal_data, ntime, tmp_path):
     idf.save(tmp_path / "temporal_data", temporal_data)
 
@@ -369,7 +370,12 @@ def test_open_by_time__without_pattern(temporal_data, ntime, tmp_path):
     assert da.sizes["x"] == 8
     assert da.sizes["time"] == ntime
     time_chunk_shape = (1,) * ntime
-    assert da.chunks == (time_chunk_shape, (da.sizes["layer"],), (da.sizes["y"],), (da.sizes["x"],))
+    assert da.chunks == (
+        time_chunk_shape,
+        (da.sizes["layer"],),
+        (da.sizes["y"],),
+        (da.sizes["x"],),
+    )
 
     # Compute and see if no error is thrown
     da = da.compute()

From 3636410351b909ad9451a5ccdf93ad2db7099253 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Tue, 6 Oct 2026 14:37:40 +0200
Subject: [PATCH 11/27] Add extra assert statments to verify shape and
 chunksizes

---
 imod/tests/test_formats/test_prj.py | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/imod/tests/test_formats/test_prj.py b/imod/tests/test_formats/test_prj.py
index e98fe913c..7cb08aed3 100644
--- a/imod/tests/test_formats/test_prj.py
+++ b/imod/tests/test_formats/test_prj.py
@@ -593,6 +593,26 @@ def test_open_projectfile_data(self, projectfile, request):
         assert isinstance(content["pcg"], dict)
         assert set(repeats["rch"]) == {datetime(1899, 4, 1), datetime(1899, 10, 1)}
 
+        # Test if chunking of topsystem IDF data went correctly, should be
+        # chunked per timestep and have only 1 layer.
+        assert content["ghb"]["conductance"].shape == (3, 1, 2, 2)
+        assert content["ghb"]["conductance"].chunksizes["time"] == (1, 1, 1)
+        assert content["ghb"]["conductance"].chunksizes["layer"] == (1,)
+        assert content["ghb"]["conductance"].chunksizes["x"] == (2,)
+        assert content["ghb"]["conductance"].chunksizes["y"] == (2,)
+        # Load the data to ensure that dask arrays are correctly formed and
+        # can be evaluated.
+        content["ghb"]["conductance"].load()
+    
+        # Test if chunking BND data is expected, should be one chunk including
+        # the two layers.
+        assert content["bnd"]["ibound"].shape == (2, 2, 2)
+        assert content["bnd"]["ibound"].chunksizes["layer"] == (2,)
+        assert content["bnd"]["ibound"].chunksizes["x"] == (2,)
+        assert content["bnd"]["ibound"].chunksizes["y"] == (2,)
+        content["bnd"]["ibound"].load()
+
+
     def test_open_projectfile_data__faulty_well(self, projectfile):
         basepath = self.basepath
         # Setup faulty well

From b9b4f31e53b74e87cf1842bda5e4cc86072f9635 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Tue, 6 Oct 2026 14:38:09 +0200
Subject: [PATCH 12/27] Format

---
 imod/formats/idf.py                 | 1 +
 imod/tests/test_formats/test_prj.py | 3 +--
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/imod/formats/idf.py b/imod/formats/idf.py
index 94fdba984..a1013baf6 100644
--- a/imod/formats/idf.py
+++ b/imod/formats/idf.py
@@ -290,6 +290,7 @@ def check_subdomain_consistency(
             f"Each subdomain must have the same number of IDF files, found: {n_idf_per_subdomain}"
         )
 
+
 def _open_idf_data_from_group(
     group: tuple[list[str], list[dict[str, Any]]],
     use_cftime: bool,
diff --git a/imod/tests/test_formats/test_prj.py b/imod/tests/test_formats/test_prj.py
index 7cb08aed3..1bca91c17 100644
--- a/imod/tests/test_formats/test_prj.py
+++ b/imod/tests/test_formats/test_prj.py
@@ -603,7 +603,7 @@ def test_open_projectfile_data(self, projectfile, request):
         # Load the data to ensure that dask arrays are correctly formed and
         # can be evaluated.
         content["ghb"]["conductance"].load()
-    
+
         # Test if chunking BND data is expected, should be one chunk including
         # the two layers.
         assert content["bnd"]["ibound"].shape == (2, 2, 2)
@@ -612,7 +612,6 @@ def test_open_projectfile_data(self, projectfile, request):
         assert content["bnd"]["ibound"].chunksizes["y"] == (2,)
         content["bnd"]["ibound"].load()
 
-
     def test_open_projectfile_data__faulty_well(self, projectfile):
         basepath = self.basepath
         # Setup faulty well

From 6c22705e7064ba6591791546bed8c77e631d1b0c Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Tue, 6 Oct 2026 14:47:54 +0200
Subject: [PATCH 13/27] Rename function to idf.open_chunked_by_time

---
 imod/formats/idf.py                 | 2 +-
 imod/tests/test_formats/test_idf.py | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/imod/formats/idf.py b/imod/formats/idf.py
index a1013baf6..5a3199229 100644
--- a/imod/formats/idf.py
+++ b/imod/formats/idf.py
@@ -393,7 +393,7 @@ def _open_idf_chunked_by_time(
     return xr.DataArray(data, coords, dims, name=template.name, attrs=template.attrs)
 
 
-def open_by_time(path, use_cftime=False, pattern=None, headers=None):
+def open_chunked_by_time(path, use_cftime=False, pattern=None, headers=None):
     """
     Open IDF files grouped by time.
 
diff --git a/imod/tests/test_formats/test_idf.py b/imod/tests/test_formats/test_idf.py
index 86aff53ec..93afe694e 100644
--- a/imod/tests/test_formats/test_idf.py
+++ b/imod/tests/test_formats/test_idf.py
@@ -334,7 +334,7 @@ def test_open_by_time__with_pattern(temporal_data, ntime, tmp_path):
     # Test with pattern
     pattern = r"{name}_{time}_l{layer}"
 
-    da = idf.open_by_time(tmp_path / "temporal_data_*.idf", pattern=pattern)
+    da = idf.open_chunked_by_time(tmp_path / "temporal_data_*.idf", pattern=pattern)
 
     assert da.dims == ("time", "layer", "y", "x")
     assert da.name == "temporal_data"
@@ -361,7 +361,7 @@ def test_open_by_time__without_pattern(temporal_data, ntime, tmp_path):
     idf.save(tmp_path / "temporal_data", temporal_data)
 
     # Test without pattern
-    da = idf.open_by_time(tmp_path / "temporal_data_*.idf")
+    da = idf.open_chunked_by_time(tmp_path / "temporal_data_*.idf")
 
     assert da.dims == ("time", "layer", "y", "x")
     assert da.name == "temporal_data"

From 2635c25c54e637339716662831af26fbfb20a546 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Tue, 6 Oct 2026 16:03:42 +0200
Subject: [PATCH 14/27] Include mypy typing for idf.py

---
 imod/formats/array_io/reading.py |  4 +-
 imod/formats/idf.py              | 97 +++++++++++++++++++++++---------
 imod/util/path.py                |  3 +-
 imod/util/time.py                |  3 +-
 pyproject.toml                   |  6 +-
 5 files changed, 81 insertions(+), 32 deletions(-)

diff --git a/imod/formats/array_io/reading.py b/imod/formats/array_io/reading.py
index fccb6e8ad..6ba82329f 100644
--- a/imod/formats/array_io/reading.py
+++ b/imod/formats/array_io/reading.py
@@ -364,7 +364,9 @@ def _load(paths, use_cftime, _read, headers):
     return out
 
 
-def handle_path(path: pathlib.Path | str | list[pathlib.Path]) -> list[pathlib.Path]:
+def handle_path(
+    path: pathlib.Path | str | list[pathlib.Path | str],
+) -> list[pathlib.Path]:
     """
     Handle a path input and return a list of pathlib.Path objects.
 
diff --git a/imod/formats/idf.py b/imod/formats/idf.py
index 5a3199229..49ee3d8c8 100644
--- a/imod/formats/idf.py
+++ b/imod/formats/idf.py
@@ -12,12 +12,13 @@
 from collections.abc import Iterable
 from pathlib import Path
 from re import Pattern
-from typing import Any, Callable, DefaultDict
+from typing import Any, Callable, DefaultDict, Sequence, TypeAlias
 
 import dask
 import dask.array
 import numpy as np
 import xarray as xr
+from numpy.typing import DTypeLike
 
 import imod
 from imod.formats import array_io
@@ -26,8 +27,11 @@
 # Make sure we can still use the built-in function...
 f_open = open
 
+PatternType: TypeAlias = Pattern[str] | str
+PathType: TypeAlias = str | Path
 
-def header(path: str | pathlib.Path, pattern: Pattern | str) -> dict[str, Any]:
+
+def header(path: PathType, pattern: PatternType | None) -> dict[str, Any]:
     """Read the IDF header information into a dictionary"""
     attrs = imod.util.path.decompose(path, pattern)
     with f_open(path, "rb") as f:
@@ -105,7 +109,14 @@ def header(path: str | pathlib.Path, pattern: Pattern | str) -> dict[str, Any]:
     return attrs
 
 
-def _read(path, headersize, nrow, ncol, nodata, dtype):
+def _read(
+    path: PathType,
+    headersize: int,
+    nrow: int,
+    ncol: int,
+    nodata: float,
+    dtype: DTypeLike,
+):
     """
     Read a single IDF file to a numpy.ndarray
 
@@ -133,7 +144,11 @@ def _read(path, headersize, nrow, ncol, nodata, dtype):
 
 
 # Open IDFs for multiple times and/or layers into one DataArray
-def open(path, use_cftime=False, pattern=None):
+def open(
+    path: PathType | list[PathType],
+    use_cftime=False,
+    pattern: PatternType | None = None,
+):
     r"""
     Open one or more IDF files as an xarray.DataArray.
 
@@ -223,9 +238,9 @@ def _more_than_one_unique_value(values: Iterable[Any]):
 
 
 def _merge_subdomains(
-    paths_per_subdomain: DefaultDict[Any, list[str]],
+    paths_per_subdomain: DefaultDict[Any, list[PathType]],
     use_cftime: bool,
-    pattern: str | Pattern,
+    pattern: PatternType | None,
 ):
     """
     Open and spatially merge all subdomain IDF files for one timestep.
@@ -244,9 +259,9 @@ def _merge_subdomains(
 
 
 def _merge_subdomains_values(
-    paths_per_subdomain: DefaultDict[Any, list[str]],
+    paths_per_subdomain: DefaultDict[Any, list[PathType]],
     use_cftime: bool,
-    pattern: str | Pattern,
+    pattern: PatternType,
 ):
     """Wraps ``_merge_subdomains`` to return just a numpy array for ``dask.array.from_delayed``."""
     data, _, _, _ = _merge_subdomains(paths_per_subdomain, use_cftime, pattern)
@@ -254,9 +269,9 @@ def _merge_subdomains_values(
 
 
 def _merge_subdomains_to_dataarray(
-    paths_per_subdomain: DefaultDict[Any, list[str]],
+    paths_per_subdomain: DefaultDict[Any, list[PathType]],
     use_cftime: bool,
-    pattern: str | Pattern,
+    pattern: PatternType,
 ) -> xr.DataArray:
     """Wraps ``_merge_subdomains`` to return a DataArray for coordinate template."""
     data, coords, dims, name = _merge_subdomains(
@@ -271,7 +286,7 @@ def _merge_subdomains_to_dataarray(
 
 
 def check_subdomain_consistency(
-    parsed: list[dict[str, Any]], paths: list[str], pattern: str | Pattern
+    parsed: list[dict[str, Any]], paths: Sequence[PathType], pattern: PatternType
 ):
     """Check that each subdomain has the same number of IDF files."""
     grouped = defaultdict(list)
@@ -307,7 +322,7 @@ def _open_idf_data_from_group(
 
 
 def _open_idf_from_group(
-    group: tuple[list[str], list[dict[str, Any]]],
+    group: tuple[list[PathType], list[dict[str, Any]]],
     use_cftime: bool,
 ) -> xr.DataArray:
     """Open IDF files and return a DataArray template with header-derived dims/coords."""
@@ -321,7 +336,7 @@ def _open_idf_from_group(
 
 
 def _open_idf_chunked_by_time(
-    open_func: Callable,
+    open_func: Callable[..., Any],
     grouped_by_time: dict[Any, Any],
     template: xr.DataArray,
     has_time: bool,
@@ -384,6 +399,7 @@ def _open_idf_chunked_by_time(
     # Build the full time coordinate
     coords = dict(template.coords)
     if has_time:
+        time_coord: xr.CFTimeIndex | np.ndarray
         if use_cftime:
             time_coord = xr.CFTimeIndex(converted_times)
         else:
@@ -393,14 +409,19 @@ def _open_idf_chunked_by_time(
     return xr.DataArray(data, coords, dims, name=template.name, attrs=template.attrs)
 
 
-def open_chunked_by_time(path, use_cftime=False, pattern=None, headers=None):
+def open_chunked_by_time(
+    path: PathType | list[PathType],
+    use_cftime: bool = False,
+    pattern: PatternType | None = None,
+    headers: list[dict[str, Any]] | None = None,
+):
     """
-    Open IDF files grouped by time.
+    Open IDF files chunked by time.
 
     Parameters
     ----------
-    path : str or Path
-        Global path.
+    path : str or Path or list of str or Path
+        Global path(s).
     use_cftime : bool, optional
         Whether to use cftime for time coordinates.
     pattern : str, regex pattern, optional
@@ -426,8 +447,8 @@ def open_chunked_by_time(path, use_cftime=False, pattern=None, headers=None):
     # Prefer header metadata for grouping and template construction.
     has_time = "time" in headers[0].get("dims", [])
 
-    grouped_paths_by_time: DefaultDict[Any, list] = defaultdict(list)
-    grouped_headers_by_time: DefaultDict[Any, list] = defaultdict(list)
+    grouped_paths_by_time: DefaultDict[Any, list[PathType]] = defaultdict(list)
+    grouped_headers_by_time: DefaultDict[Any, list[dict[str, Any]]] = defaultdict(list)
     for p, h in zip(paths, headers):
         if has_time:
             time_key = h["time"]
@@ -440,7 +461,10 @@ def open_chunked_by_time(path, use_cftime=False, pattern=None, headers=None):
 
     raw_times_sorted = sorted(grouped_paths_by_time.keys())
 
-    grouped_by_time: dict[Any, tuple[list[str], list[dict[str, Any]]]] = {
+    grouped_by_time: dict[
+        Any,
+        tuple[list[PathType], list[dict[str, Any]]],
+    ] = {
         key: (grouped_paths_by_time[key], grouped_headers_by_time[key])
         for key in raw_times_sorted
     }
@@ -461,14 +485,16 @@ def open_chunked_by_time(path, use_cftime=False, pattern=None, headers=None):
 
 
 def open_subdomains(
-    path: str | Path, use_cftime: bool = False, pattern: str | Pattern = None
+    path: PathType | list[PathType],
+    use_cftime: bool = False,
+    pattern: PatternType | None = None,
 ) -> xr.DataArray:
     """
     Combine IDF files of multiple subdomains.
 
     Parameters
     ----------
-    path : str, Path or list
+    path : PathType or list of PathType
         Global path or list of paths to open.
     use_cftime : bool, optional
     pattern : str, regex pattern, optional
@@ -504,7 +530,7 @@ def open_subdomains(
 
     # Group by time (datetime.datetime from decompose), then by subdomain.
     # Each delayed task processes one timestep, keeping the outer graph at O(n_time).
-    grouped_by_time: DefaultDict[Any, DefaultDict[Any, list]] = defaultdict(
+    grouped_by_time: DefaultDict[Any, DefaultDict[Any, list[PathType]]] = defaultdict(
         lambda: defaultdict(list)
     )
 
@@ -541,7 +567,11 @@ def open_subdomains(
     )
 
 
-def open_dataset(globpath, use_cftime=False, pattern=None):
+def open_dataset(
+    globpath: str | pathlib.Path,
+    use_cftime: bool = False,
+    pattern: PatternType | None = None,
+):
     """
     Open a set of IDFs to a dict of xarray.DataArrays.
 
@@ -590,7 +620,7 @@ def open_dataset(globpath, use_cftime=False, pattern=None):
     # group the DataArrays together using their name
     # note that directory names are ignored, and in case of duplicates, the last one wins
     names = [imod.util.path.decompose(path, pattern)["name"] for path in paths]
-    paths_by_name = {name: [] for name in np.unique(names)}
+    paths_by_name: dict[str, list[PathType]] = {name: [] for name in np.unique(names)}
     for path, name in zip(paths, names, strict=True):
         paths_by_name[name].append(path)
     # load each group into a DataArray
@@ -615,7 +645,12 @@ def open_dataset(globpath, use_cftime=False, pattern=None):
     return dataset_dict
 
 
-def write(path, a, nodata=1.0e20, dtype=np.float32):
+def write(
+    path: PathType,
+    a: xr.DataArray,
+    nodata: float = 1.0e20,
+    dtype: DTypeLike = np.float32,
+):
     """
     Write a 2D xarray.DataArray to a IDF file
 
@@ -725,7 +760,7 @@ def write(path, a, nodata=1.0e20, dtype=np.float32):
         a.values.tofile(f)
 
 
-def _as_voxeldata(a):
+def _as_voxeldata(a: xr.DataArray):
     """
     If "z" is present as a dimension, generate layer if necessary. Ensure that
     layer is the dimension (via swap_dims). Infer "dz" if necessary, and if
@@ -772,7 +807,13 @@ def _as_voxeldata(a):
     return a
 
 
-def save(path, a, nodata=1.0e20, pattern=None, dtype=np.float32):
+def save(
+    path: PathType,
+    a: xr.DataArray,
+    nodata: float = 1.0e20,
+    pattern: PatternType | None = None,
+    dtype: DTypeLike = np.float32,
+):
     """
     Write a xarray.DataArray to one or more IDF files
 
diff --git a/imod/util/path.py b/imod/util/path.py
index 23499b255..f6cacc6e4 100644
--- a/imod/util/path.py
+++ b/imod/util/path.py
@@ -8,6 +8,7 @@
 import pathlib
 import re
 import tempfile
+from re import Pattern
 from typing import Any, Optional
 
 import cftime
@@ -76,7 +77,7 @@ def _groupdict(stem: str, pattern: Optional[str | re.Pattern[str]]) -> dict[str,
     return d
 
 
-def decompose(path, pattern: Optional[str] = None) -> dict[str, Any]:
+def decompose(path, pattern: str | Pattern[str] | None = None) -> dict[str, Any]:
     r"""
     Parse a path, returning a dict of the parts, following the iMOD conventions.
 
diff --git a/imod/util/time.py b/imod/util/time.py
index 8525fcb6c..75b52bccb 100644
--- a/imod/util/time.py
+++ b/imod/util/time.py
@@ -1,5 +1,6 @@
 import datetime
 import warnings
+from typing import Any
 
 import cftime
 import dateutil
@@ -155,7 +156,7 @@ def forcing_starts_ends(package_times: np.ndarray, globaltimes: np.ndarray):
     return starts_ends
 
 
-def _convert_datetimes(times: np.ndarray, use_cftime: bool):
+def _convert_datetimes(times: list[Any], use_cftime: bool):
     """
     Return times as np.datetime64[ns] or cftime.DatetimeProlepticGregorian
     depending on whether the dates fall within the inclusive bounds of
diff --git a/pyproject.toml b/pyproject.toml
index 366b87efe..f74a2b29f 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -118,7 +118,11 @@ warn_return_any = false
 module = [
     "imod.data.*",
     "imod.evaluate.*",
-    "imod.formats.*",
+    "imod.formats.array_io.*",
+    "imod.formats.gen.*",
+    "imod.formats.prj.*",
+    "imod.formats.ipf",
+    "imod.formats.rasterio",
     "imod.schemata.*",
     "imod.templates.*",
     "imod.tests.*",

From 987ff658717eef6b3636f9446c5a3269af1a2b3d Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 7 Oct 2026 09:40:26 +0200
Subject: [PATCH 15/27] Finish faulty renaming of function

---
 imod/formats/prj/prj.py             | 2 +-
 imod/tests/test_formats/test_idf.py | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/imod/formats/prj/prj.py b/imod/formats/prj/prj.py
index 96c136c41..e5d9a7ae6 100644
--- a/imod/formats/prj/prj.py
+++ b/imod/formats/prj/prj.py
@@ -568,7 +568,7 @@ def _create_dataarray_from_paths(
     factor = _get_array_transformation_parameters(headers, "factor", dim)
     addition = _get_array_transformation_parameters(headers, "addition", dim)
     da = _try_read_with_func(
-        imod.formats.idf.open_by_time,
+        imod.formats.idf.open_chunked_by_time,
         paths,
         use_cftime=False,
         pattern=None,
diff --git a/imod/tests/test_formats/test_idf.py b/imod/tests/test_formats/test_idf.py
index 93afe694e..07b6da944 100644
--- a/imod/tests/test_formats/test_idf.py
+++ b/imod/tests/test_formats/test_idf.py
@@ -328,7 +328,7 @@ def case_multiple_times(self):
 
 
 @parametrize_with_cases("temporal_data,ntime", cases=TemporalCases)
-def test_open_by_time__with_pattern(temporal_data, ntime, tmp_path):
+def test_open_chunked_by_time__with_pattern(temporal_data, ntime, tmp_path):
     idf.save(tmp_path / "temporal_data", temporal_data)
 
     # Test with pattern
@@ -357,7 +357,7 @@ def test_open_by_time__with_pattern(temporal_data, ntime, tmp_path):
 
 
 @parametrize_with_cases("temporal_data,ntime", cases=TemporalCases)
-def test_open_by_time__without_pattern(temporal_data, ntime, tmp_path):
+def test_open_chunked_by_time__without_pattern(temporal_data, ntime, tmp_path):
     idf.save(tmp_path / "temporal_data", temporal_data)
 
     # Test without pattern

From 4b6602629b872a22dc3b7e35894fa565b469a3d5 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 7 Oct 2026 13:05:34 +0200
Subject: [PATCH 16/27] Force to string for regex pattern

---
 imod/formats/idf.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/imod/formats/idf.py b/imod/formats/idf.py
index 49ee3d8c8..c372b24d7 100644
--- a/imod/formats/idf.py
+++ b/imod/formats/idf.py
@@ -519,7 +519,7 @@ def open_subdomains(
         # If no pattern provided test if
         pattern = "{name}_c{species}_{time}_l{layer}_p{subdomain}"
         re_pattern_species = imod.util.path._custom_pattern_to_regex_pattern(pattern)
-        has_species = re_pattern_species.search(paths[0])
+        has_species = re_pattern_species.search(str(paths[0]))
         if not has_species:
             pattern = "{name}_{time}_l{layer}_p{subdomain}"
 

From 59b45c39e2040b8c4af043c4c2a915a0844d3983 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 7 Oct 2026 13:05:47 +0200
Subject: [PATCH 17/27] Add user acceptance test

---
 imod/tests/test_formats/test_idf.py | 56 +++++++++++++++++++++++++++++
 1 file changed, 56 insertions(+)

diff --git a/imod/tests/test_formats/test_idf.py b/imod/tests/test_formats/test_idf.py
index 07b6da944..e086d3ce3 100644
--- a/imod/tests/test_formats/test_idf.py
+++ b/imod/tests/test_formats/test_idf.py
@@ -638,3 +638,59 @@ def test_open_subdomains_large_scale(tmp_path):
     assert result.sizes["x"] == grid_cols * ncol
     assert result.values.dtype == np.float32
     assert np.all(result.values == 1.0)
+
+
+@pytest.mark.timeout(600, method="thread")  # 10 minutes
+@pytest.mark.user_acceptance
+def test_open_chunked_by_time_large_scale(tmp_path):
+    """
+    Test the performance and correctness of opening large-scale IDF files
+    chunked by time.
+    """
+    n_layers = 40
+    n_times = 800
+    nrow, ncol = 20, 20  # cells per subdomain
+    dx, dy = 1.0, -1.0
+
+    assert n_layers * n_times == 32_000
+
+    # Generate 50 daily timesteps starting 2000-01-01
+    start = datetime.date(2000, 1, 1)
+    date_strs = [
+        (start + datetime.timedelta(days=i)).strftime("%Y%m%d") for i in range(n_times)
+    ]
+    # Create the DataArray for each timestep
+    data = np.ones((nrow, ncol), dtype=np.float32)
+    xmin, xmax = 0, ncol * dx
+    ymin, ymax = 0, nrow * abs(dy)
+    coords = util.spatial._xycoords((xmin, xmax, ymin, ymax), (dx, dy))
+    da = xr.DataArray(data, dims=("y", "x"), coords=coords, name="head")
+
+    idf_dir = tmp_path / "idf_files"
+    idf_dir.mkdir(exist_ok=True)
+
+    # Write all 32,000 IDF files
+    for layer in range(1, n_layers + 1):
+        for date_str in date_strs:
+            idf.write(
+                idf_dir / f"head_{date_str}_l{layer}.idf",
+                da,
+            )
+
+    # Open subdomains and load into memory
+    t0 = datetime.datetime.now()
+    result = idf.open_chunked_by_time(idf_dir / "head_*.idf").load()
+    elapsed = datetime.datetime.now() - t0
+    with open(tmp_path / "open_chunked_by_time_large_scale.log", "w") as f:
+        f.write(
+            f"Elapsed time for open_chunked_by_time + load: {elapsed.total_seconds():.2f}s\n"
+        )
+
+    assert isinstance(result, xr.DataArray)
+    assert result.dims == ("time", "layer", "y", "x")
+    assert result.sizes["time"] == n_times
+    assert result.sizes["layer"] == n_layers
+    assert result.sizes["y"] == nrow
+    assert result.sizes["x"] == ncol
+    assert result.values.dtype == np.float32
+    # assert np.all(result.isel(time=slice(None, 100)) == 1.0)

From d354f0182e1ad12af3630d7294ad642f51292a1e Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 7 Oct 2026 13:26:08 +0200
Subject: [PATCH 18/27] Update changelog

---
 docs/api/changelog.rst | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst
index 88cd2a0a2..a509e7489 100644
--- a/docs/api/changelog.rst
+++ b/docs/api/changelog.rst
@@ -73,6 +73,9 @@ Changed
   (:class:`imod.mf6.River`, :class:`imod.mf6.GeneralHeadBoundary`,
   :class:`imod.mf6.Drainage`, :class:`imod.mf6.Recharge`) will also be masked
   where constant head cells are placed. 
+- :func:`imod.formats.prj.open_projectfile_data` now chunks the data by time
+  only instead of by both time and layer. This should enhance performance when
+  working with large datasets.
 
 [1.1.0] - 2026-08-03
 --------------------

From e38d98b93d817bb5e4b378a4f096a73205a9d364 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 7 Oct 2026 15:45:58 +0200
Subject: [PATCH 19/27] Call compute and remove values check that would take
 long and is not that important.

---
 imod/tests/test_formats/test_idf.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/imod/tests/test_formats/test_idf.py b/imod/tests/test_formats/test_idf.py
index e086d3ce3..5cc329d62 100644
--- a/imod/tests/test_formats/test_idf.py
+++ b/imod/tests/test_formats/test_idf.py
@@ -679,7 +679,7 @@ def test_open_chunked_by_time_large_scale(tmp_path):
 
     # Open subdomains and load into memory
     t0 = datetime.datetime.now()
-    result = idf.open_chunked_by_time(idf_dir / "head_*.idf").load()
+    result = idf.open_chunked_by_time(idf_dir / "head_*.idf").compute()
     elapsed = datetime.datetime.now() - t0
     with open(tmp_path / "open_chunked_by_time_large_scale.log", "w") as f:
         f.write(
@@ -693,4 +693,3 @@ def test_open_chunked_by_time_large_scale(tmp_path):
     assert result.sizes["y"] == nrow
     assert result.sizes["x"] == ncol
     assert result.values.dtype == np.float32
-    # assert np.all(result.isel(time=slice(None, 100)) == 1.0)

From e232bc9f83e4229de280fd887f979eac150a1a3d Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 7 Oct 2026 17:44:44 +0200
Subject: [PATCH 20/27] Enforce numpy array in delayed task

---
 imod/formats/idf.py | 14 ++++++++++++--
 1 file changed, 12 insertions(+), 2 deletions(-)

diff --git a/imod/formats/idf.py b/imod/formats/idf.py
index c372b24d7..ff513fd8e 100644
--- a/imod/formats/idf.py
+++ b/imod/formats/idf.py
@@ -262,9 +262,14 @@ def _merge_subdomains_values(
     paths_per_subdomain: DefaultDict[Any, list[PathType]],
     use_cftime: bool,
     pattern: PatternType,
-):
+) -> np.ndarray:
     """Wraps ``_merge_subdomains`` to return just a numpy array for ``dask.array.from_delayed``."""
     data, _, _, _ = _merge_subdomains(paths_per_subdomain, use_cftime, pattern)
+    # ``from_delayed`` expects each task to resolve to a concrete ndarray chunk.
+    # Returning a dask array here creates nested lazy arrays, which can leak
+    # ``dask.array.Array`` values into downstream numba-compiled code paths.
+    if isinstance(data, dask.array.Array):
+        data = data.compute()
     return data
 
 
@@ -309,7 +314,7 @@ def check_subdomain_consistency(
 def _open_idf_data_from_group(
     group: tuple[list[str], list[dict[str, Any]]],
     use_cftime: bool,
-):
+) -> np.ndarray:
     """Open IDF files from a group and return as numpy array."""
     group_paths, group_headers = group
     # load_as_components mutates header dictionaries while preparing lazy reads.
@@ -318,6 +323,11 @@ def _open_idf_data_from_group(
     data, _, _, _ = array_io.reading.load_as_components(
         group_paths, use_cftime=use_cftime, _read=_read, headers=headers
     )
+    # ``from_delayed`` expects each task to resolve to a concrete ndarray chunk.
+    # Returning a dask array here creates nested lazy arrays, which can leak
+    # ``dask.array.Array`` values into downstream numba-compiled code paths.
+    if isinstance(data, dask.array.Array):
+        data = data.compute()
     return data
 
 

From 9ece88ea3400768c9247dd66ef9791f203556f7b Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 7 Oct 2026 17:46:15 +0200
Subject: [PATCH 21/27] Add test

---
 imod/tests/test_formats/test_idf.py | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/imod/tests/test_formats/test_idf.py b/imod/tests/test_formats/test_idf.py
index 5cc329d62..7a6e5cb16 100644
--- a/imod/tests/test_formats/test_idf.py
+++ b/imod/tests/test_formats/test_idf.py
@@ -1,5 +1,7 @@
 import datetime
 
+import dask
+import dask.array
 import numpy as np
 import pandas as pd
 import pytest
@@ -383,6 +385,18 @@ def test_open_chunked_by_time__without_pattern(temporal_data, ntime, tmp_path):
     np.testing.assert_allclose(da, temporal_data)
 
 
+def test_open_chunked_by_time__delayed_chunk_returns_numpy(tmp_path):
+    temporal_data = TemporalCases().create_da(ntime=2)
+    idf.save(tmp_path / "temporal_data", temporal_data)
+
+    da = idf.open_chunked_by_time(tmp_path / "temporal_data_*.idf")
+    first_chunk_delayed = da.data.to_delayed().ravel()[0]
+    first_chunk = dask.compute(first_chunk_delayed)[0]
+
+    assert isinstance(first_chunk, np.ndarray)
+    assert not isinstance(first_chunk, dask.array.Array)
+
+
 def test_xycoords_equidistant():
     dx, dy = 1.0, -1.0
     xmin, xmax = 0.0, 4.0

From 7f8f7c6ba749b200d7d772156e32b5d849ecd321 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Mon, 12 Oct 2026 11:37:56 +0200
Subject: [PATCH 22/27] Add option to idf.open "chunk_by" to chunk by time or
 idf. idf is the default option, which keeps the old behavior.

---
 imod/formats/idf.py | 59 +++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 54 insertions(+), 5 deletions(-)

diff --git a/imod/formats/idf.py b/imod/formats/idf.py
index ff513fd8e..8f2e3eaa7 100644
--- a/imod/formats/idf.py
+++ b/imod/formats/idf.py
@@ -12,7 +12,7 @@
 from collections.abc import Iterable
 from pathlib import Path
 from re import Pattern
-from typing import Any, Callable, DefaultDict, Sequence, TypeAlias
+from typing import Any, Callable, DefaultDict, Literal, Sequence, TypeAlias
 
 import dask
 import dask.array
@@ -116,7 +116,7 @@ def _read(
     ncol: int,
     nodata: float,
     dtype: DTypeLike,
-):
+) -> np.ndarray:
     """
     Read a single IDF file to a numpy.ndarray
 
@@ -143,12 +143,12 @@ def _read(
     return array_io.reading._to_nan(a, nodata)
 
 
-# Open IDFs for multiple times and/or layers into one DataArray
 def open(
     path: PathType | list[PathType],
     use_cftime=False,
     pattern: PatternType | None = None,
-):
+    chunk_by: Literal["idf", "time"] = "idf",
+) -> xr.DataArray:
     r"""
     Open one or more IDF files as an xarray.DataArray.
 
@@ -178,6 +178,11 @@ def open(
         {name}_{time}_l{layer}, a custom pattern can be defined here either
         as a string, or as a compiled regular expression pattern. See the
         examples below.
+    chunk_by : {"idf", "time"}, optional
+        Determines how the data should be chunked when loaded into memory. "idf"
+        chunks by individual IDF files, while "time" chunks by the time
+        dimension. This can help optimize memory usage and computational
+        efficiency. The default is "idf".
 
     Returns
     -------
@@ -228,6 +233,50 @@ def open(
     generally a fiddly process. Regex notation is also impossible to
     remember. The website https://regex101.com is a nice help. Alternatively,
     the most pragmatic solution may be to just rename your files.
+
+    To increase performance when opening multiple IDF files, consider setting
+    the ``chunk_by`` parameter to "time" to control how the data is chunked in
+    memory. Data will be chunked along the time dimension, instead of per idf
+    file (so most often per time, per layer), which is the default.
+
+    >>> head = imod.idf.open("head_2001*_l*.idf", chunk_by="time")
+
+    """
+
+    if chunk_by == "idf":
+        return _open(path, use_cftime=use_cftime, pattern=pattern)
+    elif chunk_by == "time":
+        return open_chunked_by_time(
+            path, use_cftime=use_cftime, pattern=pattern, headers=None
+        )
+    else:
+        raise ValueError(
+            f"Unsupported chunk_by value: {chunk_by}, supported are 'idf' or 'time'"
+        )
+
+
+# Open IDFs for multiple times and/or layers into one DataArray
+def _open(
+    path: PathType | list[PathType],
+    use_cftime=False,
+    pattern: PatternType | None = None,
+) -> xr.DataArray:
+    """
+    Open an IDF file or list of IDF files and return as an xarray DataArray.
+
+    Parameters
+    ----------
+    path : PathType | list[PathType]
+        Path or list of paths to IDF files.
+    use_cftime : bool, optional
+        Whether to use cftime for datetime objects. Default is False.
+    pattern : PatternType | None, optional
+        Pattern to parse the filename into xarray coordinates. Default is None.
+
+    Returns
+    -------
+    xr.DataArray
+        The opened IDF data as an xarray DataArray.
     """
     return array_io.reading._open(path, use_cftime, pattern, header, _read)
 
@@ -251,7 +300,7 @@ def _merge_subdomains(
     """
     das = []
     for pathlist in paths_per_subdomain.values():
-        da = open(pathlist, use_cftime=use_cftime, pattern=pattern)
+        da = _open(pathlist, use_cftime=use_cftime, pattern=pattern)
         if "subdomain" in da.dims:
             da = da.isel(subdomain=0, drop=True)
         das.append(da)

From 1a851982cf72d5c19ace25c2885fa54476ef6130 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Mon, 12 Oct 2026 11:38:43 +0200
Subject: [PATCH 23/27] Reduce test duplication, and add test for idf.open
 option

---
 imod/tests/test_formats/test_idf.py | 48 +++++++++++------------------
 1 file changed, 18 insertions(+), 30 deletions(-)

diff --git a/imod/tests/test_formats/test_idf.py b/imod/tests/test_formats/test_idf.py
index 7a6e5cb16..cdc655adf 100644
--- a/imod/tests/test_formats/test_idf.py
+++ b/imod/tests/test_formats/test_idf.py
@@ -330,12 +330,11 @@ def case_multiple_times(self):
 
 
 @parametrize_with_cases("temporal_data,ntime", cases=TemporalCases)
-def test_open_chunked_by_time__with_pattern(temporal_data, ntime, tmp_path):
+@parametrize("pattern", [None, r"{name}_{time}_l{layer}"])
+def test_open_chunked_by_time(temporal_data, ntime, pattern, tmp_path):
     idf.save(tmp_path / "temporal_data", temporal_data)
 
     # Test with pattern
-    pattern = r"{name}_{time}_l{layer}"
-
     da = idf.open_chunked_by_time(tmp_path / "temporal_data_*.idf", pattern=pattern)
 
     assert da.dims == ("time", "layer", "y", "x")
@@ -358,34 +357,12 @@ def test_open_chunked_by_time__with_pattern(temporal_data, ntime, tmp_path):
     np.testing.assert_allclose(da, temporal_data)
 
 
-@parametrize_with_cases("temporal_data,ntime", cases=TemporalCases)
-def test_open_chunked_by_time__without_pattern(temporal_data, ntime, tmp_path):
-    idf.save(tmp_path / "temporal_data", temporal_data)
-
-    # Test without pattern
-    da = idf.open_chunked_by_time(tmp_path / "temporal_data_*.idf")
-
-    assert da.dims == ("time", "layer", "y", "x")
-    assert da.name == "temporal_data"
-    assert da.sizes["layer"] == 3
-    assert da.sizes["y"] == 6
-    assert da.sizes["x"] == 8
-    assert da.sizes["time"] == ntime
-    time_chunk_shape = (1,) * ntime
-    assert da.chunks == (
-        time_chunk_shape,
-        (da.sizes["layer"],),
-        (da.sizes["y"],),
-        (da.sizes["x"],),
-    )
-
-    # Compute and see if no error is thrown
-    da = da.compute()
-
-    np.testing.assert_allclose(da, temporal_data)
-
-
 def test_open_chunked_by_time__delayed_chunk_returns_numpy(tmp_path):
+    """
+    Test that the dask chunks returned by open_chunked_by_time is a numpy array,
+    and not a dask array, as that would mean there are nested dask arrays, which
+    are problematic.
+    """
     temporal_data = TemporalCases().create_da(ntime=2)
     idf.save(tmp_path / "temporal_data", temporal_data)
 
@@ -455,6 +432,17 @@ def test_save_topbot__single_layer(test_da, tmp_path):
     assert da_read["dz"] == approx(1.0)
 
 
+@parametrize("chunk_by", ["idf", "time"])
+def test_saveopen_layer_chunking(test_layerda, chunk_by, tmp_path):
+    da = test_layerda
+    idf.save(tmp_path / "layer", da)
+    actual = idf.open(tmp_path / "layer_l*.idf", chunk_by=chunk_by)
+    if chunk_by == "idf":
+        assert actual.chunks[0] == (1,) * actual.sizes["layer"]
+    elif chunk_by == "time":
+        assert actual.chunks[0] == (actual.sizes["layer"],)
+
+
 def test_save_topbot__layers(test_layerda, tmp_path):
     da = test_layerda
     da = da.assign_coords(z=("layer", np.arange(1.0, 6.0) - 0.5))

From 275a462b0808249706eb198729da46969f2cca53 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Mon, 12 Oct 2026 12:02:35 +0200
Subject: [PATCH 24/27] Clarify chunking behavior with better names for
 argument and values

---
 docs/api/changelog.rst              |  4 ++++
 imod/formats/idf.py                 | 25 ++++++++++++++-----------
 imod/tests/test_formats/test_idf.py | 10 +++++-----
 3 files changed, 23 insertions(+), 16 deletions(-)

diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst
index a509e7489..ab4d24974 100644
--- a/docs/api/changelog.rst
+++ b/docs/api/changelog.rst
@@ -39,6 +39,10 @@ Added
   the final package has fully empty layers trimmed off afterwards, so this
   does not affect computed values. Set to False to keep the previous
   full-layer-coordinate behaviour.
+- Added option ``chunk_layout`` to :func:`imod.idf.open` to control how the data is
+  loaded in chunks, which can improve performance and memory usage for large datasets.
+  The default option ``per_idf`` keeps the old behavior, by chunking per IDF file.
+  The alternative option ``per_timestep`` chunks per timestep.
 
 Fixed
 ~~~~~
diff --git a/imod/formats/idf.py b/imod/formats/idf.py
index 8f2e3eaa7..d669a5243 100644
--- a/imod/formats/idf.py
+++ b/imod/formats/idf.py
@@ -147,7 +147,7 @@ def open(
     path: PathType | list[PathType],
     use_cftime=False,
     pattern: PatternType | None = None,
-    chunk_by: Literal["idf", "time"] = "idf",
+    chunk_layout: Literal["per_idf", "per_timestep"] = "per_idf",
 ) -> xr.DataArray:
     r"""
     Open one or more IDF files as an xarray.DataArray.
@@ -178,11 +178,12 @@ def open(
         {name}_{time}_l{layer}, a custom pattern can be defined here either
         as a string, or as a compiled regular expression pattern. See the
         examples below.
-    chunk_by : {"idf", "time"}, optional
-        Determines how the data should be chunked when loaded into memory. "idf"
-        chunks by individual IDF files, while "time" chunks by the time
-        dimension. This can help optimize memory usage and computational
-        efficiency. The default is "idf".
+    chunk_layout : {"per_idf", "per_timestep"}, optional
+        Determines how the data should be chunked when loaded into memory.
+        ``"per_idf"`` chunks by individual IDF files, while
+        ``"per_timestep"`` chunks by timestep, grouping all layers within a
+        timestep into one chunk. This can help optimize memory usage and
+        computational efficiency. The default is ``"per_idf"``.
 
     Returns
     -------
@@ -235,23 +236,25 @@ def open(
     the most pragmatic solution may be to just rename your files.
 
     To increase performance when opening multiple IDF files, consider setting
-    the ``chunk_by`` parameter to "time" to control how the data is chunked in
+    the ``chunk_layout`` parameter to ``"per_timestep"`` to control how the
+    data is chunked in
     memory. Data will be chunked along the time dimension, instead of per idf
     file (so most often per time, per layer), which is the default.
 
-    >>> head = imod.idf.open("head_2001*_l*.idf", chunk_by="time")
+    >>> head = imod.idf.open("head_2001*_l*.idf", chunk_layout="per_timestep")
 
     """
 
-    if chunk_by == "idf":
+    if chunk_layout == "per_idf":
         return _open(path, use_cftime=use_cftime, pattern=pattern)
-    elif chunk_by == "time":
+    elif chunk_layout == "per_timestep":
         return open_chunked_by_time(
             path, use_cftime=use_cftime, pattern=pattern, headers=None
         )
     else:
         raise ValueError(
-            f"Unsupported chunk_by value: {chunk_by}, supported are 'idf' or 'time'"
+            "Unsupported chunk_layout value: "
+            f"{chunk_layout}, supported are 'per_idf' or 'per_timestep'"
         )
 
 
diff --git a/imod/tests/test_formats/test_idf.py b/imod/tests/test_formats/test_idf.py
index cdc655adf..6bcac1fd2 100644
--- a/imod/tests/test_formats/test_idf.py
+++ b/imod/tests/test_formats/test_idf.py
@@ -432,14 +432,14 @@ def test_save_topbot__single_layer(test_da, tmp_path):
     assert da_read["dz"] == approx(1.0)
 
 
-@parametrize("chunk_by", ["idf", "time"])
-def test_saveopen_layer_chunking(test_layerda, chunk_by, tmp_path):
+@parametrize("chunk_layout", ["per_idf", "per_timestep"])
+def test_saveopen_layer_chunking(test_layerda, chunk_layout, tmp_path):
     da = test_layerda
     idf.save(tmp_path / "layer", da)
-    actual = idf.open(tmp_path / "layer_l*.idf", chunk_by=chunk_by)
-    if chunk_by == "idf":
+    actual = idf.open(tmp_path / "layer_l*.idf", chunk_layout=chunk_layout)
+    if chunk_layout == "per_idf":
         assert actual.chunks[0] == (1,) * actual.sizes["layer"]
-    elif chunk_by == "time":
+    elif chunk_layout == "per_timestep":
         assert actual.chunks[0] == (actual.sizes["layer"],)
 
 

From 02218e0d5d22d14eeabea66124174bb84fe51c7e Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Mon, 12 Oct 2026 12:05:25 +0200
Subject: [PATCH 25/27] Rename open_chunked_by_time to
 open_chunked_per_timestep

---
 imod/formats/idf.py                 |  6 +++---
 imod/formats/prj/prj.py             |  2 +-
 imod/tests/test_formats/test_idf.py | 18 +++++++++---------
 3 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/imod/formats/idf.py b/imod/formats/idf.py
index d669a5243..f3b62fc3e 100644
--- a/imod/formats/idf.py
+++ b/imod/formats/idf.py
@@ -248,7 +248,7 @@ def open(
     if chunk_layout == "per_idf":
         return _open(path, use_cftime=use_cftime, pattern=pattern)
     elif chunk_layout == "per_timestep":
-        return open_chunked_by_time(
+        return open_chunked_per_timestep(
             path, use_cftime=use_cftime, pattern=pattern, headers=None
         )
     else:
@@ -471,14 +471,14 @@ def _open_idf_chunked_by_time(
     return xr.DataArray(data, coords, dims, name=template.name, attrs=template.attrs)
 
 
-def open_chunked_by_time(
+def open_chunked_per_timestep(
     path: PathType | list[PathType],
     use_cftime: bool = False,
     pattern: PatternType | None = None,
     headers: list[dict[str, Any]] | None = None,
 ):
     """
-    Open IDF files chunked by time.
+    Open IDF files chunked per timestep.
 
     Parameters
     ----------
diff --git a/imod/formats/prj/prj.py b/imod/formats/prj/prj.py
index e5d9a7ae6..6eda34e95 100644
--- a/imod/formats/prj/prj.py
+++ b/imod/formats/prj/prj.py
@@ -568,7 +568,7 @@ def _create_dataarray_from_paths(
     factor = _get_array_transformation_parameters(headers, "factor", dim)
     addition = _get_array_transformation_parameters(headers, "addition", dim)
     da = _try_read_with_func(
-        imod.formats.idf.open_chunked_by_time,
+        imod.formats.idf.open_chunked_per_timestep,
         paths,
         use_cftime=False,
         pattern=None,
diff --git a/imod/tests/test_formats/test_idf.py b/imod/tests/test_formats/test_idf.py
index 6bcac1fd2..a69275db5 100644
--- a/imod/tests/test_formats/test_idf.py
+++ b/imod/tests/test_formats/test_idf.py
@@ -331,11 +331,11 @@ def case_multiple_times(self):
 
 @parametrize_with_cases("temporal_data,ntime", cases=TemporalCases)
 @parametrize("pattern", [None, r"{name}_{time}_l{layer}"])
-def test_open_chunked_by_time(temporal_data, ntime, pattern, tmp_path):
+def test_open_chunked_per_timestep(temporal_data, ntime, pattern, tmp_path):
     idf.save(tmp_path / "temporal_data", temporal_data)
 
     # Test with pattern
-    da = idf.open_chunked_by_time(tmp_path / "temporal_data_*.idf", pattern=pattern)
+    da = idf.open_chunked_per_timestep(tmp_path / "temporal_data_*.idf", pattern=pattern)
 
     assert da.dims == ("time", "layer", "y", "x")
     assert da.name == "temporal_data"
@@ -357,16 +357,16 @@ def test_open_chunked_by_time(temporal_data, ntime, pattern, tmp_path):
     np.testing.assert_allclose(da, temporal_data)
 
 
-def test_open_chunked_by_time__delayed_chunk_returns_numpy(tmp_path):
+def test_open_chunked_per_timestep__delayed_chunk_returns_numpy(tmp_path):
     """
-    Test that the dask chunks returned by open_chunked_by_time is a numpy array,
+    Test that the dask chunks returned by open_chunked_per_timestep is a numpy array,
     and not a dask array, as that would mean there are nested dask arrays, which
     are problematic.
     """
     temporal_data = TemporalCases().create_da(ntime=2)
     idf.save(tmp_path / "temporal_data", temporal_data)
 
-    da = idf.open_chunked_by_time(tmp_path / "temporal_data_*.idf")
+    da = idf.open_chunked_per_timestep(tmp_path / "temporal_data_*.idf")
     first_chunk_delayed = da.data.to_delayed().ravel()[0]
     first_chunk = dask.compute(first_chunk_delayed)[0]
 
@@ -644,7 +644,7 @@ def test_open_subdomains_large_scale(tmp_path):
 
 @pytest.mark.timeout(600, method="thread")  # 10 minutes
 @pytest.mark.user_acceptance
-def test_open_chunked_by_time_large_scale(tmp_path):
+def test_open_chunked_per_timestep_large_scale(tmp_path):
     """
     Test the performance and correctness of opening large-scale IDF files
     chunked by time.
@@ -681,11 +681,11 @@ def test_open_chunked_by_time_large_scale(tmp_path):
 
     # Open subdomains and load into memory
     t0 = datetime.datetime.now()
-    result = idf.open_chunked_by_time(idf_dir / "head_*.idf").compute()
+    result = idf.open_chunked_per_timestep(idf_dir / "head_*.idf").compute()
     elapsed = datetime.datetime.now() - t0
-    with open(tmp_path / "open_chunked_by_time_large_scale.log", "w") as f:
+    with open(tmp_path / "open_chunked_per_timestep_large_scale.log", "w") as f:
         f.write(
-            f"Elapsed time for open_chunked_by_time + load: {elapsed.total_seconds():.2f}s\n"
+            f"Elapsed time for open_chunked_per_timestep + load: {elapsed.total_seconds():.2f}s\n"
         )
 
     assert isinstance(result, xr.DataArray)

From d73af2fa1f96154b73d71c15ac423cb8497abf3c Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Mon, 12 Oct 2026 12:07:53 +0200
Subject: [PATCH 26/27] Clarify private functions as well

---
 imod/formats/idf.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/imod/formats/idf.py b/imod/formats/idf.py
index f3b62fc3e..b797b0540 100644
--- a/imod/formats/idf.py
+++ b/imod/formats/idf.py
@@ -397,7 +397,7 @@ def _open_idf_from_group(
     return xr.DataArray(data, coords, dims, name=name)
 
 
-def _open_idf_chunked_by_time(
+def _open_idf_chunked_per_timestep(
     open_func: Callable[..., Any],
     grouped_by_time: dict[Any, Any],
     template: xr.DataArray,
@@ -406,7 +406,7 @@ def _open_idf_chunked_by_time(
     *args,
 ):
     """
-    Open IDF files chunked by time using a template DataArray.
+    Open IDF files chunked per timestep using a template DataArray.
 
     Parameters
     ----------
@@ -427,7 +427,7 @@ def _open_idf_chunked_by_time(
     Returns
     -------
     xr.DataArray
-        Combined DataArray with data from all subdomains, chunked by time.
+        Combined DataArray with data from all subdomains, chunked per timestep.
 
 
     """
@@ -537,7 +537,7 @@ def open_chunked_per_timestep(
         use_cftime=use_cftime,
     )
 
-    return _open_idf_chunked_by_time(
+    return _open_idf_chunked_per_timestep(
         _open_idf_data_from_group,
         grouped_by_time,
         template,
@@ -619,7 +619,7 @@ def open_subdomains(
     # Collect additional arguments to pass to the open function.
     args_for_func = (pattern,)
 
-    return _open_idf_chunked_by_time(
+    return _open_idf_chunked_per_timestep(
         _merge_subdomains_values,
         grouped_by_time,
         template,

From df2954716f7c276e219c32f77bc6e890d8fc72e3 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Mon, 12 Oct 2026 12:08:04 +0200
Subject: [PATCH 27/27] Format

---
 imod/tests/test_formats/test_idf.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/imod/tests/test_formats/test_idf.py b/imod/tests/test_formats/test_idf.py
index a69275db5..be8bddeac 100644
--- a/imod/tests/test_formats/test_idf.py
+++ b/imod/tests/test_formats/test_idf.py
@@ -335,7 +335,9 @@ def test_open_chunked_per_timestep(temporal_data, ntime, pattern, tmp_path):
     idf.save(tmp_path / "temporal_data", temporal_data)
 
     # Test with pattern
-    da = idf.open_chunked_per_timestep(tmp_path / "temporal_data_*.idf", pattern=pattern)
+    da = idf.open_chunked_per_timestep(
+        tmp_path / "temporal_data_*.idf", pattern=pattern
+    )
 
     assert da.dims == ("time", "layer", "y", "x")
     assert da.name == "temporal_data"