Repository: DHI/mikeio

Issue:
mikeio.read string time argument not working as intended
mikeio.read("ts.dfs0", time="2018") gives different results for equidistant and non-equidistant files. 

Reading non-equidistant files works as intended (i.e. gives all data in 2018 in the above example). If ts.dfs0 is equidistant, however, only one timestep i returned: 2018-1-1 00:00:00 !

In other cases, an error is thrown when time is a string! The below lines should give the same Dataset, but the second throws an error: 

    mikeio.read("../tests/testdata/HD2D.dfsu")["1985-08-06"]
    mikeio.read("../tests/testdata/HD2D.dfsu", time="1985-08-06")

Discussion:

Patch:
diff --git a/mikeio/__init__.py b/mikeio/__init__.py
index f4d95f3c..507c0c6b 100644
--- a/mikeio/__init__.py
+++ b/mikeio/__init__.py
@@ -38,7 +38,9 @@ from .spatial.grid_geometry import Grid1D, Grid2D, Grid3D
 from .eum import ItemInfo, EUMType, EUMUnit
 
 
-def read(filename, *, items=None, time_steps=None, time=None, **kwargs) -> Dataset:
+def read(
+    filename, *, items=None, time_steps=None, time=None, keepdims=False, **kwargs
+) -> Dataset:
     """Read data from a dfs file
 
     Parameters
@@ -63,7 +65,9 @@ def read(filename, *, items=None, time_steps=None, time=None, **kwargs) -> Datas
 
     dfs = open(filename)
 
-    return dfs.read(items=items, time_steps=time_steps, time=time, **kwargs)
+    return dfs.read(
+        items=items, time_steps=time_steps, time=time, keepdims=keepdims, **kwargs
+    )
 
 
 def open(filename: str, **kwargs):
diff --git a/mikeio/data_utils.py b/mikeio/data_utils.py
index 4aef5ec4..d032835a 100644
--- a/mikeio/data_utils.py
+++ b/mikeio/data_utils.py
@@ -24,7 +24,19 @@ class DataUtilsMixin:
     @staticmethod
     def _get_time_idx_list(time: pd.DatetimeIndex, steps):
         """Find list of idx in DatetimeIndex"""
-        # TODO: allow steps to be other DateTimeAxis
+
+        if isinstance(steps, str):
+            parts = steps.split(",")
+            if len(parts) == 1:
+                parts.append(parts[0])  # end=start
+
+            if parts[0] == "":
+                steps = slice(parts[1])  # stop only
+            elif parts[1] == "":
+                steps = slice(parts[0], None)  # start only
+            else:
+                steps = slice(parts[0], parts[1])
+
         if (isinstance(steps, Iterable) and not isinstance(steps, str)) and isinstance(
             steps[0], (str, datetime, np.datetime64, pd.Timestamp)
         ):
diff --git a/mikeio/dfs.py b/mikeio/dfs.py
index e96155f9..213d9e22 100644
--- a/mikeio/dfs.py
+++ b/mikeio/dfs.py
@@ -40,7 +40,9 @@ class _Dfs123(TimeSeries):
         self._dfs = None
         self._source = None
 
-    def read(self, *, items=None, time=None, time_steps=None) -> Dataset:
+    def read(
+        self, *, items=None, time=None, time_steps=None, keepdims=False
+    ) -> Dataset:
         """
         Read data from a dfs file
 
@@ -68,9 +70,8 @@ class _Dfs123(TimeSeries):
         item_numbers = _valid_item_numbers(self._dfs.ItemInfo, items)
         n_items = len(item_numbers)
 
-        time_steps = _valid_timesteps(self._dfs.FileInfo, time)
-        nt = len(time_steps)
-        single_time_selected = np.isscalar(time) if time is not None else False
+        single_time_selected, time_steps = _valid_timesteps(self._dfs.FileInfo, time)
+        nt = len(time_steps) if not single_time_selected else 1
 
         if self._ndim == 1:
             shape = (nt, self._nx)
@@ -79,7 +80,7 @@ class _Dfs123(TimeSeries):
         else:
             shape = (nt, self._nz, self._ny, self._nx)
 
-        if single_time_selected:
+        if single_time_selected and not keepdims:
             shape = shape[1:]
 
         data_list = [
diff --git a/mikeio/dfs0.py b/mikeio/dfs0.py
index 1cfefb60..084a3e9f 100644
--- a/mikeio/dfs0.py
+++ b/mikeio/dfs0.py
@@ -154,7 +154,7 @@ class Dfs0(TimeSeries):
 
         dfs.Close()
 
-    def read(self, items=None, time=None, time_steps=None) -> Dataset:
+    def read(self, items=None, time=None, time_steps=None, keepdims=False) -> Dataset:
         """
         Read data from a dfs0 file.
 
@@ -197,7 +197,7 @@ class Dfs0(TimeSeries):
             time_steps = range(self._n_timesteps)
         else:
             sel_time_step_str = None
-            time_steps = _valid_timesteps(dfs.FileInfo, time)
+            _, time_steps = _valid_timesteps(dfs.FileInfo, time)
 
         dfs.Close()
 
diff --git a/mikeio/dfs2.py b/mikeio/dfs2.py
index 0b10171e..3fe839f1 100644
--- a/mikeio/dfs2.py
+++ b/mikeio/dfs2.py
@@ -221,7 +221,9 @@ class Dfs2(_Dfs123):
 
         self._read_header()
 
-    def read(self, *, items=None, time=None, area=None, time_steps=None) -> Dataset:
+    def read(
+        self, *, items=None, time=None, area=None, time_steps=None, keepdims=False
+    ) -> Dataset:
         """
         Read data from a dfs2 file
 
@@ -253,9 +255,8 @@ class Dfs2(_Dfs123):
         n_items = len(item_numbers)
         items = _get_item_info(self._dfs.ItemInfo, item_numbers)
 
-        time_steps = _valid_timesteps(self._dfs.FileInfo, time)
-        nt = len(time_steps)
-        single_time_selected = np.isscalar(time) if time is not None else False
+        single_time_selected, time_steps = _valid_timesteps(self._dfs.FileInfo, time)
+        nt = len(time_steps) if not single_time_selected else 1
 
         if area is not None:
             take_subset = True
@@ -267,7 +268,7 @@ class Dfs2(_Dfs123):
             shape = (nt, self._ny, self._nx)
             geometry = self.geometry
 
-        if single_time_selected:
+        if single_time_selected and not keepdims:
             shape = shape[1:]
 
         data_list = [
@@ -288,7 +289,7 @@ class Dfs2(_Dfs123):
                 if take_subset:
                     d = np.take(np.take(d, jj, axis=0), ii, axis=-1)
 
-                if single_time_selected:
+                if single_time_selected and not keepdims:
                     data_list[item] = d
                 else:
                     data_list[item][i] = d
diff --git a/mikeio/dfs3.py b/mikeio/dfs3.py
index 9e3b9da2..b9ebaeb1 100644
--- a/mikeio/dfs3.py
+++ b/mikeio/dfs3.py
@@ -203,6 +203,7 @@ class Dfs3(_Dfs123):
         time_steps=None,
         area=None,
         layers=None,
+        keepdims=False,
     ) -> Dataset:
 
         if area is not None:
@@ -221,8 +222,8 @@ class Dfs3(_Dfs123):
                 )
             )
             time = time_steps
-        time_steps = _valid_timesteps(dfs.FileInfo, time)
-        nt = len(time_steps)
+        single_time_selected, time_steps = _valid_timesteps(dfs.FileInfo, time)
+        nt = len(time_steps) if not single_time_selected else 1
 
         # Determine the size of the grid
         zNum = self.geometry.nz
diff --git a/mikeio/dfsu.py b/mikeio/dfsu.py
index 16a3fef6..2259e70a 100644
--- a/mikeio/dfsu.py
+++ b/mikeio/dfsu.py
@@ -778,6 +778,7 @@ class _Dfsu(_UnstructuredFile, EquidistantTimeSeries):
         area=None,
         x=None,
         y=None,
+        keepdims=False,
     ) -> Dataset:
         """
         Read data from a dfsu file
@@ -842,8 +843,7 @@ class _Dfsu(_UnstructuredFile, EquidistantTimeSeries):
             )
             time = time_steps
 
-        single_time_selected = np.isscalar(time) if time is not None else False
-        time_steps = _valid_timesteps(dfs, time)
+        single_time_selected, time_steps = _valid_timesteps(dfs, time)
 
         self._validate_elements_and_geometry_sel(elements, area=area, x=x, y=y)
         if elements is None:
@@ -868,7 +868,11 @@ class _Dfsu(_UnstructuredFile, EquidistantTimeSeries):
         data_list = []
 
         n_steps = len(time_steps)
-        shape = (n_elems,) if single_time_selected else (n_steps, n_elems)
+        shape = (
+            (n_elems,)
+            if (single_time_selected and not keepdims)
+            else (n_steps, n_elems)
+        )
         for item in range(n_items):
             # Initialize an empty data block
             data = np.ndarray(shape=shape, dtype=self._dtype)
@@ -887,7 +891,7 @@ class _Dfsu(_UnstructuredFile, EquidistantTimeSeries):
                 if elements is not None:
                     d = d[elements]
 
-                if single_time_selected:
+                if single_time_selected and not keepdims:
                     data_list[item] = d
                 else:
                     data_list[item][i] = d
@@ -898,7 +902,10 @@ class _Dfsu(_UnstructuredFile, EquidistantTimeSeries):
 
         dfs.Close()
 
-        dims = ("time", "element") if not single_time_selected else ("element",)
+        dims = ("time", "element")
+
+        if single_time_selected and not keepdims:
+            dims = ("element",)
 
         if elements is not None and len(elements) == 1:
             # squeeze point data
@@ -1295,7 +1302,7 @@ class Dfsu2DH(_Dfsu):
         n_items = len(item_numbers)
 
         self._n_timesteps = dfs.NumberOfTimeSteps
-        time_steps = _valid_timesteps(dfs, time_steps=None)
+        _, time_steps = _valid_timesteps(dfs, time_steps=None)
 
         deletevalue = self.deletevalue
 
diff --git a/mikeio/dfsu_layered.py b/mikeio/dfsu_layered.py
index b72d4eef..63f140b1 100644
--- a/mikeio/dfsu_layered.py
+++ b/mikeio/dfsu_layered.py
@@ -105,6 +105,7 @@ class DfsuLayered(_Dfsu):
         y=None,
         z=None,
         layers=None,
+        keepdims=False,
     ) -> Dataset:
         """
         Read data from a dfsu file
@@ -171,8 +172,7 @@ class DfsuLayered(_Dfsu):
             )
             time = time_steps
 
-        single_time_selected = np.isscalar(time) if time is not None else False
-        time_steps = _valid_timesteps(dfs, time)
+        single_time_selected, time_steps = _valid_timesteps(dfs, time)
 
         self._validate_elements_and_geometry_sel(
             elements, area=area, layers=layers, x=x, y=y, z=z
@@ -221,7 +221,7 @@ class DfsuLayered(_Dfsu):
 
         t_seconds = np.zeros(n_steps, dtype=float)
 
-        if single_time_selected:
+        if single_time_selected and not keepdims:
             data = data[0]
 
         for i in trange(n_steps, disable=not self.show_progress):
@@ -238,7 +238,7 @@ class DfsuLayered(_Dfsu):
                     else:
                         d = d[elements]
 
-                if single_time_selected:
+                if single_time_selected and not keepdims:
                     data_list[item] = d
                 else:
                     data_list[item][i] = d
@@ -371,7 +371,7 @@ class Dfsu3D(DfsuLayered):
             elem3d = self.geometry.e2_e3_table[elem2d]
             return elem3d
 
-    def extract_surface_elevation_from_3d(self, filename=None, time=None, n_nearest=4):
+    def extract_surface_elevation_from_3d(self, filename=None, n_nearest=4):
         """
         Extract surface elevation from a 3d dfsu file (based on zn)
         to a new 2d dfsu file with a surface elevation item.
@@ -380,14 +380,12 @@ class Dfsu3D(DfsuLayered):
         ---------
         filename: str
             Output file name
-        time: str, int or list[int], optional
-            Extract only selected time_steps
         n_nearest: int, optional
             number of points for spatial interpolation (inverse_distance), default=4
 
         Examples
         --------
-        >>> dfsu.extract_surface_elevation_from_3d('ex_surf.dfsu', time='2018-1-1,2018-2-1')
+        >>> dfsu.extract_surface_elevation_from_3d('ex_surf.dfsu')
         """
         # validate input
         assert (
@@ -395,7 +393,6 @@ class Dfsu3D(DfsuLayered):
             or self._type == DfsuFileType.Dfsu3DSigmaZ
         )
         assert n_nearest > 0
-        time_steps = _valid_timesteps(self._source, time)
 
         # make 2d nodes-to-elements interpolator
         top_el = self.top_elements
@@ -410,11 +407,11 @@ class Dfsu3D(DfsuLayered):
             weights = get_idw_interpolant(dist)
 
         # read zn from 3d file and interpolate to element centers
-        ds = self.read(items=0, time_steps=time_steps)  # read only zn
+        ds = self.read(items=0, keepdims=True)  # read only zn
         node_ids_surf, _ = self.geometry._get_nodes_and_table_for_elements(
             top_el, node_layers="top"
         )
-        zn_surf = ds[0].values[:, node_ids_surf]  # surface
+        zn_surf = ds[0]._zn[:, node_ids_surf]  # surface
         surf2d = interp2d(zn_surf, node_ids, weights)
 
         # create output
diff --git a/mikeio/dfsu_spectral.py b/mikeio/dfsu_spectral.py
index acbce0be..b9656c91 100644
--- a/mikeio/dfsu_spectral.py
+++ b/mikeio/dfsu_spectral.py
@@ -80,6 +80,7 @@ class DfsuSpectral(_Dfsu):
         area=None,
         x=None,
         y=None,
+        keepdims=False,
     ) -> Dataset:
         """
         Read data from a spectral dfsu file
@@ -139,8 +140,7 @@ class DfsuSpectral(_Dfsu):
             )
             time = time_steps
 
-        single_time_selected = np.isscalar(time) if time is not None else False
-        time_steps = _valid_timesteps(dfs, time)
+        single_time_selected, time_steps = _valid_timesteps(dfs, time)
 
         if self._type == DfsuFileType.DfsuSpectral2D:
             self._validate_elements_and_geometry_sel(elements, area=area, x=x, y=y)
@@ -173,7 +173,7 @@ class DfsuSpectral(_Dfsu):
 
         t_seconds = np.zeros(n_steps, dtype=float)
 
-        if single_time_selected:
+        if single_time_selected and not keepdims:
             data = data[0]
 
         for i in trange(n_steps, disable=not self.show_progress):
@@ -191,7 +191,7 @@ class DfsuSpectral(_Dfsu):
                 if pts is not None:
                     d = d[pts, ...]
 
-                if single_time_selected:
+                if single_time_selected and not keepdims:
                     data_list[item] = d
                 else:
                     data_list[item][i] = d
diff --git a/mikeio/dfsutil.py b/mikeio/dfsutil.py
index 2c069dbe..2d8a8967 100644
--- a/mikeio/dfsutil.py
+++ b/mikeio/dfsutil.py
@@ -1,4 +1,6 @@
-from typing import Iterable, List, Union
+from datetime import datetime
+from multiprocessing.sharedctypes import Value
+from typing import Iterable, List, Tuple, Union
 import numpy as np
 import pandas as pd
 from .eum import EUMType, EUMUnit, ItemInfo, TimeAxisType
@@ -38,12 +40,15 @@ def _valid_item_numbers(
     return items
 
 
-def _valid_timesteps(dfsFileInfo: DfsFileInfo, time_steps):
-    # TODO: naming: time_steps or timesteps?
-    n_steps_file = dfsFileInfo.TimeAxis.NumberOfTimeSteps
+def _valid_timesteps(dfsFileInfo: DfsFileInfo, time_steps) -> Tuple[bool, List[int]]:
+
+    single_time_selected = False
+    if isinstance(time_steps, int) and np.isscalar(time_steps):
+        single_time_selected = True
 
+    n_steps_file = dfsFileInfo.TimeAxis.NumberOfTimeSteps
     if time_steps is None:
-        return list(range(n_steps_file))
+        return single_time_selected, list(range(n_steps_file))
 
     if isinstance(time_steps, int):
         time_steps = [time_steps]
@@ -60,29 +65,23 @@ def _valid_timesteps(dfsFileInfo: DfsFileInfo, time_steps):
         else:
             time_steps = slice(parts[0], parts[1])
 
-    if isinstance(time_steps, slice):
+    if isinstance(time_steps, (slice, pd.Timestamp, datetime, pd.DatetimeIndex)):
         if dfsFileInfo.TimeAxis.TimeAxisType != TimeAxisType.EquidistantCalendar:
             # TODO: handle non-equidistant calendar
             raise ValueError(
                 "Only equidistant calendar files are supported for this type of time_step argument"
             )
+
         start_time_file = dfsFileInfo.TimeAxis.StartDateTime
         time_step_file = dfsFileInfo.TimeAxis.TimeStep
-
-        freq = pd.tseries.offsets.DateOffset(seconds=time_step_file)
+        freq = pd.Timedelta(seconds=time_step_file)
         time = pd.date_range(start_time_file, periods=n_steps_file, freq=freq)
-        if time_steps.start is None:
-            time_steps_start = time[0]
-        else:
-            time_steps_start = pd.Timestamp(time_steps.start)
-        if time_steps.stop is None:
-            time_steps_stop = time[-1]
-        else:
-            time_steps_stop = pd.Timestamp(time_steps.stop)
 
-        s = time.slice_indexer(time_steps_start, time_steps_stop)
+    if isinstance(time_steps, slice):
+
+        s = time.slice_indexer(time_steps.start, time_steps.stop)
         time_steps = list(range(s.start, s.stop))
-    elif isinstance(time_steps[0], int):
+    elif isinstance(time_steps, Iterable) and isinstance(time_steps[0], int):
         time_steps = np.array(time_steps)
         time_steps[time_steps < 0] = n_steps_file + time_steps[time_steps < 0]
         time_steps = list(time_steps)
@@ -91,8 +90,25 @@ def _valid_timesteps(dfsFileInfo: DfsFileInfo, time_steps):
             raise IndexError(f"Timestep cannot be larger than {n_steps_file}")
         if min(time_steps) < 0:
             raise IndexError(f"Timestep cannot be less than {-n_steps_file}")
-
-    return time_steps
+    elif isinstance(time_steps, Iterable):
+        steps = []
+        for t in time_steps:
+            _, step = _valid_timesteps(dfsFileInfo, t)
+            steps.append(step[0])
+        single_time_selected = len(steps) == 1
+        time_steps = steps
+
+    elif isinstance(time_steps, (pd.Timestamp, datetime)):
+        s = time.slice_indexer(time_steps, time_steps)
+        time_steps = list(range(s.start, s.stop))
+    elif isinstance(time_steps, pd.DatetimeIndex):
+        time_steps = list(time.get_indexer(time_steps))
+
+    else:
+        raise TypeError(f"Indexing is not possible with {type(time_steps)}")
+    if len(time_steps) == 1:
+        single_time_selected = True
+    return single_time_selected, time_steps
 
 
 def _item_numbers_by_name(dfsItemInfo, item_names, ignore_first=False):


Test patch:
diff --git a/tests/test_consistency.py b/tests/test_consistency.py
index 4724a82c..1bcbeedc 100644
--- a/tests/test_consistency.py
+++ b/tests/test_consistency.py
@@ -1,6 +1,7 @@
+import pandas as pd
+from datetime import datetime
 import pytest
 
-import numpy as np
 import mikeio
 from mikeio.dataarray import DataArray
 from mikeio.spatial.geometry import GeometryUndefined
@@ -137,8 +138,7 @@ def test_read_dfs2_single_time():
     assert "time" not in ds.dims
 
     ds = mikeio.read(
-        "tests/testdata/consistency/oresundHD.dfs2",
-        time=[-1],  # time as array, forces time dimension to be kept
+        "tests/testdata/consistency/oresundHD.dfs2", time=[-1], keepdims=True
     )
 
     assert ds.n_timesteps == 1
@@ -207,7 +207,246 @@ def test_read_dfsu2d_single_time():
     ds = mikeio.read(
         "tests/testdata/consistency/oresundHD.dfsu",
         time=[-1],
+        keepdims=True,
     )
 
     assert ds.n_timesteps == 1
     assert "time" in ds.dims
+
+
+def test_read_dfs_time_selection_str():
+
+    extensions = ["dfsu", "dfs2", "dfs1", "dfs0"]
+    for ext in extensions:
+        filename = f"tests/testdata/consistency/oresundHD.{ext}"
+        time = "2018-03"
+        ds = mikeio.read(filename=filename)
+        dssel = ds.sel(time=time)
+        assert dssel.n_timesteps == 5
+
+        dsr = mikeio.read(filename=filename, time=time)
+        assert all(dsr.time == dssel.time)
+        assert dsr.shape == dssel.shape
+
+        dsgetitem = ds[time]
+        assert all(dsr.time == dsgetitem.time)
+        assert dsr.shape == dsgetitem.shape
+
+
+def test_read_dfs_time_selection_str_specific():
+
+    extensions = ["dfsu", "dfs2", "dfs1", "dfs0"]
+    for ext in extensions:
+        filename = f"tests/testdata/consistency/oresundHD.{ext}"
+        time = "2018-03-08 00:00:00"
+        ds = mikeio.read(filename=filename)
+        dssel = ds.sel(time=time)
+        assert dssel.n_timesteps == 1
+        assert "time" not in dssel.dims
+
+        dsr = mikeio.read(filename=filename, time=time)
+        assert all(dsr.time == dssel.time)
+        assert dsr.shape == dssel.shape
+
+        dsgetitem = ds[time]
+        assert all(dsr.time == dsgetitem.time)
+        assert dsr.shape == dsgetitem.shape
+
+
+def test_read_dfs_time_selection_list_str():
+
+    extensions = ["dfsu", "dfs2", "dfs1", "dfs0"]
+    for ext in extensions:
+        filename = f"tests/testdata/consistency/oresundHD.{ext}"
+        time = ["2018-03-08 00:00", "2018-03-10 00:00"]
+        ds = mikeio.read(filename=filename)
+        dssel = ds.sel(time=time)
+        assert dssel.n_timesteps == 2
+
+        dsr = mikeio.read(filename=filename, time=time)
+        assert all(dsr.time == dssel.time)
+        assert dsr.shape == dssel.shape
+
+        dsgetitem = ds[time]
+        assert all(dsr.time == dsgetitem.time)
+        assert dsr.shape == dsgetitem.shape
+
+
+def test_read_dfs_time_selection_pdTimestamp():
+
+    extensions = ["dfsu", "dfs2", "dfs1", "dfs0"]
+    for ext in extensions:
+        filename = f"tests/testdata/consistency/oresundHD.{ext}"
+        time = pd.Timestamp("2018-03-08 00:00:00")
+        ds = mikeio.read(filename=filename)
+        dssel = ds.sel(time=time)
+        assert dssel.n_timesteps == 1
+        assert "time" not in dssel.dims
+
+        dsr = mikeio.read(filename=filename, time=time)
+        assert all(dsr.time == dssel.time)
+        assert dsr.shape == dssel.shape
+
+        dsgetitem = ds[time]
+        assert all(dsr.time == dsgetitem.time)
+        assert dsr.shape == dsgetitem.shape
+
+
+def test_read_dfs_time_selection_pdDatetimeIndex():
+
+    extensions = ["dfsu", "dfs2", "dfs1", "dfs0"]
+    for ext in extensions:
+        filename = f"tests/testdata/consistency/oresundHD.{ext}"
+        time = pd.date_range("2018-03-08", end="2018-03-10", freq="D")
+        ds = mikeio.read(filename=filename)
+        dssel = ds.sel(time=time)
+        assert dssel.n_timesteps == 3
+
+        dsr = mikeio.read(filename=filename, time=time)
+        assert all(dsr.time == dssel.time)
+        assert dsr.shape == dssel.shape
+
+        dsgetitem = ds[time]
+        assert all(dsr.time == dsgetitem.time)
+        assert dsr.shape == dsgetitem.shape
+
+
+def test_read_dfs_time_selection_datetime():
+
+    extensions = ["dfsu", "dfs2", "dfs1", "dfs0"]
+    for ext in extensions:
+        filename = f"tests/testdata/consistency/oresundHD.{ext}"
+        time = datetime(2018, 3, 8, 0, 0, 0)
+        ds = mikeio.read(filename=filename)
+        dssel = ds.sel(time=time)
+        assert dssel.n_timesteps == 1
+        assert "time" not in dssel.dims
+
+        dsr = mikeio.read(filename=filename, time=time)
+        assert all(dsr.time == dssel.time)
+        assert dsr.shape == dssel.shape
+
+        dsgetitem = ds[time]
+        assert all(dsr.time == dsgetitem.time)
+        assert dsr.shape == dsgetitem.shape
+
+        dsr2 = mikeio.read(filename=filename, time=pd.Timestamp(time))
+        assert all(dsr2.time == dsr.time)
+        assert dsr2.shape == dsr.shape
+
+
+def test_read_dfs_time_list_datetime():
+
+    extensions = ["dfsu", "dfs2", "dfs1", "dfs0"]
+    for ext in extensions:
+        filename = f"tests/testdata/consistency/oresundHD.{ext}"
+        time = [datetime(2018, 3, 8), datetime(2018, 3, 10)]
+        ds = mikeio.read(filename=filename)
+        dssel = ds.sel(time=time)
+        assert dssel.n_timesteps == 2
+
+        dsr = mikeio.read(filename=filename, time=time)
+        assert all(dsr.time == dssel.time)
+        assert dsr.shape == dssel.shape
+
+        dsgetitem = ds[time]
+        assert all(dsr.time == dsgetitem.time)
+        assert dsr.shape == dsgetitem.shape
+
+
+def test_read_dfs_time_slice_datetime():
+
+    extensions = ["dfsu", "dfs2", "dfs1", "dfs0"]
+    for ext in extensions:
+        filename = f"tests/testdata/consistency/oresundHD.{ext}"
+        time = slice(datetime(2018, 3, 8), datetime(2018, 3, 10))
+        ds = mikeio.read(filename=filename)
+        dssel = ds.sel(time=time)
+        assert dssel.n_timesteps == 3
+
+        dsr = mikeio.read(filename=filename, time=time)
+        assert all(dsr.time == dssel.time)
+        assert dsr.shape == dssel.shape
+
+        dsgetitem = ds[time]
+        assert all(dsr.time == dsgetitem.time)
+        assert dsr.shape == dsgetitem.shape
+
+
+def test_read_dfs_time_slice_str():
+
+    extensions = ["dfsu", "dfs2", "dfs1", "dfs0"]
+    for ext in extensions:
+        filename = f"tests/testdata/consistency/oresundHD.{ext}"
+        time = slice("2018-03-08", "2018-03-10")
+        ds = mikeio.read(filename=filename)
+        dssel = ds.sel(time=time)
+        assert dssel.n_timesteps == 3
+
+        dsr = mikeio.read(filename=filename, time=time)
+        assert all(dsr.time == dssel.time)
+        assert dsr.shape == dssel.shape
+
+        dsgetitem = ds[time]
+        assert all(dsr.time == dsgetitem.time)
+        assert dsr.shape == dsgetitem.shape
+
+
+def test_read_dfs_time_selection_str_comma():
+
+    extensions = ["dfs0", "dfs2", "dfs1", "dfs0"]
+    for ext in extensions:
+        filename = f"tests/testdata/consistency/oresundHD.{ext}"
+        time = "2018-03-08,2018-03-10"
+        ds = mikeio.read(filename=filename)
+        dssel = ds.sel(time=time)
+        assert dssel.n_timesteps == 3
+
+        dsr = mikeio.read(filename=filename, time=time)
+        assert all(dsr.time == dssel.time)
+        assert dsr.shape == dssel.shape
+
+        dsgetitem = ds[time]
+        assert all(dsr.time == dsgetitem.time)
+        assert dsr.shape == dsgetitem.shape
+
+
+def test_read_dfs_time_int():
+
+    extensions = ["dfsu", "dfs2", "dfs1", "dfs0"]
+    for ext in extensions:
+        filename = f"tests/testdata/consistency/oresundHD.{ext}"
+        time = -1
+        ds = mikeio.read(filename=filename)
+        dssel = ds.isel(time=time)
+        assert dssel.n_timesteps == 1
+        assert "time" not in dssel.dims
+
+        dsr = mikeio.read(filename=filename, time=time)
+        assert all(dsr.time == dssel.time)
+        assert dsr.shape == dssel.shape
+
+        # integer time selection for DataArray (not Dataset)
+        dsgetitem = ds[0][time]
+        assert all(dsr[0].time == dsgetitem.time)
+        assert dsr[0].shape == dsgetitem.shape
+
+
+def test_read_dfs_time_list_int():
+
+    extensions = ["dfsu", "dfs2", "dfs1", "dfs0"]
+    for ext in extensions:
+        filename = f"tests/testdata/consistency/oresundHD.{ext}"
+        time = [0, 1]
+        ds = mikeio.read(filename=filename)
+        dssel = ds.isel(time=time)
+        assert dssel.n_timesteps == 2
+
+        dsr = mikeio.read(filename=filename, time=time)
+        assert all(dsr.time == dssel.time)
+        assert dsr.shape == dssel.shape
+
+        # integer time selection for DataArray (not Dataset)
+        dsgetitem = ds[0][time]
+        assert all(dsr[0].time == dsgetitem.time)
+        assert dsr[0].shape == dsgetitem.shape
diff --git a/tests/test_dfs2.py b/tests/test_dfs2.py
index 00b3a68a..53ab4c6b 100644
--- a/tests/test_dfs2.py
+++ b/tests/test_dfs2.py
@@ -660,13 +660,13 @@ def test_incremental_write_from_dfs2(tmpdir):
 
     nt = dfs.n_timesteps
 
-    ds = dfs.read(time=[0])
+    ds = dfs.read(time=[0], keepdims=True)
 
     dfs_to_write = Dfs2()
     dfs_to_write.write(outfilename, ds, dt=dfs.timestep, keep_open=True)
 
     for i in range(1, nt):
-        ds = dfs.read(time=[i])
+        ds = dfs.read(time=[i], keepdims=True)
         dfs_to_write.append(ds)
 
     dfs_to_write.close()
@@ -686,13 +686,13 @@ def test_incremental_write_from_dfs2_context_manager(tmpdir):
 
     nt = dfs.n_timesteps
 
-    ds = dfs.read(time=[0])
+    ds = dfs.read(time=[0], keepdims=True)
 
     dfs_to_write = Dfs2()
     with dfs_to_write.write(outfilename, ds, dt=dfs.timestep, keep_open=True) as f:
 
         for i in range(1, nt):
-            ds = dfs.read(time=[i])
+            ds = dfs.read(time=[i], keepdims=True)
             f.append(ds)
 
         # dfs_to_write.close() # called automagically by context manager
@@ -708,7 +708,7 @@ def test_read_concat_write_dfs2(tmp_path):
 
     ds1 = mikeio.read("tests/testdata/waves.dfs2", time=[0, 1])
     # ds2 = mikeio.read("tests/testdata/waves.dfs2", time=2)  # dont do this, it will not work, since reading a single time step removes the time dimension
-    ds2 = mikeio.read("tests/testdata/waves.dfs2", time=[2])
+    ds2 = mikeio.read("tests/testdata/waves.dfs2", time=[2], keepdims=True)
     dsc = mikeio.Dataset.concat([ds1, ds2])
     assert dsc.n_timesteps == 3
     assert dsc.end_time == ds2.end_time
diff --git a/tests/test_dfsu.py b/tests/test_dfsu.py
index 634f5361..6612ea11 100644
--- a/tests/test_dfsu.py
+++ b/tests/test_dfsu.py
@@ -191,7 +191,7 @@ def test_read_single_time_step():
     ds = dfs.read(items=[0, 3], time=1)
     assert "time" not in ds.dims
 
-    ds = dfs.read(items=[0, 3], time=[1])  # this forces time dimension to be kept
+    ds = dfs.read(items=[0, 3], time=[1], keepdims=True)
     assert "time" in ds.dims
 
 
@@ -471,12 +471,12 @@ def test_incremental_write_from_dfsu(tmpdir):
 
     nt = dfs.n_timesteps
 
-    ds = dfs.read(time=[0])
+    ds = dfs.read(time=[0], keepdims=True)
 
     dfs.write(outfilename, ds, keep_open=True)
 
     for i in range(1, nt):
-        ds = dfs.read(time=[i])
+        ds = dfs.read(time=[i], keepdims=True)
         dfs.append(ds)
 
     dfs.close()
@@ -495,11 +495,11 @@ def test_incremental_write_from_dfsu_context_manager(tmpdir):
 
     nt = dfs.n_timesteps
 
-    ds = dfs.read(time=[0])
+    ds = dfs.read(time=[0], keepdims=True)
 
     with dfs.write(outfilename, ds, keep_open=True) as f:
         for i in range(1, nt):
-            ds = dfs.read(time=[i])
+            ds = dfs.read(time=[i], keepdims=True)
             f.append(ds)
 
         # dfs.close() # should be called automagically by context manager
diff --git a/tests/test_dfsu_layered.py b/tests/test_dfsu_layered.py
index 85a2639d..ce67bc77 100644
--- a/tests/test_dfsu_layered.py
+++ b/tests/test_dfsu_layered.py
@@ -494,7 +494,7 @@ def test_extract_surface_elevation_from_3d():
     outputfile = "tests/testdata/oresund_surface_elev_extracted.dfsu"
     n_top1 = len(dfs.top_elements)
 
-    dfs.extract_surface_elevation_from_3d(outputfile, time=-1)
+    dfs.extract_surface_elevation_from_3d(outputfile)
 
     dfs2 = Dfsu(outputfile)
     assert dfs2.n_elements == n_top1